什么情况下应该使用`sparse`？

Question

我一直在查阅Matlab的sparse文档，试图找到何时使用稀疏表示而不是完整表示的指南。

例如，我有一个矩阵data，其中约30％的条目为非零值。我可以检查所使用的内存。

whos data
  Name             Size                 Bytes  Class     Attributes

  data      84143929x11            4394073488  double    sparse    

---

data = full(data);
whos data
  Name             Size                 Bytes  Class     Attributes

  data      84143929x11            7404665752  double              

在这里，我明显节省了内存，但是对于任何30%非零条目的矩阵都是如此吗？50%呢？有没有一个经验法则来决定在什么百分比时应该切换到完整矩阵？

那么计算方面呢？使用稀疏矩阵进行矩阵乘法通常更快还是更慢？Sparse Matrix Operations说：

“稀疏操作的计算复杂度与矩阵中非零元素的数量nnz成正比。计算复杂度还线性依赖于矩阵的行大小m和列大小n，但与零和非零元素总数m*n的乘积无关。”

这很难与完整矩阵进行比较，因为我们需要了解更多细节。

Scipy的稀疏矩阵库解释了每种稀疏格式的优缺点。例如，对于csc_matrix。

CSC 格式的优点

• 高效的算术运算，如CSC + CSC、CSC * CSC等
• 高效的列切片操作
• 快速的矩阵向量乘积（CSR、BSR可能更快）

CSC 格式的缺点

• 行切片操作较慢（考虑使用 CSR）
• 稀疏结构的更改较为昂贵（考虑使用 LIL 或 DOK）

是否存在关于 Matlab 的 sparse 实现的类似信息？如果有，在哪里可以找到它？

- Cecilia

---

如果您想详细了解我的任何答案部分，请告诉我。 - krisdestruction

我想为这个问题提供赏金，但我的声望非常低。@LuisMendo 也许你想要？ - krisdestruction

从稀疏矩阵文档中，选定参考书目：Gilbert, John R., Cleve Moler 和 Robert Schreiber，“Sparse Matrices in MATLAB: Design and Implementation," SIAM J. Matrix Anal. Appl.，Vol. 13，No.1，1992年1月，pp.333-356。这可能是一个很好的起点。 - Sam Roberts

Answer 1

许多针对全矩阵的操作都使用了BLAS/LAPACK库调用，这些代码经过了极度优化并且很难被超越。实际上，在专门能够充分利用矩阵的稀疏性和特殊结构的情况下，对稀疏矩阵进行的操作才可能胜过对全矩阵的操作。

随意使用稀疏矩阵很可能会让你变得更差。例如：哪个更快，将一个大小为10000x10000的全矩阵加到另一个10000x10000的全矩阵中？还是将一个大小完全为零（即全部为空）的10000x10000矩阵与一个10000x10000的全矩阵相加？试一试！在我的系统上，全+全的速度更快！

有哪些情况下稀疏矩阵能完胜全矩阵呢？

例1：解线性方程组A*x=b，其中A是5000x5000的块对角矩阵，由500个5x5的块构成。建立代码：

As = sparse(rand(5, 5));
for(i=1:999)
   As = blkdiag(As, sparse(rand(5,5))); 
end;                         %As is made up of 500 5x5 blocks along diagonal
Af = full(As); b = rand(5000, 1);

然后，您可以测试速度差异：

As \ b % operation on sparse As takes .0012 seconds
Af \ b % solving with full Af takes about 2.3 seconds

总的来说，一个含有5000个变量的线性系统是有些困难的，但是1000个包含5个变量的分离的线性系统是微不足道的。后者基本上是在稀疏情况下得到解决的。

总体故事是，如果你有特殊的矩阵结构并且能够巧妙地利用稀疏性，那么就有可能解决原本难以处理的超大问题。如果你有一个足够大的专用问题，并且有一个足够稀疏的矩阵，并且对线性代数很聪明（以保留稀疏性），那么一种稀疏类型的矩阵就会变得非常强大。

另一方面，随意投入稀疏性而没有深入、仔细思考几乎肯定会使你的代码变慢。

- Matthew Gunn

Answer 2

我不是使用稀疏矩阵的专家，但Mathworks有一些与操作和计算效率相关的文档。

它们的计算复杂度描述如下：

稀疏操作的计算复杂度与矩阵中非零元素数量nnz成比例。计算复杂度也与矩阵的行大小m和列大小n成线性关系，但与总数m*n（包括零和非零元素）无关。

相当复杂的操作的复杂度，如解稀疏线性方程组，涉及到像排序和填充等因素，这些因素在前面的部分中讨论过。然而，通常情况下，对非零量进行算术运算所需的计算机时间与稀疏矩阵操作的数量成比例。

不想让您感到无聊，另一个答案建议，如果数组只有25％的非零值，则不应使用稀疏矩阵。他们提供了一些代码供您测试。请查看他们的帖子以获取详细信息。

A = sprand(2000,2000,0.25);
tic,B = A*A;toc
Elapsed time is 1.771668 seconds.

Af = full(A);
tic,B = Af*Af;toc
Elapsed time is 0.499045 seconds.

- krisdestruction

---

我发现了相同复杂度的文档。但是你如何将其与完整矩阵操作复杂度进行比较呢？例如，假设一个nxm矩阵和一个mxp矩阵的完整矩阵乘法是O(nmp)。如果稀疏矩阵的文档表示复杂度与m和n成线性关系，那么我能得出结论说乘法是O(m + n)还是O(nnz(m) + nnz(n))? - Cecilia

