稀疏线性代数 (GNU Octave(版本11.1.0))

下一个：迭代技术在稀疏矩阵中的应用，前一个：稀疏矩阵的创建与运算，上级：稀疏矩阵 [目录][索引]

22.2 稀疏矩阵的线性代数 ¶

Octave 包含一个稀疏矩阵的多态求解器，用于分解矩阵的具体求解器取决于稀疏矩阵本身的特性。通常，确定矩阵类型的成本相对于分解矩阵本身的成本较小，但在任何情况下，矩阵类型一旦计算出来就会被缓存，因此每次在线性方程中使用时不会重新确定。

线性方程求解的选择树如下：

如果矩阵是对角矩阵，直接求解并跳转到 8
如果矩阵是置换对角矩阵，考虑置换关系直接求解。跳转到 8
如果矩阵是方阵、带状矩阵，且带状密度小于 spparms ("bandden") 给出的值，则继续，否则跳转到 4。
1. 如果矩阵是三对角矩阵且右侧不是稀疏矩阵，则继续，否则跳转到 3b。
  1. 如果矩阵是埃尔米特矩阵且具有正实对角元，尝试使用 LAPACK xPTSV 进行 Cholesky 分解。
  2. 如果上述分解失败或矩阵不是具有正实对角元的埃尔米特矩阵，则使用带主元的的高斯消去法（LAPACK xGTSV），并跳转到 8。
2. 如果矩阵是埃尔米特矩阵且具有正实对角元，尝试使用 LAPACK xPBTRF 进行 Cholesky 分解。
3. 如果上述分解失败或矩阵不是具有正实对角元的埃尔米特矩阵，则使用带主元的高斯消去法（LAPACK xGBTRF），并跳转到 8。
如果矩阵是上三角或下三角矩阵，执行稀疏前向或后向替换，并跳转到 8
如果矩阵是带列置换的上三角矩阵或带行置换的下三角矩阵，执行稀疏前向或后向替换，并跳转到 8
如果矩阵是方阵、埃尔米特矩阵且具有正实对角元，尝试使用 CHOLMOD 进行稀疏 Cholesky 分解。
如果稀疏 Cholesky 分解失败，或矩阵不是具有正实对角元的埃尔米特矩阵，且矩阵是方阵，则使用 UMFPACK 进行分解、求解并执行一次精化迭代。
如果矩阵不是方阵，或者之前的任何求解器标记了奇异或近似奇异矩阵，则使用 CXSPARSE¹⁰ 寻找最小范数解。

带状密度定义为带状区域中非零值的数量除以整个带状区域中值的总数。可以通过使用 spparms 将 bandden 设置为 1（即 spparms ("bandden", 1)）来完全禁用带状矩阵求解器。

QR 求解器使用 Dulmage-Mendelsohn 分解将问题分解为多个块，这些块可以处理为超定块、多个良定块以及最后一个超定块。对于具有强连通节点块的矩阵，这是一个巨大的优势，因为可以对许多块使用 LU 分解。它还显著提高了找到超定问题解的机会，而不是仅仅返回一个 NaN 向量。

上述所有求解器都可以计算条件数的估计值。这可用于检测解中的数值稳定性问题，并强制使用最小范数解。然而，对于窄带状、三角或对角矩阵，计算条件数的成本很高，实际上可能超过分解矩阵的成本。因此，在这些情况下不会计算条件数，Octave 依赖更简单的技术来检测奇异矩阵，或者在带状矩阵的情况下依赖底层的 LAPACK 代码。

用户可以使用 matrix_type 函数强制指定矩阵的类型。这避免了确定矩阵类型的开销。但需要注意的是，错误地识别矩阵类型会导致不可预测的结果，因此应谨慎使用 matrix_type。

nest = normest (A) ¶

nest = normest (A, tol) ¶

[nest, iter] = normest (…) ¶

使用幂级数分析估计矩阵 A 的 2-范数。

这通常用于大型矩阵，此时计算 norm (A) 的成本过高，而 2-范数的近似值是可以接受的。

