equilibrate

缩放矩阵以改善条件

语法

[P,R,C] = equilibrate(A)

[P,R,C] = equilibrate(A,outputForm)

说明

[P,R,C] = equilibrate(A) 置换和重新缩放矩阵 A，使新矩阵 B = R*P*A*C 的对角线上元素的模为 1，非对角线上元素的模不大于 1。

[P,R,C] = equilibrate(A,outputForm) 以 outputForm 指定的形式返回输出 P、R 和 C。例如，您可以将 outputForm 指定为 "vector" 以便以列向量形式返回输出。

示例

平衡矩阵以改善条件

平衡具有较大的条件数的矩阵，以提高用迭代求解器 gmres 求解线性方程组的效率和稳定性。

加载 west0479 矩阵，这是一个实数值 479×479 稀疏矩阵。使用 condest 计算矩阵的估计条件数。

load west0479
A = west0479;
c1 = condest(A)

c1 = 1.4244e+12

尝试使用 gmres 求解线性方程组 $\mathrm{Ax}=\mathit{b}$，迭代 450 次，容差为 1e-11。指定五个输出，以便 gmres 在每次迭代中返回解的残差范数（使用 ~ 隐藏不需要的输出）。在半对数图中绘制残差范数。该图显示，gmres 无法实现合理的残差范数，因此对 $\mathit{x}$ 计算得出的解不可靠。

b = ones(size(A,1),1);
tol = 1e-11;
maxit = 450;
[x,flx,~,~,rvx] = gmres(A,b,[],tol,maxit);
semilogy(rvx)
title('Residual Norm at Each Iteration')

使用 equilibrate 置换和重新缩放 A。创建新矩阵 B = R*P*A*C，它具有更好的条件数且对角线上元素只有 1 和 -1。

[P,R,C] = equilibrate(A);
B = R*P*A*C;
c2 = condest(B)

c2 = 5.1036e+04

使用 equilibrate 返回的输出，您可以将问题 $\mathrm{Ax}=\mathit{b}$ 重新表示为 $\mathrm{By}=\mathit{d}$，其中 $\mathit{B}=\mathrm{RPAC}$ 且 $\mathit{d}=\mathrm{RPb}$。在此形式中，您可以使用 $\mathit{x}=\mathrm{Cy}$ 还原原始方程组的解。然而，还原的解可能不具有原始方程组所需的残差。有关详细信息，请参阅重新缩放以求解线性方程组。

使用 gmres 求解 $\mathrm{By}=\mathit{d}$ 以得出 $\mathit{y}$，然后重新绘制每次迭代的残差范数。该图显示，对矩阵进行平衡处理提高了问题的稳定性，gmres 在不到 200 次迭代中收敛于期望容差 1e-11。

d = R*P*b;
[y,fly,~,~,rvy] = gmres(B,d,[],tol,maxit);
hold on
semilogy(rvy)
legend('Original', 'Equilibrated', 'Location', 'southeast')
title('Relative Residual Norms (No Preconditioner)')
hold off

使用预条件子改进解

获得矩阵 B 后，要进一步提高问题的稳定性，您可以计算预条件子并将其与 gmres 结合使用。B 的数值属性优于原始矩阵 A 的数值属性，因此您应该使用经平衡处理的矩阵来计算预条件子。

用 ilu 算出两个不同的预条件子，并将它们用作 gmres 的输入来再次求解问题。在绘制了平衡后范数的同一张图上，叠加绘制应用预条件子后每次迭代的残差范数，以进行比较。该图显示，基于经平衡处理的矩阵算出的预条件子极大地提高了问题的稳定性，gmres 在不到 30 次迭代中达到了所需的容差。

semilogy(rvy)
hold on

[L1,U1] = ilu(B,struct('type','ilutp','droptol',1e-1,'thresh',0));
[yp1,flyp1,~,~,rvyp1] = gmres(B,d,[],tol,maxit,L1,U1);
semilogy(rvyp1)

[L2,U2] = ilu(B,struct('type','ilutp','droptol',1e-2,'thresh',0));
[yp2,flyp2,~,~,rvyp2] = gmres(B,d,[],tol,maxit,L2,U2);
semilogy(rvyp2)

legend('No preconditioner', 'ILUTP(1e-1)', 'ILUTP(1e-2)')
title('Relative Residual Norms with ILU Preconditioner (Equilibrated)')
hold off

指定输出格式

创建一个 6×6 幻方矩阵，然后使用 equilibrate 对该矩阵进行分解。默认情况下，equilibrate 以矩阵形式返回置换因子和缩放因子。

A = magic(6)

A = 6×6

    35     1     6    26    19    24
     3    32     7    21    23    25
    31     9     2    22    27    20
     8    28    33    17    10    15
    30     5    34    12    14    16
     4    36    29    13    18    11

[P,R,C] = equilibrate(A)

P = 6×6

     0     0     0     0     1     0
     0     0     0     0     0     1
     0     0     0     1     0     0
     1     0     0     0     0     0
     0     0     1     0     0     0
     0     1     0     0     0     0

