我应该什么时候使用 `sparse`?

When should I be using `sparse`?

我一直在浏览 Matlab 的 sparse documentation,试图找出是否有任何指导方针来说明何时使用稀疏表示而不是完整表示是有意义的。

例如,我有一个包含大约 30% 非零条目的矩阵 data。我可以检查使用的内存。

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 表示

The computational complexity of sparse operations is proportional to nnz, the number of nonzero elements in the matrix. Computational complexity also depends linearly on the row size m and column size n of the matrix, but is independent of the product m*n, the total number of zero and nonzero elements.

如果不了解更多细节,很难将其与完整矩阵进行比较。

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

Advantages of the CSC format

  • efficient arithmetic operations CSC + CSC, CSC * CSC, etc.
  • efficient column slicing
  • fast matrix vector products (CSR, BSR may be faster)

Disadvantages of the CSC format

  • slow row slicing operations (consider CSR)
  • changes to the sparsity structure are expensive (consider LIL or DOK)

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

我不是使用 sparse 矩阵的专家,但是 Mathworks 确实有 some documentation 与运算和计算效率有关。

他们的计算复杂度描述:

The computational complexity of sparse operations is proportional to nnz, the number of nonzero elements in the matrix. Computational complexity also depends linearly on the row size m and column size n of the matrix, but is independent of the product m*n, the total number of zero and nonzero elements.

The complexity of fairly complicated operations, such as the solution of sparse linear equations, involves factors like ordering and fill-in, which are discussed in the previous section. In general, however, the computer time required for a sparse matrix operation is proportional to the number of arithmetic operations on nonzero quantities.

为了不让您厌烦算法细节,another answer 建议您不必为只有 25% 非零的数组操心稀疏。他们提供了一些代码供您测试。有关详细信息,请参阅他们的 post。

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.

全矩阵上的许多操作使用 BLAS/LAPACK 库调用,这些库调用经过疯狂优化且难以击败。实际上,在可以充分利用 (i) 稀疏性和 (ii) 特殊矩阵结构的特殊情况下,稀疏矩阵上的操作只会优于完整矩阵上的操作。

只是随机使用稀疏可能会让你变得更糟。 示例:哪个更快,将 10000x10000 全矩阵添加到 10000x10000 全矩阵?或者将一个 10000x10000 的完整矩阵添加到一个完全稀疏的(即一切都为零)10000x10000 矩阵?尝试一下!在我的系统上,full + full 更快!

有哪些稀疏 CRUSHES 满的情况示例?

示例 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 变量线性系统是微不足道的。后者基本上是在稀疏情况下得到解决的问题。

总的来说,如果你有特殊的矩阵结构并且可以巧妙地利用稀疏性,就有可能解决非常大的问题,否则这些问题将是棘手的。如果你有一个足够大的特殊问题,有一个足够稀疏的矩阵,并且精通线性代数(以保持稀疏性),稀疏类型矩阵可能非常强大。

另一方面,在没有经过深思熟虑的情况下随意使用稀疏代码几乎肯定会使您的代码变慢。

如果您有一个固定维度的矩阵,那么建立可靠答案的最佳方法就是反复试验。但是,如果您不知道 matrices/vectors 的尺寸,则经验法则是

Your sparse vectors should have effectively constant number of nonzero entries

对于矩阵将意味着

Your N x N sparse matrix should have <= c * N nonzero entries, where c is a constant "much less" than 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的回答中有解释。

UPD: In fact, in a for cycle, you risk to choose a wrong branch only once, at the end of your cycle, since by definition the default branch to choose will be going into the cycle. This amounts in at most 1 pipeline restart per scalar product operation.

现在让我们看看稀疏向量在BLAS中是如何实现的。在那里,每个向量都由两个数组编码:一个值和一个相应的索引,比如

1.7    -0.8    3.6
171     83     215

(加上一个整数告诉假设的长度N)。在 BLAS 文档中指出,索引的排序在这里不起作用,因此 data

-0.8    3.6    1.7
 83     215    171

编码相同的向量。这条评论提供了足够的信息来重构标量积的算法。给定由数据 int[] indices, double[] valuesint[] jndices, double[] walues 编码的两个稀疏向量,将在 "code":

的行中计算它们的标量积
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 的向量已经完全毁了。

UPD. In this code snippet, you have an if branch, and the probability to take the wrong path is 50%. Thus, a scalar product of two sparse vectors will amount in about 0.5 to 1 times the average number of nonzero entries per sparse vector) pipeline restarts, depending on how sparse your vectors are . The numbers are to be adjusted: in an if statement without else, the shortest instruction will be taken first, which is "do nothing", but still.

请注意,最有效的运算是稀疏向量和密集向量的标量积。给定 indicesvalues 的稀疏向量以及 dense 的密集向量,您的代码将看起来像

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

即你会有 indices.length 个条件分支,这很好。

UPD. Once again, I bet you'll have at most one pipeline restart per operation. Note also that afaik in modern multicore processors both alternatives are performed in parallel on two different cores, so that in alternative branches you only need to wait for the longest one to finish.

现在,当矩阵与向量相乘时,您基本上需要向量的#rows 标量积。通过向量乘法将矩阵与矩阵量相乘,在第二个矩阵中采用#((非零)列)矩阵。欢迎您自行解决复杂性。

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

如果将由稀疏行构建的矩阵与密集矩阵或密集列矩阵相乘,您将 "always" 在速度上有优势。您可能希望将稀疏矩阵数据存储为对角线的密集向量 - 所以在 convolution neural networks - and then you'll need completely different algorithms. You may want to make your matrix a block matrix - so does BLAS - and get a reasonable computation boost. You may want to store your data as two matrices - say, a diagonal and a sparse, which is the case for finite element method. You could make use of sparsity for general neural networks (like. fast forward, extreme learning machine or echo state network) 的情况下,如果您总是将行存储矩阵乘以列向量,但要避免矩阵相乘。而且,如果您遵循经验法则,您将 "always" 通过使用稀疏矩阵获得优势 - 它适用于有限元和卷积网络,但不适用于水库计算。