抽象数组的泛化:机遇与挑战

2016年3月27日 | Tim Holy

引言:使用抽象数组的通用算法

这篇博文有些不同寻常,它是面向未来的:它主要关注尚不存在的事物。其目的是为社区关于AbstractArray核心 API 可能更改的讨论奠定基础,并作为更专注的“julep”(Julia 增强提案)的背景阅读和参考材料。这里,我通常使用“数组”作为AbstractArray的简写,如果我明确指的是 Julia 的具体Array类型,则使用Array

读者可能已经知道,在 Julia 中,可以编写算法,其中一个或多个输入仅假定为AbstractArray。这是“泛型”代码,意味着它应该适用于任何特定的具体数组类型(即产生正确的结果)。在一个理想的世界中——Julia 在很多情况下都很好地接近了这个理想——代码的通用性不应该对其性能产生负面影响:泛型实现应该与仅限于特定数组类型的实现一样快。这意味着泛型算法应该使用低级操作来编写,这些操作在各种数组类型中都能提供良好的性能。

提供高效的低级操作是一项与“向量化”所有内容的编程语言所面临的设计挑战不同的挑战。如果成功,它将极大地促进代码的重用,因为高效、泛型的低级部分允许您编写各种高效、泛型的更高级函数。

自然地,随着数组类型的多样性增加,我们必须更加谨慎地对待这些低级操作的抽象。

数组示例

在讨论数组上的通用操作时,心中有一个多样的具体数组集合是有用的。

在 Julia 核心库中,我们相当好地支持的一些类型是

Base 中另一类重要的数组类型是稀疏数组:SparseMatrixCSCSparseVector,以及其他稀疏表示,如 DiagonalBidiagonalTridiagonal。这些是数组类型的很好的例子,其中访问模式值得仔细考虑。值得注意的是,尽管我们拥有的 5 种左右的稀疏参数化在“策略”方面有很多共同点,但核心算法(例如矩阵乘法)的实现针对每种稀疏类型都是专门的——换句话说,这些模仿了其他语言中常见的“高级向量化函数”策略。我们缺少的是一个“稀疏迭代 API”,它允许您以通用的方式高效地编写稀疏线性代数的主要算法。我们目前的模型对于SparseLike*Dense操作可能还不错,但如果您想高效地计算例如Bidiagonal*SparseMatrixCSC,则变得难以管理:您需要支持的可能组合的数量随着稀疏类型的增加而迅速增长,因此代表了开发高效、泛型低级操作的强大动力。

在 Base 之外,还有一些其他令人脑洞大开的数组示例,包括

SubArrays:一个案例研究

对于固定维度的数组,可以编写算法来索引数组,如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)。能够使用通用基础设施做到这一点是使不同的优化数组类型在泛型编程中发挥作用的部分原因。

也值得指出一些问题

不幸的是,我怀疑如果我们想添加对某些新操作或类型的支持(下面有具体示例),它将迫使我们点燃后一个问题。

具有挑战性的例子

一些可能的新的AbstractArray类型带来了新的挑战。

