22.1.4.3 数学考虑¶

已尽力使稀疏矩阵的行为与其对应的满矩阵完全相同。然而，确实存在一些差异，尤其是与其他产品的稀疏实现相比存在差异。

首先，"./" 和 ".^" 运算符必须谨慎使用。考虑以下示例会产生什么结果：

  s = speye (4);
  a1 = s .^ 2;
  a2 = s .^ s;
  a3 = s .^ -2;
  a4 = s ./ 2;
  a5 = 2 ./ s;
  a6 = s ./ s;

第一个示例中，s 被提升到2的幂没有任何问题。然而，s 逐元素地自乘涉及大量 0 .^ 0 项，其结果为1。因此 s .^ s 是一个满矩阵。

同样地，s .^ -2 涉及诸如 0 .^ -2 这样的项，结果为无穷大，所以 s .^ -2 同样是一个满矩阵。

对于 "./" 运算符，s ./ 2 没有问题，但是 2 ./ s 也涉及大量无穷大项，同样是一个满矩阵。s ./ s 涉及诸如 0 ./ 0 这样的项，结果为 NaN，因此这也是一个满矩阵，其中 s 的零元素被填充为 NaN 值。

上述行为与满矩阵一致，但与其他产品中的稀疏实现不一致。

稀疏矩阵的一个特殊问题在于：由于零没有被存储，这些零的符号位也同样没有被存储。在某些情况下，零的符号位是重要的。例如：

 a = 0 ./ [-1, 1; 1, -1];
 b = 1 ./ a
 ⇒   -Inf            Inf
     Inf           -Inf
 c = 1 ./ sparse (a)
 ⇒   Inf            Inf
     Inf            Inf

要纠正这种行为，就需要将符号位为负的零元素存储在矩阵中，以确保它们的符号位得到保留。出于效率原因，目前尚未这样做，因此警告用户：在零的符号位很重要的计算中，不能使用稀疏矩阵。

一般来说，在稀疏矩阵上使用的任何函数或运算符，其结果稀疏矩阵的非零元素数量与原始矩阵相同或更多。这对于稀疏矩阵分解这一重要情况尤其如此。通常的处理方法是对矩阵进行重新排序，使其分解比原始矩阵的分解更稀疏。即，分解 L * U = P * S * Q 的项 L 和 U 比等价分解 L * U = S 更稀疏。

根据待分解矩阵的类型，有几个函数可用于重新排序。如果矩阵是对称正定的，则应使用 symamd 或 csymamd。否则应使用 amd、colamd 或 ccolamd。为完整起见，还提供了重排序函数 colperm 和 randperm。

请参阅 图22.3，这是一个简单正定矩阵结构的示例。

图22.3：简单稀疏矩阵的结构。

该矩阵的标准 Cholesky 分解可以通过与满矩阵相同的命令获得：r = chol (A); spy (r);。请参阅 图22.4。原始矩阵有 598 个非零项，而此 Cholesky 分解有 10200 个非零项，且仅存储了对称矩阵的一半。这是一个显著的填充水平，虽然对于这样一个小测试用例来说不是问题，但在处理其他稀疏矩阵时可能代表巨大的开销。

原始矩阵的适当保稀疏排列由 symamd 给出，使用此重排序的分解可以通过命令 q = symamd (A); r = chol (A(q,q)); spy (r) 可视化。这给出了 399 个非零项，是一个显著的改进。

Cholesky 分解本身可用于在分解过程中确定矩阵的适当保稀疏重排序。在这种情况下，可以通过三个返回值获得：[r, p, q] = chol (A); spy (r)。

图22.4：上述矩阵未经置换的 Cholesky 分解的结构。

图22.5：上述矩阵经置换的 Cholesky 分解的结构。

在非对称矩阵的情况下，适当的保稀疏排列是 colamd，使用此重排序的分解可以通过命令 q = colamd (A); [l, u, p] = lu (A(:,q)); spy (l+u) 可视化。

