行优先顺序对于矩阵向量乘法是否更有效？

Question

如果 M 是一个 n x m 矩阵，而 v 和 u 是向量，那么就索引而言，矩阵向量乘法看起来像 u[i] = sum(M[i,j] v_j, 1 <= j <= m)。由于 v 是一个向量，它的元素大概存储在面向数值计算的语言的连续内存位置中。如果 M 以行优先顺序存储（如在 C、Mathematica 和 Pascal 中），则总和中的后续 M[i,j] 也存储在连续的内存位置，因为 j 递增，使迭代非常有效。如果它以列优先顺序存储（如在 Fortran、Matlab、R 和 Julia 中），则递增 j 需要移动等于外部矩阵步幅的内存位置，在本例中等于 n。对于多行的矩阵，这看起来效率很低。（对于矩阵 - 矩阵乘法，问题不会出现，因为在任何一种排序约定下，递增求和索引需要在一个矩阵或另一个矩阵的内存中移动主要步幅。）

与乘法和加法运算相比，在大多数计算机体系结构中，在内存中移动一个单位和移动多个单位之间的区别是明显的还是可以忽略不计？（我猜 "negligible"，因为实际上 Fortran 通常至少和 C 一样快，但谁能详细说明原因？）

Answer 1

在大多数计算机体系结构中，这种差异预计会很大，至少在原则上是这样。

矩阵向量乘法是一种内存限制计算，因为内存的重复使用率很低。 v 的所有 (N) 个分量都被重新用于计算 u 的每个元素，但矩阵 (N^2) 的每个元素仅使用一次。如果我们将典型内存的延迟（参见 https://gist.github.com/hellerbarde/2843375）视为（小于）100ns 与执行浮点运算所需的时间（小于 1ns）相比，我们会发现大部分时间都花在了加载和存储值 from/to 个数组。

我们仍然可以实现缓存友好的，即尽可能具有数据局部性。由于内存是按行加载到缓存中的，所以我们必须尽可能使用加载的缓存行。这就是访问连续内存区域减少从内存加载数据所花费的时间的原因。

为了支持这一点，让我们尝试一个非常简单的代码：

program mv
integer, parameter :: n=10000
real, allocatable :: M(:,:), v(:), u(:)
real :: start, finish
integer :: i, j
allocate(M(n,n),v(n),u(n))
call random_number(M)
call random_number(v)
u(:)=0.
call cpu_time(start)
do i=1,n
do j=1,n
    ! non-contiguous order
    u(i)=u(i)+M(i,j)*v(j)
    ! contiguous order
    ! u(i)=u(i)+M(j,i)*v(j)
enddo
enddo
call cpu_time(finish)
print*,'elapsed time: ',finish-start
end program mv

部分结果：

               non-contiguous order   contiguous order
gfortran -O0            1.                 0.5
gfortran -O3           0.3                 0.1
ifort -O0              1.5                0.85
ifort -O3            0.037               0.035

如您所见，区别在于编译时没有优化。启用优化 gfortran 仍然显示出显着差异，而使用 ifort 只有很小的差异。查看编译器报告，似乎编译器交换了循环，从而导致对内部循环的连续访问。

但是，我们可以说具有行优先排序的语言对于矩阵向量计算更有效吗？不，我不能那样说。不仅仅是因为编译器可以补偿差异。代码本身并不知道关于 M 的行和列的所有信息：它基本上知道 M 有两个索引，其中一个——取决于语言——在内存中是连续的。对于矩阵向量，最好的数据局部性是将 "fast" 索引映射到矩阵行索引。您可以使用 "row-major" 和 "column-major" 两种语言来实现。您只需要根据此存储 M 的值。例如，如果您有 "algebraic" 矩阵

     [ M11 M12 ]
M =  [         ]
     [ M21 M22 ]

您将其存储为 "computational matrix"

C       ==> M[1,1] = M11 ; M[1,2] = M12 ; M[2,1] = M21 ; M[2,2] = M22 
Fortran ==> M[1,1] = M11 ; M[2,1] = M12 ; M[1,2] = M21 ; M[2,2] = M22

这样您就可以始终在 "algebraic matrix" 行中保持连续。计算机对初始矩阵一无所知，但我们知道计算矩阵是代数矩阵的转置版本。在这两种情况下，我都会让内部循环遍历连续索引，最终结果将是相同的向量。

在复杂的代码中，如果我已经分配并用值填充了矩阵，但我无法决定存储转置矩阵，则 "row-major" 语言可能会提供最佳性能。但是，交换循环（参见 https://en.wikipedia.org/wiki/Loop_interchange) as automatically done by intel compilers and as done by BLAS implementations (see http://www.netlib.org/lapack/explore-html/db/d58/sgemv_8f_source.html），将差异减少到非常小的差异值。因此，使用 Fortran 你可以更喜欢：

do j=1,n
    do i=1,n
        u(i)=u(i)+M(i,j)*v(j)
    enddo
enddo

行优先顺序对于矩阵向量乘法是否更有效？

Is row-major ordering more efficient for matrix-vector multiplication?

language-agnostic

matrix

linear-algebra

matrix-multiplication