就内存使用而言，我可以扩展该代码片段以适应非零元素变化百分比的函数。这可能很有用，但似乎这些特性也应该是实现固有的，并且我们应该能够证明矩阵实现将使用一定数量的内存，给定一些非零条目数。 - Cecilia

基本上，我想要无聊的算法细节，拜托了。 :) - Cecilia

@Cecilia，我很欣赏喜欢细节的人，但不幸的是我不知道它（正如我所说，我不是专家）。现在，“排列和重新排序”部分似乎描述了它如何对乘法进行排列。也许有人能更有用？ :/ - krisdestruction

Answer 3

如果您有一个固定维度的矩阵，则建立可靠答案的最佳方法就是试错法。然而，如果您不知道矩阵/向量的维度，则有以下经验法则：

您的稀疏向量应该具有有效的非零条目数量

对于矩阵来说，这意味着

您的N x N稀疏矩阵应该具有<= c * N个非零条目，其中c是一个“远小于”N的常数。

让我们为这个规则提供一个伪理论解释。我们将考虑一个相当简单的任务，即使用双值坐标的两个向量的标量（或点）积。现在，如果您有两个相同长度N的密集向量，则您的代码将如下：

//define vectors vector, wector as double arrays of length N 
double sum = 0;
for (int i = 0; i < N; i++)
{
    sum += vector[i] * wector[i];
}

这个过程涉及 N 次加法，N 次乘法和 N 次条件分支（循环操作）。最耗费时间的是条件分支，以至于我们可以忽略乘法和加法。为什么它如此昂贵，这个问题在this question的答案中有所解释。

更新：实际上，在 for 循环中，你只有在循环结束时才会决定选择哪个分支，因为按照定义，进入循环的默认分支将会被选择。这导致每个标量积运算最多只会有一次 pipeline 重启。

现在让我们看一下 BLAS 中如何实现稀疏向量。在那里，每个向量由两个数组编码：一个是值，另一个是相应索引，类似于

1.7    -0.8    3.6
171     83     215

（再加上一个整数，表示所述长度N）。在BLAS文档中指出，索引的排序在此处不起作用，因此数据

-0.8    3.6    1.7
 83     215    171

编码相同的向量。这个备注提供了足够的信息来重构标量积算法。给定由数据int[] indices, double[] values和int[] jndices, double[] walues编码的两个稀疏向量，将按照以下“代码”中的行计算它们的标量积：

double sum = 0;
for (int i = 0; i < indices.length; i++)
{
    for (int j = 0; j < jndices.length; j++)
    {
        if(indices[i] == jndices[j])
        {
            sum += values[indices[i]] * walues[jndices[j]];
        }
    }
}

这给我们提供了一个总数为indices.length * jndices.length * 2 + indices.length的条件分支。这意味着为了应对密集算法，您的向量最多只能有sqrt(N)个非零条目。关键在于依赖于N已经是非线性的，因此没有必要问您需要填充1%、10%还是25%。10%对于长度为10的向量来说是完美的，对于长度为50的向量来说还算可以，但对于长度为100的向量来说就已经彻底毁掉了。

更新：在这段代码片段中，您有一个if分支，错误路径的概率为50%。因此，两个稀疏向量的标量积将导致大约0.5到1倍（每个稀疏向量的平均非零条目数）的管道重启，具体取决于您的向量有多稀疏。这些数字需要进行调整：在没有else的if语句中，将首先执行最短的指令，即"什么也不做"，但仍然会执行。

请注意，最有效率的操作是一个稀疏向量和一个密集向量的标量积。给定一个稀疏向量的 indices 和 values 以及一个密集向量 dense，您的代码将如下所示：

double sum = 0;
for (int i = 0; i < indices.length; i++)
{
    sum += values[indices[i]] * dense[indices[i]];
}

例如，您将拥有indices.length个条件分支，这是很好的。

更新。再次说明，我打赌每个操作最多只会有一个管道重启。请注意，据我所知，在现代多核处理器中，两种替代方案都在两个不同的核心上并行执行，因此在替代分支中，您只需要等待最长的分支完成即可。

现在，当矩阵与向量相乘时，您基本上会取#行标量向量积。矩阵与矩阵相乘意味着取第二个矩阵中的#（非零）列的矩阵乘以向量乘积。欢迎您自己计算复杂度。

这就是所有关于不同矩阵存储的深层理论的黑魔法开始的地方。您可以将稀疏矩阵存储为稠密的稀疏行数组，稀疏的稠密行数组或稀疏的稀疏行数组。列也是一样。问题中引用的Scipy中的所有有趣缩写都与此有关。

如果您用稀疏行构建一个矩阵并与一个密集矩阵或密集列的矩阵相乘，那么您将“始终”具有速度优势。您可能希望将稀疏矩阵数据存储为对角线的密集向量 - 因此在卷积神经网络的情况下 - 然后您将需要完全不同的算法。您可能希望将矩阵转换为块矩阵 - 就像BLAS一样 - 并获得合理的计算提升。您可能希望将数据存储为两个矩阵 - 比如一个对角线和一个稀疏矩阵，这是有限元法的情况。如果您始终将一行存储矩阵乘以列向量，而避免矩阵乘法，则可以利用疏密性进行常规神经网络（例如快速前馈、极限学习机或回声状态网络），但是请避免矩阵乘法。并且，如果您遵循经验法则，始终使用稀疏矩阵，那么您将“始终”获得优势 - 对于有限元和卷积网络有效，但对于储备计算而言则无效。

- Kolya Ivankov