tol 是计算 2-范数时的容差。默认情况下，tol 为 1e-6。

可选的输出 iter 返回 normest 收敛所需的迭代次数。

另请参阅： normest1，norm，cond，condest。

nest = normest1 (A) ¶

nest = normest1 (A, t) ¶

nest = normest1 (A, t, x0) ¶

nest = normest1 (Afcn, t, x0, p1, p2, …) ¶

[nest, v] = normest1 (A, …) ¶

[nest, v, w] = normest1 (A, …) ¶

[nest, v, w, iter] = normest1 (A, …) ¶

使用块算法估计矩阵 A 的 1-范数。

normest1 最适合只需要范数近似值的大型稀疏矩阵。对于中小型矩阵，请考虑使用 norm (A, 1)。此外，当矩阵-向量乘积 A * x 和 A' * x 可以廉价计算时，normest1 可用于估计线性算子 A 的 1-范数。在这种情况下，不使用矩阵 A，而是使用函数 Afcn (flag, x)；它必须返回：

如果 flag 为 "dim"，返回 A 的维度 n
如果 flag 为 "real"，返回 A 是否为实算子
如果 flag 为 "notransp"，返回结果 A * x
如果 flag 为 "transp"，返回结果 A' * x

典型的情况是 A 由 b ^ m 定义，此时结果 A * x 可以在无需显式形成 b ^ m 的情况下计算：

y = x;
for i = 1:m
  y = b * y;
endfor

参数 p1、p2、… 是 Afcn (flag, x, p1, p2, …) 的参数。

t 的默认值为 2。该算法需要大小为 n x n 和 n x t 的矩阵-矩阵乘积。

初始矩阵 x0 的列应具有单位 1-范数。默认的初始矩阵 x0 的第一列为 ones (n, 1) / n，并且如果 t > 1，其余列包含随机元素 -1 / n、1 / n，并除以 n。

输出时，nest 是所需的估计值，v 和 w 是满足 w = A * v 且 norm (w, 1) = c * norm (v, 1) 的向量。iter 中 iter(1) 包含迭代次数（最大值为硬编码的 5），iter(2) 包含算法执行的乘积 A * x 或 A' * x 的总次数。

算法说明：normest1 在求值过程中使用随机数。因此，如果需要一致的结果，应在调用 normest1 之前固定随机生成器的 "state"。

参考文献：N. J. Higham and F. Tisseur， "A block algorithm for matrix 1-norm estimation, with and application to 1-norm pseudospectra"， SIAM J. Matrix Anal. Appl.，Vol. 21, No. 4， pp. 1185–1201, 2000。

另请参阅： normest，norm，cond，condest。

cest = condest (A) ¶

cest = condest (A, t) ¶

cest = condest (A, Ainvfcn) ¶

cest = condest (A, Ainvfcn, t) ¶

cest = condest (A, Ainvfcn, t, p1, p2, …) ¶

cest = condest (Afcn, Ainvfcn) ¶

cest = condest (Afcn, Ainvfcn, t) ¶

cest = condest (Afcn, Ainvfcn, t, p1, p2, …) ¶

[cest, v] = condest (…) ¶

使用 t 个测试向量和随机 1-范数估计器，估计方阵 A 的 1-范数条件数。

可选输入 t 指定测试向量的数量（默认为 5）。

输入可以是矩阵 A（该算法特别适合大型稀疏矩阵）。或者，矩阵的行为可以通过函数隐式定义。使用隐式定义时，condest 需要以下函数：

Afcn (flag, x) 必须返回：
- 如果 flag 是 "dim"，返回 A 的维度 n
- 如果 flag 是 "real"，返回 A 是否为实算子
- 如果 flag 是 "notransp"，返回结果 A * x
- 如果 flag 是 "transp"，返回结果 A' * x
Ainvfcn (flag, x) 必须返回：
- 如果 flag 是 "dim"，返回 inv (A) 的维度 n
- 如果 flag 是 "real"，返回 inv (A) 是否为实算子
- 如果 flag 是 "notransp"，返回结果 inv (A) * x
- 如果 flag 是 "transp"，返回结果 inv (A)' * x