最后，当使用 div (/) 和 ldiv (\) 运算符时，Octave 会隐式地对矩阵进行重排序，因此用户无需显式地重排矩阵以最大化性能。

p = amd (S) ¶

p = amd (S, opts) ¶

返回矩阵的近似最小度（Approximate Minimum Degree）置换。

这是一个置换，使得 S (p, p) 的 Cholesky 分解往往比 S 本身的 Cholesky 分解更稀疏。amd 通常比 symamd 更快，但具有类似的目的。

可选参数 opts 是一个控制 amd 行为的结构体。该结构体的字段为：

opts.dense

确定 amd 认为输入矩阵的密集行或列的标准。超过 max (16, (dense * sqrt (n))) 个元素的行或列（其中 n 是矩阵 S 的阶数）在 amd 计算置换时会被忽略。密集值必须是正标量，默认值为 10.0。

opts.aggressive

如果该值为非零标量，则 amd 会进行激进吸收。默认情况是不进行激进吸收。

代码本身的作者是 Timothy A. Davis（请参阅 http://faculty.cse.tamu.edu/davis/suitesparse.html）。

另请参阅： symamd, colamd.

p = ccolamd (S) ¶

p = ccolamd (S, knobs) ¶

p = ccolamd (S, knobs, cmember) ¶

[p, stats] = ccolamd (S) ¶

[p, stats] = ccolamd (S, knobs) ¶

[p, stats] = ccolamd (S, knobs, cmember) ¶

返回约束的近似最小度（Approximate Minimum Degree）置换。

对于矩阵 S 的 ccolamd 工作方式与 colamd 大致相同，但还可以选择使用 cmember 向量限制每列的排序。如果希望某些列组在其他列组之前进行排序，这非常有用。cmember 可以用于在排序过程中部分地保持矩阵结构。

代码本身的作者是 Stefan I. Larimore 和 Timothy A. Davis。该算法是与施乐 PARC 的 John Gilbert 和橡树岭国家实验室的 Esmond Ng 合作开发的。（请参阅 http://faculty.cse.tamu.edu/davis/suitesparse.html）

另请参阅： colamd, csymamd.

p = colamd (S) ¶

p = colamd (S, knobs) ¶

[p, stats] = colamd (S) ¶

[p, stats] = colamd (S, knobs) ¶

返回列近似最小度（Column Approximate Minimum Degree）置换。

p 是列交换后的置换向量，使得矩阵 S (:, p) 的 LU 分解倾向于比 S 的 LU 分解更稀疏。colamd 还可以用于对称矩阵的 Cholesky 分解。

代码本身的作者是 Stefan I. Larimore 和 Timothy A. Davis。该算法是与施乐 PARC 的 John Gilbert 和橡树岭国家实验室的 Esmond Ng 合作开发的。（请参阅 http://faculty.cse.tamu.edu/davis/suitesparse.html）

另请参阅： colperm, symamd, ccolamd.

p = csymamd (S) ¶

p = csymamd (S, knobs) ¶

[p, stats] = csymamd (S) ¶

[p, stats] = csymamd (S, knobs) ¶

返回约束的对称近似最小度（Constrained Symmetric Approximate Minimum Degree）置换。

与 symamd 类似，但 csymamd 允许将排序约束为与唯一的整数变量数组 S 的相等节点集。如果 S = speye (4)，则没有约束；csymamd 的行为与 symamd 完全相同。

代码本身的作者是 Stefan I. Larimore 和 Timothy A. Davis。该算法是与施乐 PARC 的 John Gilbert 和橡树岭国家实验室的 Esmond Ng 合作开发的。（请参阅 http://faculty.cse.tamu.edu/davis/suitesparse.html）

另请参阅： symamd, colamd.

p = colperm (S) ¶

返回矩阵的列置换，该置换按非零元素计数对列进行排序。

该函数返回列置换向量 p，使得 S (:, p) 的列按非零元素数量升序排列。

另请参阅： colamd, symamd.

p = randperm (n) ¶

p = randperm (n, m) ¶