ReshapedArrays(#15449

这些是这篇文章的中心动机。这些动机来自确保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 的大多数性能问题——只有在您需要“随机访问”时,您才需要缓慢的(基于整数的)操作。

但是,一个大问题是,与“朴素”实现相比,这相当丑陋。

行主序矩阵、PermutedDimensionArrays 和“认真对待转置”

Julia 的 Array 类型以列优先顺序存储其元素,这意味着 A[i,j]A[i+1,j] 位于相邻的内存位置。对于某些应用——或与某些外部代码库交互——支持行优先数组可能很方便,在这种情况下,A[i,j]A[i,j+1] 位于相邻的内存位置。从更根本的角度来说,这与 Julia 开发历史上讨论最多的问题之一有关,称为“认真对待转置”或 #4774。至少已经尝试过两种实现方式,#6837mb/transpose 分支,对于后者,其优势和挑战的总结已 发布

提到的最大挑战之一是需要支持的方法数量急剧增加。泛型代码可以在这里提供帮助吗?有两个相关的考虑因素。第一个是线性索引:通常将其与“存储顺序”混淆,即,对于同一数组的两个线性索引 ij,内存中的偏移量与 i-j 成正比。对于行优先数组,这个概念不可行,因为否则循环

function copy!(dest, src)
    for i = 1:length(src)
        dest[i] = src[i]  # trouble if `i` means "memory offset"
    end
    dest
end

如果 srcdest 使用不同的存储顺序,最终将进行转置。因此,必须根据相应的笛卡尔(全维度)索引来定义线性索引。这不是什么大问题,因为我们知道如何解决:在必要时使用 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 声称返回一个“高效迭代器”,毫无疑问,对于行优先数组,RowMajorRangeCartesianRange 迭代器更有效。所以让我们假设 eachindex 做到了它所说的,并返回一个 RowMajorRange 迭代器。使用此策略,上述两种算法可以合并到一个泛型实现中

for I in eachindex(B)
    dest[I[1]] += B[I]*v[I[2]]
end

耶!高效泛型实现得分。

但我们的胜利是短暂的。让我们回到上面 copy! 的示例,并意识到 destsrc 可能属于两种不同的数组类型,因此可能最有效地使用不同的迭代器类型进行索引。我们很想这样写:

function copy!(dest, src)
    for (idest, isrc) in zip(eachindex(dest), eachindex(src))
        dest[idest] = src[isrc]
    end
    dest
end

在我们为 eachindex 引入 RowMajorRange 返回类型之前,此实现是可以的。但我们刚刚破坏了它,因为现在在某些情况下,这将错误地进行转置。

换句话说,如果没有仔细的设计,“最大效率迭代”和“保持访问同步”的目标就会发生冲突。

OffsetArrays 和 AbstractArray 的含义

Julia 的数组从 1 开始索引,而其他一些语言从 0 开始编号。如果你根据各种博客文章的评论来判断,那么有大量的程序员迫切希望采用 Julia,但由于索引的这种差异,他们甚至不会尝试。由于招募这些军队将导致世界统治,这显然是一个极其紧迫的问题。

更严肃地说,确实存在一些算法,如果可以索引超出 1:size(A,d) 的范围,则这些算法会得到简化。在我们自己的实验室内部代码中,我们长期以来一直在使用 CenterIndexedArray 类型,其中此类数组(所有数组的大小都是奇数)在 -n:n 范围内进行索引,其中 0 指的是“中心”元素。OffsetArrays 包推广了这个概念。不幸的是,在实践中,对于 Julia 的许多核心函数,这两种数组类型都会产生段错误(由于对 @inbounds 何时合适做出了内置假设);随着时间的推移,我们的实验室不得不为相当多的 Julia 函数编写专门针对 CenterIndexedArrays 的实现。

OffsetArrays 说明了另一个概念上的挑战,可以通过 copy! 轻松地演示出来。当 dest 是一个一维 OffsetArraysrc 是一个标准 Vector 时,copy! 应该做什么?特别是,src[1] 放在哪里?它放在 dest 的第一个元素中,还是存储在 dest[1] 中(这可能不是第一个元素)。

此类示例迫使我们更深入地思考数组的真正含义。似乎有两个潜在的概念。一个是数组是列表,多维数组是列表的列表的列表的……在这种世界观中,正确的方法是将 src[1] 放入 dest 的第一个槽中,因为 1 只是 first 的同义词。但是,这种世界观并没有真正赋予数组的索引元组任何“意义”,从某种意义上说,它甚至没有包含 OffsetArray 传达的区别。换句话说,在这个世界中,OffsetArray 根本没有意义,不应该存在。

如果相反地认为 OffsetArray 应该存在,这实质上迫使人们采用不同的世界观:数组实际上是关联容器,其中每个索引元组都是检索值的“键”。在这种思维方式下,src[1] 应该存储在 dest[1] 中。

形式化 AbstractArray

这些示例建议对 AbstractArray 进行形式化

寻找前进的道路

解决这些相互冲突的需求并非易事。一种方法可能是规定某些数组类型根本无法使用泛型代码支持。这可能是正确的策略。或者,可以尝试设计一个处理所有这些类型(并希望更多类型)的数组 API。

在 GitHub 问题 #15648 中,我们正在讨论可能解决这些挑战的 API。鼓励读者参与此讨论。