在行上迭代矩阵的推荐方法是啥？

Posted 2023-02-16

【中文标题】在行上迭代矩阵的推荐方法是啥？

【英文标题】：What is the recommended way to iterate a matrix over rows?在行上迭代矩阵的推荐方法是什么？

【发布时间】：2014-03-13 18:08:35

【问题描述】：

给定一个矩阵m = [10i+j for i=1:3, j=1:4]，我可以通过对矩阵进行切片来迭代它的行：

for i=1:size(m,1)
    print(m[i,:])
end

这是唯一的可能吗？这是推荐的方式吗？

那么理解呢？切片是迭代矩阵行的唯一可能性吗？

[ sum(m[i,:]) for i=1:size(m,1) ]

【问题讨论】：

地图切片？ mapslices(sum, m, 2) 是后者

@jverzani mapslices 完成了这项工作，尽管在某些情况下它需要我定义一个匿名函数。感谢您的建议。

【参考方案1】：

您自己列出的解决方案以及mapslices 都可以正常工作。但是，如果“推荐”您真正的意思是“高性能”，那么最好的答案是：不要迭代行。

问题在于，由于数组是以列优先顺序存储的，因此对于小矩阵以外的任何内容，如果您以行优先顺序遍历数组，您最终会得到一个糟糕的cache hit ratio。

正如excellent blog post 中指出的那样，如果您想对行求和，最好的办法是执行以下操作：

msum = zeros(eltype(m), size(m, 1))
for j = 1:size(m,2)
    for i = 1:size(m,1)
        msum[i] += m[i,j]
    end
end

我们以原始存储顺序遍历m 和msum，因此每次加载缓存行时，我们都会使用所有值，从而产生缓存命中率 1。您可能天真地认为遍历它会更好按行优先顺序并将结果累积到 tmp 变量，但在任何现代机器上，缓存未命中比 msum[i] 查找要昂贵得多。

许多采用 region 参数的 Julia 内部算法（例如 sum(m, 2)）会为您处理这个问题。

【讨论】：

我认为这回答了我的问题，我将再等一天接受答案。我非常喜欢这个答案，因为它让我意识到，由于 Julia 是列优先的，我最好将数据向量排列为列而不是行。

您链接到的博文已不存在。请参阅docs.julialang.org/en/release-0.4/manual/performance-tips/…。

404 似乎是由于斜杠造成的。此网址有效：julialang.org/blog/2013/09/fast-numeric

那是正确的。但是 3 维数组呢？例如：A[i,j,k]。顺序是什么？ k --> j --> i 或 i --> j --> k ?

【参考方案2】：

从 Julia 1.1 开始，有用于迭代矩阵的列或行的迭代器实用程序。遍历行：

M = [1 2 3; 4 5 6; 7 8 9]

for row in eachrow(af)
    println(row)
end

将输出：

[1, 2, 3]
[4, 5, 6]
[7, 8, 9]

【讨论】：

有没有办法也可以使用这种方法获取行索引？也就是说，不仅是每一行本身，还有它的索引。

Skumin，你可以使用for (i, row) in enumerate(eachrow(M))。实际上，您可以将enumerate 应用于任何迭代器

使用eachrow(df) 对性能有何影响？它可以与幼稚循环相提并论吗？

【参考方案3】：

根据我的经验，显式迭代比理解要快得多。

遍历列也是一个很好的建议。

此外，您还可以使用新的宏 @simd 和 @inbounds 来进一步加速。

【参考方案4】：

在我的例子中，我不能使用 eachrow 迭代器或嵌套循环，因为我需要用其他东西压缩 eachindex，然后迭代那个 zip 迭代器。因此，我写道：

ncols = size(m, 2)
for i in eachindex(m)
    rowi, coli = fldmod1(i, ncols)
    elem = m[rowi, coli]
end

请注意，这仅适用于 eachindex 返回线性索引的情况。如果eachindex 返回笛卡尔坐标的迭代器，您可能需要迭代1:prod(size(m))。