为什么稀疏矩阵(0.5)与密集向量相乘比正常乘法慢？

Question

我生成了一个稀疏度为0.5的矩阵,我想知道它的速度乘以密集向量与正常方式的乘法相比.

function sparseTest()
    A = sprand(1000, 1000, 0.5);
    test(A);
    test(full(A));
end
function test(A)
    tic;
    b = rand(size(A, 2), 1);
    for i = 1:10000
        c = A * b;
    end
    toc;
end

结果是

Elapsed time is 10.968797 seconds.
Elapsed time is 10.194440 seconds.

令我困惑的是,由于只有一半的非零数字,我认为稀疏乘法应该更快.

Answer 1

不,那不是稀疏矩阵.当稀疏矩阵有用时,你会误解.我们这些大量使用它们的人会说,除非密度小于1%,否则矩阵并不是很稀疏.即使密度不是很稀疏!

稀疏矩阵在使用时会有额外的簿记,因为现在MATLAB基本上需要担心非零的位置.因此,如果你创建一个基本上是密集矩阵的矩阵,并假装它是稀疏的,那么你就是在挫败工具的价值.

让我们做一两个比较.首先,一个完全密集的矩阵.我将使用文件交换中的Steve Eddins timeit工具.这在估计所花费的时间方面做得非常好,并且不需要显式循环.

A = rand(1000);
As = sparse(A);
b = rand(1000,1);

whos A As
  Name         Size                 Bytes  Class     Attributes
  A         1000x1000             8000000  double              
  As        1000x1000            16008008  double    sparse    

看到完全密集的矩阵比稀疏形式占用的空间要小得多.

timeit(@() A*b)
ans =
   0.00039474

timeit(@() As*b)
ans =
    0.0023663

简单乘法中稀疏矩阵的开销太大.

As = sprand(1000,1000,.5);
A = full(As);
b = rand(1000,1);

whos A As
Name         Size                Bytes  Class     Attributes
  A         1000x1000            8000000  double              
  As        1000x1000            6303448  double    sparse    

这两个矩阵几乎没有显示稀疏版本的大小减少.因此,将其称为稀疏是浪费时间的浪费.同样,正如您所发现的那样,乘法所消耗的时间仍然不是增益.

timeit(@() A*b)
ans =
   0.00027594

timeit(@() As*b)
ans =
    0.0005011

然而,当矩阵更接近我在稀疏性方面的建议时会发生什么？

As = sprand(1000,1000,.01);
A = full(As);
b = rand(1000,1);

whos A As
  Name         Size                Bytes  Class     Attributes

  A         1000x1000            8000000  double              
  As        1000x1000             167176  double    sparse   

timeit(@() A*b)
ans =
   0.00023629

timeit(@() As*b)
ans =
   5.3167e-05

好的,所以现在这样可以节省一些空间.还有一些严肃的时间.

但实际上,乘法的时间并不是我们使用稀疏矩阵的原因.这是解决线性方程组所需的时间,这使得我们使用稀疏矩阵.这可能是一个巨大的节省,使用完整的矩阵基本上无法解决问题.对我来说,一个常见问题可能涉及40000x40000矩阵,密度可能为0.02%.我会花几天时间等待完整版本完成,但在这种情况下只需要几秒钟即可完成稀疏解决.

例如,考虑gridfit,这是我在文件交换上的函数.它从分散的数据创建表面.在此示例中,该工具将需要求解具有40000个未知数的线性方程组.

timeit(@() gridfit(rand(100,1),rand(100,1),rand(100,1),linspace(0,1,200),linspace(0,1,200)))
ans =
      0.83624

因此,对于所有事情,包括设置方程组,它花了不到1秒钟.

它解决的(在这种情况下,非方形)线性系统是非常稀疏的.

whos A
  Name          Size                 Bytes  Class     Attributes
  A         80100x40000            4126408  double    sparse

它有多稀疏？

nnz(A)
ans =
      237900

nnz(A)/40000
ans =
       5.9475

所以平均来说,每行只有5.9475个非零,有40000列.我们没有时间对该系统进行全线性求解.至少我没有任何尝试它的愿望.