这篇博文有些不同寻常,它是面向未来的:它主要关注尚不存在的事物。其目的是为社区关于AbstractArray
核心 API 可能更改的讨论奠定基础,并作为更专注的“julep”(Julia 增强提案)的背景阅读和参考材料。这里,我通常使用“数组”作为AbstractArray
的简写,如果我明确指的是 Julia 的具体Array
类型,则使用Array
。
读者可能已经知道,在 Julia 中,可以编写算法,其中一个或多个输入仅假定为AbstractArray
。这是“泛型”代码,意味着它应该适用于任何特定的具体数组类型(即产生正确的结果)。在一个理想的世界中——Julia 在很多情况下都很好地接近了这个理想——代码的通用性不应该对其性能产生负面影响:泛型实现应该与仅限于特定数组类型的实现一样快。这意味着泛型算法应该使用低级操作来编写,这些操作在各种数组类型中都能提供良好的性能。
提供高效的低级操作是一项与“向量化”所有内容的编程语言所面临的设计挑战不同的挑战。如果成功,它将极大地促进代码的重用,因为高效、泛型的低级部分允许您编写各种高效、泛型的更高级函数。
自然地,随着数组类型的多样性增加,我们必须更加谨慎地对待这些低级操作的抽象。
在讨论数组上的通用操作时,心中有一个多样的具体数组集合是有用的。
在 Julia 核心库中,我们相当好地支持的一些类型是
Array
:所有数组的原型
Range
:一个很好的“计算”数组的例子,其中实际上没有任何值存储在内存中。由于没有存储,因此这些是不可变的容器:您不能设置各个槽中的值。
BitArray
:只能存储 0 或 1(false
或 true
)的数组,并且其内部存储被压缩,以便每个条目仅需要一个比特。
SubArray
:这种类型引入的问题以及我们采用的解决方案可能是此处所考虑的泛化的最佳模型。因此,此案例将在下面更详细地讨论。
Base 中另一类重要的数组类型是稀疏数组:SparseMatrixCSC
和 SparseVector
,以及其他稀疏表示,如 Diagonal
、Bidiagonal
和 Tridiagonal
。这些是数组类型的很好的例子,其中访问模式值得仔细考虑。值得注意的是,尽管我们拥有的 5 种左右的稀疏参数化在“策略”方面有很多共同点,但核心算法(例如矩阵乘法)的实现针对每种稀疏类型都是专门的——换句话说,这些模仿了其他语言中常见的“高级向量化函数”策略。我们缺少的是一个“稀疏迭代 API”,它允许您以通用的方式高效地编写稀疏线性代数的主要算法。我们目前的模型对于SparseLike*Dense
操作可能还不错,但如果您想高效地计算例如Bidiagonal*SparseMatrixCSC
,则变得难以管理:您需要支持的可能组合的数量随着稀疏类型的增加而迅速增长,因此代表了开发高效、泛型低级操作的强大动力。
在 Base 之外,还有一些其他令人脑洞大开的数组示例,包括
DataFrames
:使用符号而不是整数来索引数组。其他相关类型包括NamedArrays
、AxisArrays
。
Interpolations
:使用非整数浮点数来索引数组
DistributedArrays
:另一个需要仔细考虑访问模式的极佳示例
对于固定维度的数组,可以编写算法来索引数组,如A[i,j,k,...]
(可以在我们的线性代数代码中找到很好的例子,其中所有内容都是向量或矩阵)。对于必须支持任意维度的算法,长期以来我们的后备方案是线性索引,对于整数i
使用A[i]
。但是,一般来说,SubArrays 无法通过线性索引高效地访问,因为这会导致对div
的调用,而div
很慢。这是一个 CPU 问题,而不是 Julia 特定的问题。尽管最近添加了基础设施以使其速度快得多,但div
的缓慢仍然是事实——现在可以使其仅仅是“非常糟糕”而不是"糟糕透顶"。
我们(很大程度上)解决此问题的方法是,使对任意维度的数组进行笛卡尔索引A[i,j,k,...]
成为可能(CartesianIndex
类型)。为了在实际代码中利用这一点,我们还必须使用for I in eachindex(A)
结构扩展我们的迭代器。这允许您选择一个迭代器,该迭代器优化对A
元素的访问效率。在泛型算法中,性能提升并不小,有时可达十到五十倍。这些类型在之前的博文中进行了描述。
据我所知,这种方法使 Julia 拥有了任何编程语言中最灵活且高效的“数组视图”类型之一。许多语言将视图基于数组的步长,这意味着内存偏移量沿每个维度都是规则的。除此之外,这还要求底层数组是稠密的。相反,在 Julia 中,我们可以轻松处理非步长数组(例如,沿一个维度在[1,3,17,428,...]
处采样,或创建SparseMatrixCSC
的视图)。我们还可以处理没有底层存储的数组(例如Range
)。能够使用通用基础设施做到这一点是使不同的优化数组类型在泛型编程中发挥作用的部分原因。
也值得指出一些问题
最重要的是,它要求采用略微不同的编程风格。尽管已经进入另一个发布周期,但这种转变甚至在 Base 中都尚未完成。
对于涉及两个或多个数组的算法,它们的“最佳”迭代器可能属于不同的类型。原则上,这是一个大问题。考虑矩阵向量乘法A[i,j]*v[j]
,其中j
需要对A
和v
保持同步,但您也希望所有这些访问都尽可能高效。实际上,现在这不是一个紧迫的问题:即使我们的数组都不具有高效的线性索引,据我所知,我们所有的(稠密)数组类型都具有高效的笛卡尔索引。由于使用N
个整数(其中N
等于数组的维数)进行索引始终具有良好的性能,因此这成为泛型代码的可靠默认值。(值得注意的是,对于稀疏数组,情况并非如此,并且缺少相应的泛型解决方案可能是我们缺少编写稀疏算法的泛型 API 的主要原因。)
不幸的是,我怀疑如果我们想添加对某些新操作或类型的支持(下面有具体示例),它将迫使我们点燃后一个问题。
一些可能的新的AbstractArray
类型带来了新的挑战。
这些是这篇文章的中心动机。这些动机来自确保reshape(A, dims)
始终返回A
的“视图”而不是分配A
的副本。(如果我们决定放弃这个目标,那么这个 julep 的大部分紧迫性都会消失,在这种情况下,为了保持一致性,我们应该始终返回A
的副本。)值得注意的是,除了显式的reshape
之外,我们还有一些机制会导致创建副本的重塑,特别是应用于 3D 数组的A[:]
或A[:, :]
。
与SubArrays
类似,ReshapedArrays
的主要挑战是获得良好的性能。如果A
是一个 3D 数组,并且您将其重塑为一个 2D 数组B
,那么B[i,j]
必须扩展为A[k,l,m]
。问题在于计算正确的k,l,m
可能会导致对div
的调用。因此,ReshapedArrays 违反了我们当前生态系统的支柱,即使用N
个整数进行索引可能不是访问B
元素的最快方法。从性能的角度来看,这个问题很大(参见#15449,大约五到十倍)。
在简单的情况下,有一种简单的方法可以规避此性能问题:定义一种新的迭代器类型,该类型(在内部)直接迭代父A
的索引。换句话说,创建一个迭代器,以便B[I]
立即扩展为A[I']
,并且后者具有“理想”的性能。
不幸的是,当您需要使两个数组保持同步时,此策略会遇到很多麻烦:如果您想采用此策略,则根本无法再为矩阵向量乘法编写B[i,j]*v[j]
。解决此问题的潜在方法是定义一类新的迭代器,这些迭代器在数组的特定维度上操作(#15459),编写B[ii,jj]*v[j]
。jj
(无论是什么)和j
需要保持同步,但它们不一定都需要是整数。使用这种策略,矩阵向量乘法
for j = 1:size(B, 2)
vj = v[j]
for i = 1:size(B, 1)
dest[i] += B[i,j] * vj
end
end
可以以更有效的方式编写,如下所示
for (jj, vj) in zip(eachindex(B, Dimension{2}), v)
for (i, ii) in zip(eachindex(dest), eachindex(B, (:, jj)))
dest[i] += B[ii,jj]*vj
end
end
弄清楚eachindex(B, Dimension{2})
和eachindex(B, (:, jj))
应该做什么并不难:例如,ii
可以是一个CartesianInnerIndex
(一种尚不存在的类型),对于B
的特定列,它从A[3,7,4]
迭代到A[5,8,4]
,其中第d
个索引分量在size(A, d)
处环绕。此策略的最大性能优势在于,您只需要计算一个div
来设置每列迭代器的边界;内部循环不需要在每次元素访问时都计算div
。毫无疑问,如果合适地定义jj
,则可以更加巧妙地完全避免计算div
。对作者而言,此策略似乎有望解决有关 ReshapedArrays 的大多数性能问题——只有在您需要“随机访问”时,您才需要缓慢的(基于整数的)操作。
但是,一个大问题是,与“朴素”实现相比,这相当丑陋。
Julia 的 Array
类型以列优先顺序存储其元素,这意味着 A[i,j]
和 A[i+1,j]
位于相邻的内存位置。对于某些应用——或与某些外部代码库交互——支持行优先数组可能很方便,在这种情况下,A[i,j]
和 A[i,j+1]
位于相邻的内存位置。从更根本的角度来说,这与 Julia 开发历史上讨论最多的问题之一有关,称为“认真对待转置”或 #4774。至少已经尝试过两种实现方式,#6837 和 mb/transpose
分支,对于后者,其优势和挑战的总结已 发布。
提到的最大挑战之一是需要支持的方法数量急剧增加。泛型代码可以在这里提供帮助吗?有两个相关的考虑因素。第一个是线性索引:通常将其与“存储顺序”混淆,即,对于同一数组的两个线性索引 i
和 j
,内存中的偏移量与 i-j
成正比。对于行优先数组,这个概念不可行,因为否则循环
function copy!(dest, src)
for i = 1:length(src)
dest[i] = src[i] # trouble if `i` means "memory offset"
end
dest
end
如果 src
和 dest
使用不同的存储顺序,最终将进行转置。因此,必须根据相应的笛卡尔(全维度)索引来定义线性索引。这不是什么大问题,因为我们知道如何解决:在必要时使用 ind2sub
(速度很慢),但为了效率,使行优先数组属于默认使用笛卡尔索引进行迭代的数组类别 (LinearSlow
)。这样做将确保如果使用像 eachindex(src)
这样的泛型构造而不是 1:length(src)
,则上述循环可以很快。
更具挑战性的问题涉及缓存效率:以与 存储顺序 不同的任何方式访问数组元素的速度都会慢得多。一些编写矩阵向量乘法的合理快速方法是
for j = 1:size(B, 2)
vj = v[j]
for i = 1:size(B, 1)
dest[i] += B[i,j] * vj
end
end
对于列优先矩阵 B
,以及
for i = 1:size(B, 1)
for j = 1:size(B, 2)
dest[i] += B[i,j] * v[j]
end
end
对于行优先矩阵。(通过使用标量临时累加器可以做得更好,但这里我们不考虑这一点。)需要注意的关键点是循环的顺序已切换。
可以通过定义一个 RowMajorRange
迭代器来概括这一点,它非常类似于我们的 CartesianRange
迭代器,但以行优先顺序遍历数组。eachindex
声称返回一个“高效迭代器”,毫无疑问,对于行优先数组,RowMajorRange
比 CartesianRange
迭代器更有效。所以让我们假设 eachindex
做到了它所说的,并返回一个 RowMajorRange
迭代器。使用此策略,上述两种算法可以合并到一个泛型实现中
for I in eachindex(B)
dest[I[1]] += B[I]*v[I[2]]
end
耶!高效泛型实现得分。
但我们的胜利是短暂的。让我们回到上面 copy!
的示例,并意识到 dest
和 src
可能属于两种不同的数组类型,因此可能最有效地使用不同的迭代器类型进行索引。我们很想这样写:
function copy!(dest, src)
for (idest, isrc) in zip(eachindex(dest), eachindex(src))
dest[idest] = src[isrc]
end
dest
end
在我们为 eachindex
引入 RowMajorRange
返回类型之前,此实现是可以的。但我们刚刚破坏了它,因为现在在某些情况下,这将错误地进行转置。
换句话说,如果没有仔细的设计,“最大效率迭代”和“保持访问同步”的目标就会发生冲突。
Julia 的数组从 1 开始索引,而其他一些语言从 0 开始编号。如果你根据各种博客文章的评论来判断,那么有大量的程序员迫切希望采用 Julia,但由于索引的这种差异,他们甚至不会尝试。由于招募这些军队将导致世界统治,这显然是一个极其紧迫的问题。
更严肃地说,确实存在一些算法,如果可以索引超出 1:size(A,d)
的范围,则这些算法会得到简化。在我们自己的实验室内部代码中,我们长期以来一直在使用 CenterIndexedArray
类型,其中此类数组(所有数组的大小都是奇数)在 -n:n
范围内进行索引,其中 0 指的是“中心”元素。OffsetArrays
包推广了这个概念。不幸的是,在实践中,对于 Julia 的许多核心函数,这两种数组类型都会产生段错误(由于对 @inbounds
何时合适做出了内置假设);随着时间的推移,我们的实验室不得不为相当多的 Julia 函数编写专门针对 CenterIndexedArrays
的实现。
OffsetArrays
说明了另一个概念上的挑战,可以通过 copy!
轻松地演示出来。当 dest
是一个一维 OffsetArray
且 src
是一个标准 Vector
时,copy!
应该做什么?特别是,src[1]
放在哪里?它放在 dest
的第一个元素中,还是存储在 dest[1]
中(这可能不是第一个元素)。
此类示例迫使我们更深入地思考数组的真正含义。似乎有两个潜在的概念。一个是数组是列表,多维数组是列表的列表的列表的……在这种世界观中,正确的方法是将 src[1]
放入 dest
的第一个槽中,因为 1 只是 first
的同义词。但是,这种世界观并没有真正赋予数组的索引元组任何“意义”,从某种意义上说,它甚至没有包含 OffsetArray
传达的区别。换句话说,在这个世界中,OffsetArray
根本没有意义,不应该存在。
如果相反地认为 OffsetArray
应该存在,这实质上迫使人们采用不同的世界观:数组实际上是关联容器,其中每个索引元组都是检索值的“键”。在这种思维方式下,src[1]
应该存储在 dest[1]
中。
这些示例建议对 AbstractArray
进行形式化
AbstractArray
是专门的关联容器,因为允许的“键”可能不仅受其 Julia 类型限制。具体来说,允许的键必须能够表示为一维值列表的笛卡尔积。允许的键可能不仅取决于数组类型,还取决于特定的数组(例如,其大小)。尝试使用无法转换为特定数组的允许键之一的键进行访问将导致 BoundsError
。
对于任何给定的数组,必须能够从数组本身生成有效键的完整域的有限维参数化。这可能只需要知道数组大小,或者键可能依赖于某些内部存储(例如 DataFrames
和 OffsetArrays
)。在某些情况下,仅数组类型可能就足够了(例如,FixedSizeArrays
)。根据此定义,请注意,Dict{ASCII5,Int}
(其中 ASCII5
是一种表示具有 5 个字符的 ASCII 字符串的类型)将被视为一个 5 维(稀疏)数组,但 Dict{ASCIIString,Int}
则不会(因为 ASCIIString
没有长度限制,因此没有有限维度)。
数组可以通过多种键类型进行索引(即,键可能有多种参数化)。当不同的键参数化引用给定数组的相同元素时,它们是等价的。线性索引和笛卡尔索引是可互换表示的简单示例,但专门的迭代器也可以生成其他键类型。
数组可能支持多个产生非等价键序列的迭代器。换句话说,行优先矩阵可以支持 CartesianRange
和 RowMajorRange
迭代器,它们以不同的顺序访问元素。
解决这些相互冲突的需求并非易事。一种方法可能是规定某些数组类型根本无法使用泛型代码支持。这可能是正确的策略。或者,可以尝试设计一个处理所有这些类型(并希望更多类型)的数组 API。
在 GitHub 问题 #15648 中,我们正在讨论可能解决这些挑战的 API。鼓励读者参与此讨论。