任何参数 p1、p2、… 是 Afcn (flag, x, p1, p2, …) 和 Ainvfcn (flag, x, p1, p2, …) 的额外参数。

主要输出是 1-范数条件数估计值 cest。

可选的第二个输出 v 是近似零向量；它满足方程 norm (A*v, 1) == norm (A, 1) * norm (v, 1) / cest。

算法说明：condest 使用随机化算法来近似 1-范数。因此，如果需要一致的结果，应在调用 condest 之前固定随机生成器的 "state"。

参考文献：

N.J. Higham and F. Tisseur, "A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra"， SIAM Journal on Matrix Analysis and Applications，Vol. 21， Iss. 4，pp. 1185–1201, 2000， https://dx.doi.org/10.1137/S0895479899356080。
N.J. Higham and F. Tisseur，A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra， https://citeseer.ist.psu.edu/223007.html。

另请参阅： cond，rcond，norm，normest1，normest。

vals = spparms () ¶

vals = spparms (key) ¶

vals = spparms ()（过时语法） ¶

val = spparms (key)（过时语法） ¶

key = spparms (key, val)（过时语法） ¶

vals = spparms (keys, vals) ¶

vals = spparms (key, val) ¶

查询或设置 Octave 使用的稀疏矩阵求解器的参数。

当使用多个键/值对调用时，返回一个结构体，或设置稀疏求解器的参数。

以下参数是有效的：

"spumoni": 稀疏求解器诊断信息的打印级别（0|1|2|3）。
"ths_rel": 当非零值的数量与元素总数的比值低于此阈值时，使用稀疏矩阵算法。默认为 0.01。
"ths_abs": 当矩阵的维数低于此阈值时，使用稠密矩阵算法。默认为 100。
"exact_d": 非零值的数量与元素总数的比值，用于确定是否实际需要按行或按列重新排序。默认为 0.3。
"bandden": 用于尝试使用 LAPACK 带状求解器的阈值。默认为 0.5。
"umfpack": umfpack 求解器的启用标志（0|1）。默认为 1。

另请参阅： pbasename，pcg，pcr，pbasename。

d = eigs (A) ¶

d = eigs (A, k) ¶

d = eigs (A, k, sigma) ¶

d = eigs (A, k, sigma, opts) ¶

d = eigs (A, k, sigma, opts, B) ¶

[V, D] = eigs (A, k, sigma, opts, B) ¶

[V, D, flag] = eigs (A, k, sigma, opts, B) ¶

计算矩阵 A 的几个特征值。

默认情况下，使用 ARPACK 计算模最大的 k 个特征值。

如果给出了 opts，它是一个定义了 eigs 应该使用的参数的 structure。opts 结构体的字段是：

issym: 如果给出了 Af 则此标志（true/false）确定函数 Af 是否定义了一个对称问题。如果给出了矩阵 A，则该标志由矩阵的类型确定。默认值为 false。
isreal: 如果给出了 Af 则此标志（true/false）确定函数 Af 是否定义了一个实值问题。如果给出了矩阵 A，则该标志由矩阵的类型确定。默认值为 true。
tol: 定义所需的收敛容差，计算为 tol * norm (A)。默认值为 eps。
maxit: 最大迭代次数。默认值为 300。
p: 要使用的 Lanczos 基向量数量。更多的向量会导致更快的收敛，但会使用更多内存。p 的最佳值取决于问题，应在 k + 1 到 n 的范围内。默认值为 2 * k。
v0: 算法的起始向量。接近最终解的初始向量将加快收敛速度。默认情况下，ARPACK 会随机生成一个起始向量。如果指定了，v0 必须是一个 n-by-1 向量，其中 n = rows (A)。
disp: 诊断打印输出的级别（0|1|2）。如果 disp 为 0，则禁用诊断输出。默认值为 0。
cholB: 如果正在计算广义特征值问题，此标志（true/false）指定 B 输入表示 chol (B) 还是仅仅是矩阵 B。默认值为 false。
permB: 如果 cholB 为 true，则是 B 的 Cholesky 分解的置换向量。它通过 [R, ~, permB] = chol (B, "vector") 获得。默认值为 1:n。