返回一个包含从1到n的整数的随机排列的行向量。

如果调用时有两个参数，则返回一个长度为 m 的随机排列，其中 m ≤ n。

另请参阅： colperm.

p = symamd (S) ¶

p = symamd (S, knobs) ¶

[p, stats] = symamd (S) ¶

[p, stats] = symamd (S, knobs) ¶

对于对称正定矩阵 S，返回置换向量 p，使得 S (p, p) 的 Cholesky 因子倾向于比 S 的更稀疏。

有时 symamd 也适用于对称不定矩阵。矩阵 S 被假定为对称的；仅引用严格下三角部分。S 必须是方阵。

knobs 是一个可选的一到二元素输入向量。如果 S 是 n×n 矩阵，则元素个数超过 max(16, knobs(1)*sqrt(n)) 的行和列在排序之前被删除，并在输出排列 p 中最后排序。如果 knobs(1) < 0，则 stats 和 knobs 会被打印。默认为 knobs = [10 0]。注意，knobs 与早期版本的 symamd 不同。

stats 是一个可选的20元素输出向量，提供关于输入矩阵 S 的排序和有效性的数据。排序统计数据位于 stats(1:3)。stats(1) 和 stats(2) 是 SYMAMD 忽略的密集或空行和列的数量，stats(3) 是对 SYMAMD 使用的内部数据结构执行的垃圾回收次数（大致为 8.4 * nnz (tril (S, -1)) + 9 * n 个整数）。

Octave 内置函数用于生成有效的稀疏矩阵，没有重复条目，每列中非零元的行索引按升序排列，每列中的条目数为非负数（！），依此类推。如果矩阵无效，则 SYMAMD 可能无法继续。如果存在重复条目（同一列中出现两次或多次相同的行索引），或者一列中的行索引顺序错乱，则 SYMAMD 可以通过忽略重复条目并对矩阵 S 的内部副本的每列进行排序来纠正这些错误（但是输入矩阵 S 不会被修复）。如果矩阵以其他方式无效，则 SYMAMD 无法继续，将打印错误消息，并且不返回任何输出参数（p 或 stats）。因此，SYMAMD 是检查稀疏矩阵是否有效的一种简单方法。

stats(4:7) 提供信息，指示 SYMAMD 是否能够继续。如果 stats(4) 为零，则矩阵正常；如果为 1，则无效。stats(5) 是未排序或包含重复项的最右侧列索引，如果不存在此类列，则为零。stats(6) 是在 stats(5) 给出的列索引中最后出现的重复或无序行索引，如果不存在此类行索引，则为零。stats(7) 是重复或无序行索引的数量。stats(8:20) 在 SYMAMD 的当前版本中始终为零（保留以供将来使用）。

排序之后是列消除树后排序。

代码本身的作者是 Stefan I. Larimore 和 Timothy A. Davis。该算法是与施乐 PARC 的 John Gilbert 和橡树岭国家实验室的 Esmond Ng 合作开发的。（请参阅 http://faculty.cse.tamu.edu/davis/suitesparse.html）

另请参阅： colperm, colamd.

p = symrcm (S) ¶

返回 S 的对称逆 Cuthill-McKee（Reverse Cuthill-McKee）置换。

p 是一个置换向量，使得 S (p, p) 的对角线元素往往比 S 更接近对角线。这是对来自"长而瘦"问题的矩阵进行 LU 或 Cholesky 分解的一个良好预排序。它适用于对称和非对称的 S。

该算法代表了 NP-完全带宽最小化问题的一种启发式方法。实现基于以下文献中的描述：

E. Cuthill, J. McKee. "Reducing the Bandwidth of Sparse Symmetric Matrices". Proceedings of the 24th ACM National Conference, pp. 157–172, 1969, Brandon Press, New Jersey.

A. George, J.W.H. Liu. "Computer Solution of Large Sparse Positive Definite Systems", Prentice Hall Series in Computational Mathematics, ISBN 0-13-165274-5, 1981.

另请参阅： colperm, colamd, symamd.