R = 6×6

    0.1852         0         0         0         0         0
         0    0.1749         0         0         0         0
         0         0    0.1909         0         0         0
         0         0         0    0.1588         0         0
         0         0         0         0    0.1793         0
         0         0         0         0         0    0.1966

C = 6×6

    0.1799         0         0         0         0         0
         0    0.1588         0         0         0         0
         0         0    0.1588         0         0         0
         0         0         0    0.2422         0         0
         0         0         0         0    0.2066         0
         0         0         0         0         0    0.2035

使用矩阵因子 P、R 和 C 构造经平衡处理的矩阵 Bmatrix。

Bmatrix = R*P*A*C

Bmatrix = 6×6

    1.0000    0.1471    1.0000    0.5385    0.5358    0.6031
    0.1259    1.0000    0.8056    0.5509    0.6506    0.3916
    0.2747    0.8485    1.0000    0.7859    0.3943    0.5825
    1.0000    0.0252    0.1513    1.0000    0.6233    0.7754
    1.0000    0.2562    0.0569    0.9553    1.0000    0.7295
    0.1061    0.9988    0.2185    1.0000    0.9341    1.0000

现在，指定 "vector" 选项以便以向量形式返回输出。对于大型输入矩阵，以向量形式返回输出可以节省内存并提高效率。

[p,r,c] = equilibrate(A,"vector")

p = 6×1

     5
     6
     4
     1
     3
     2

r = 6×1

    0.1852
    0.1749
    0.1909
    0.1588
    0.1793
    0.1966

c = 6×1

    0.1799
    0.1588
    0.1588
    0.2422
    0.2066
    0.2035

使用向量输出 p、r 和 c 构造经平衡处理的矩阵 Bvector。

Bvector = r.*A(p,:).*c'

Bvector = 6×6

    1.0000    0.1471    1.0000    0.5385    0.5358    0.6031
    0.1259    1.0000    0.8056    0.5509    0.6506    0.3916
    0.2747    0.8485    1.0000    0.7859    0.3943    0.5825
    1.0000    0.0252    0.1513    1.0000    0.6233    0.7754
    1.0000    0.2562    0.0569    0.9553    1.0000    0.7295
    0.1061    0.9988    0.2185    1.0000    0.9341    1.0000

比较 Bvector 和 Bmatrix。结果表明它们是相等的。

norm(Bmatrix - Bvector)

ans = 0

输入参数

A — 输入矩阵

方阵

输入矩阵，指定为方阵。A 可以是稠密矩阵，也可以是稀疏矩阵，但在结构上必须为非奇异矩阵，这由 sprank 决定。

当 A 为稀疏矩阵时，输出 P、R 和 C 也为稀疏矩阵。

数据类型: single | double
复数支持: 是

outputForm — 输出格式

"matrix" (默认) | "vector"

输出格式，指定为 "vector" 或 "matrix"。此选项允许您指定输出 P、R 和 C 是以列向量还是矩阵形式返回：

• "matrix" - P、R 和 C 是矩阵，满足 B = R*P*A*C。

• "vector" - P、R 和 C 是列向量，满足 B = R.*A(P,:).*C'。

示例: [P,R,C] = equilibrate(A,"vector") 以向量形式返回 P、R 和 C。

输出参量

P — 置换信息

矩阵 (默认) | 向量

置换信息，以矩阵或向量形式返回。P*A（或 A(P,:) 带 "vector" 选项）是 A 的置换，它使其对角线元素乘积的绝对值最大化。

R — 行缩放

对角矩阵 (默认) | 向量

行缩放，以对角矩阵或向量形式返回。R 和 C 中的非零项是正实数。

C — 列缩放

对角矩阵 (默认) | 向量

列缩放，以对角矩阵或向量形式返回。R 和 C 中的非零项是正实数。

详细信息

重新缩放以求解线性方程组

对于线性方程组解 x = A\b，A 的条件数对于计算的准确度和效率很重要。equilibrate 可以通过重新缩放基向量来改进 A 的条件数。这实际上是形成一个新坐标系以便同时在其中表示 b 和 x。

当原始方程组 x = A\b 中的 b 和 x 向量的缩放不相关时，equilibrate 最有用。但是，如果 b 和 x 的缩放是相关的，则不推荐使用 equilibrate 来重新缩放 A 以仅用于线性方程组求解。所获得的解通常不会产生原始方程组的小残差，即使使用原始基向量来表示也是如此。

参考

[1] Duff, I. S., and J. Koster. “On Algorithms For Permuting Large Entries to the Diagonal of a Sparse Matrix.” SIAM Journal on Matrix Analysis and Applications 22, no. 4 (January 2001): 973–96.

扩展功能

基于线程的环境

使用 MATLAB® backgroundPool 在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool 加快代码运行速度。

此函数完全支持基于线程的环境。有关详细信息，请参阅在基于线程的环境中运行 MATLAB 函数。

版本历史记录

在 R2019a 中推出

R2022a: 指定输出格式

将 outputForm 指定为 "vector" 或 "matrix"，以控制 equilibrate 是以向量还是矩阵形式返回输出参量。对于大型分解，将输出作为向量返回可以节省内存并提高效率。

另请参阅

matchpairs | sprank | ilu | ichol | balance | gmres