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

Question

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

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

Answer 1

至少在原则上，大多数计算机体系结构的差异预计会很大。

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

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

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

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 有两个索引，其中一个——取决于语言——在内存中是连续的。对于矩阵向量，最适合数据局部性的方法是将“快速”索引映射到矩阵行索引。您可以使用“行主要”和“列主要”语言来实现这一点。您只需要根据此存储 M 的值。例如，如果您有“代数”矩阵

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

您将其存储为“计算矩阵”

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 

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

在复杂的代码中，如果我已经分配并用值填充了矩阵并且我无法决定存储转置矩阵，则“行优先”语言可能会提供最佳性能。但是，由英特尔编译器自动完成循环（参见https://en.wikipedia.org/wiki/Loop_interchange）和由 BLAS 实现完成（参见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