转置矩阵并保持列优先的内存布局

Question

问题说明：矩阵的行范数

考虑这个玩具示例，我计算随机矩阵 M 的所有列的范数

julia> M = rand(Float64, 10000, 10000);

julia> @time map(x -> norm(x), M[:,j] for j in 1:size(M)[2]);
  0.363795 seconds (166.70 k allocations: 770.086 MiB, 27.78% gc time)

然后计算行范数

julia> @time map(x -> norm(x), M[:,i] for i in 1:size(M)[1]);
  1.288872 seconds (176.19 k allocations: 770.232 MiB, 0.37% gc time)

两次执行之间的因素是（我认为）矩阵的内存布局（主要列）。实际上，行规范的计算是对非连续数据的循环，这会导致缓存未命中的非向量化代码。 我希望两种规范计算的执行时间相同。

在计算行的范数时是否可以将M的布局转换为行主要以获得相同的速度？

我做了什么

我尝试使用 transpose 和 permutedims 没有成功，似乎在使用这些函数时，内存现在是行优先的（所以原始矩阵的主要列）。

julia> Mt = copy(transpose(M));

julia> @time map(x -> norm(x), Mt[j,:] for j in 1:size(M)[2]);
  1.581778 seconds (176.19 k allocations: 770.230 MiB)

julia> Mt = copy(permutedims(M,[2,1]));

julia> @time map(x -> norm(x), Mt[j,:] for j in 1:size(M)[2]);
  1.454153 seconds (176.19 k allocations: 770.236 MiB, 9.98% gc time)

我在这里使用copy 来尝试强制使用新布局。

如何强制转置的列优先布局，或原始矩阵的行优先布局？

编辑

正如@mcabbott 和@przemyslaw-szufel 所指出的，我显示的最后一个代码中有一个错误，我计算了Mt 行的范数而不是列的范数。

对Mt的列范数的检验改为：

julia> Mt = transpose(M);
julia> @time map(x -> norm(x), M[:,j] for j in 1:size(M)[2]);
  1.307777 seconds (204.52 k allocations: 772.032 MiB, 0.45% gc time)

julia> Mt = permutedims(M)
julia> @time map(x -> norm(x), M[:,j] for j in 1:size(M)[2]); 
  0.334047 seconds (166.53 k allocations: 770.079 MiB, 1.42% gc time)

所以最后看来permutedims 以列为主存储，正如预期的那样。事实上，Julia 数组总是以列为主。 transpose 是一种例外，因为它是列优先存储矩阵的行优先view。

Answer 1

这里有几个问题：

• 您错误地对代码进行基准测试 - 很可能在第一次运行时测试已编译的代码，而在第二次运行时测试未编译的代码（并因此测量编译时间）。您应该始终运行 @time 两次或改用 BenchmarkTools
• 您的代码效率低下 - 进行了不必要的内存复制
• M 的类型不稳定，因此测量包括找出其类型的时间，这在您通常运行 Julia 函数时不是这种情况
• 您不需要 lambda - 您可以直接解析函数。
• 正如@mcabbott 所提到的，您的代码包含一个错误，并且您测量的是同一事物的两倍。 清理后您的代码如下所示：

julia> using LinearAlgebra, BenchmarkTools

julia> const M = rand(10000, 10000);

julia> @btime map(norm, @view M[:,j] for j in 1:size(M)[2]);
  49.378 ms (2 allocations: 78.20 KiB)

julia> @btime map(norm, @view M[i, :] for i in 1:size(M)[1]);
  1.013 s (2 allocations: 78.20 KiB)

现在是关于数据布局的问题。 Julia 正在使用以列为主的内存布局。因此，在列上工作的操作将比在行上工作的更快。 一种可能的解决方法是拥有M 的转置副本：

const Mᵀ = collect(M')

这需要一些时间来复制，但可以让您稍后匹配性能：

julia> @btime map(norm, @view Mᵀ[:,j] for j in 1:size(M)[2]);
  48.455 ms (2 allocations: 78.20 KiB)

julia> map(norm, Mᵀ[:,j] for j in 1:size(M)[2]) == map(norm, M[i,:] for i in 1:size(M)[1])
true

Answer 2

在计算规范时，您会浪费大量时间来创建每列/行的副本。改用views，或者更好的是eachcol/eachrow，它们也不会分配：

julia> M = rand(1000, 1000);

julia> @btime map(norm, $M[:,j] for j in 1:size($M, 2));  # slow and ugly
  946.301 μs (1001 allocations: 7.76 MiB)

julia> @btime map(norm, eachcol($M));  # fast and nice
  223.199 μs (1 allocation: 7.94 KiB)

julia> @btime norm.(eachcol($M));  # even nicer, but allocates more for some reason.
  227.701 μs (3 allocations: 47.08 KiB)