A 也可以通过表示为 Af 来指定。Af 后面必须跟一个标量参数 n，用于定义 Af 接受的向量参数的长度。Af 可以是一个函数句柄、内联函数或字符串。当 Af 是一个字符串时，它指定了要使用的函数名称。

Af 是一个形式为 y = Af (x) 的函数，其中 Af 的返回值由 sigma 的值决定。四种可能的形式是：

A * x: 如果 sigma 未给出或是一个不是 "sm" 的字符串。
A \ x: 如果 sigma 是 0 或 "sm"。
(A - sigma * I) \ x: 如果 sigma 是一个不等于 0 的标量；I 是与 A 大小相同的单位矩阵。
(A - sigma * B) \ x: 用于广义特征值问题。

返回参数及其形式取决于请求的返回参数数量。对于单个返回参数，返回一个长度为 k 的列向量 d，其中包含已找到的 k 个特征值。对于两个返回参数，V 是一个 n-by-k 矩阵，其列是对应于返回特征值的 k 个特征向量。特征值本身以 k-by-k 矩阵的形式在 D 中返回，其中对角线上的元素就是特征值。

第三个返回参数 flag 返回收敛状态。如果 flag 为 0，则所有特征值都已收敛。任何其他值表示未能收敛。

编程说明：对于小问题（n < 500），请考虑使用 eig (full (A))。

如果 ARPACK 未能收敛，请考虑增加 Lanczos 向量的数量（opt.p）、增加迭代次数（opt.maxiter）或减小容差（opt.tol）。

参考文献：此函数基于 ARPACK 包，由 R. Lehoucq、K. Maschhoff、D. Sorensen 和 C. Yang 编写。更多信息请参见 http://www.caam.rice.edu/software/ARPACK/。

另请参阅： eig，svds。

s = svds (A) ¶

s = svds (A, k) ¶

s = svds (A, k, sigma) ¶

s = svds (A, k, sigma, opts) ¶

[u, s, v] = svds (…) ¶

[u, s, v, flag] = svds (…) ¶

求矩阵 A 的几个奇异值。

奇异值的计算方式如下：

[m, n] = size (A);
s = eigs ([sparse(m, m), A;
                     A', sparse(n, n)])

eigs 返回的特征值对应于 A 的奇异值。要计算的奇异值数量由 k 给出，默认为 6。

参数 sigma 指定要查找的奇异值。当 sigma 是字符串 'L' 时（默认情况），查找 A 的最大奇异值。否则，sigma 必须是一个实数标量，查找最接近 sigma 的奇异值。作为推论，sigma = 0 查找最小的奇异值。请注意，对于较小的 sigma 值，可能无法找到所需数量的奇异值。在这种情况下，应增大 sigma。

opts 是一个定义了 svds 将传递给 eigs 的选项的结构体。该结构体的可能字段在 eigs 中有说明。默认情况下，svds 设置以下三个字段：

tol: 奇异值所需的收敛容差。默认值为 1e-10。传递给 eigs 的是 tol / sqrt (2)。
maxit: 最大迭代次数。默认值为 300。
disp: 诊断打印输出的级别（0|1|2）。如果 disp 为 0，则禁用诊断输出。默认值为 0。

如果请求了多个输出，则 svds 将返回 A 的奇异值分解的近似值：

A_approx = u*s*v'

其中 A_approx 是一个大小与 A 相同但秩为 k 的矩阵。

flag 在算法成功收敛时返回 0，否则返回 1。收敛性测试如下：

norm (A*v - u*s, 1) <= tol * norm (A, 1)

svds 最适合从一个大型稀疏矩阵中仅查找几个奇异值。否则，svd (full (A)) 通常会更高效率。

另请参阅： svd，eigs。

脚注

(10)

CHOLMOD、UMFPACK 和 CXSPARSE 包由 Tim Davis 编写，可在 http://faculty.cse.tamu.edu/davis/suitesparse.html 获取。