ldl
埃尔米特不定矩阵的分块 LDL 分解
语法
L = ldl(A)
[L,D] = ldl(A)
[L,D,P] = ldl(A)
[L,D,p] = ldl(A,'vector')
[U,D,P] = ldl(A,'upper')
[U,D,p] = ldl(A,'upper','vector')
[L,D,P,S] = ldl(A)
[L,D,P,S] = LDL(A,THRESH)
[U,D,p,S] = LDL(A,THRESH,'upper','vector')
说明
L = ldl(A)
仅返回和双输出形式中相同的置换下三角矩阵 L
。置换信息丢失,分块对角因子 D
也丢失。默认情况下,ldl
仅引用 A
的对角和下三角,并假定上三角是下三角的复共轭转置。因此 [L,D,P] = ldl(TRIL(A))
和 [L,D,P] = ldl(A)
都返回完全相同的因子。请注意,此语法对稀疏 A
无效。
[L,D] = ldl(A)
将分块对角矩阵 D
和置换下三角矩阵存储在 L
中,使得 A = L*D*L'
。分块对角矩阵 D
在其对角上具有 1×1 和 2×2 个分块。请注意,此语法对稀疏 A
无效。
[L,D,P] = ldl(A)
返回单位下三角矩阵 L
、分块对角 D
和置换矩阵 P
,使得 P'*A*P = L*D*L'
。这与 [L,D,P] = ldl(A,'matrix')
等效。
[L,D,p] = ldl(A,'vector')
将置换信息作为向量 p
而不是矩阵返回。即,p
是行向量,使得 A(p,p) = L*D*L'
。
[U,D,P] = ldl(A,'upper')
仅引用 A
的对角和上三角,并假定下三角是上三角的复共轭转置。此语法返回单位上三角矩阵 U
,使得 P'*A*P = U'*D*U
(假定 A
是埃尔米特矩阵,而不仅仅是上三角矩阵)。[L,D,P] = ldl(A,'lower')
提供相似的默认行为。
[U,D,p] = ldl(A,'upper','vector')
将置换信息作为向量 p
返回,[L,D,p] = ldl(A,'lower','vector')
也是如此。A
为满矩阵。
[L,D,P,S] = ldl(A)
返回单位下三角矩阵 L
、分块对角 D
、置换矩阵 P
和缩放矩阵 S
,使得 P'*S*A*S*P = L*D*L'
。此语法仅适用于实数稀疏矩阵,仅引用 A
的下三角。
[L,D,P,S] = LDL(A,THRESH)
将 THRESH
用作算法中的主元容差。THRESH
必须是位于区间 [0, 0.5]
内的双精度标量。THRESH
的默认值是 0.01
。使用 THRESH
的较小值可能会提供较快的分解时间和更少的条目,但也可能会导致不太稳定的分解。此语法仅适用于实数稀疏矩阵。
如上所述,[U,D,p,S] = LDL(A,THRESH,'upper','vector')
设置主元容差并返回上三角 U
和置换向量 p
。
示例
这些示例介绍如何使用不同形式的 ldl
函数(包括单输出、双输出和三输出形式)以及如何使用 vector
和 upper
选项。涵盖的主题包括:
在运行上面任何示例之前,需要生成以下正定和不定埃尔米特矩阵:
A = full(delsq(numgrid('L', 10))); B = gallery('uniformdata',10,0); M = [eye(10) B; B' zeros(10)];
此处 M
的结构体在优化和流体流动问题中很常见,M
实际上是不定矩阵。请注意,正定矩阵 A
必须为满矩阵,因为 ldl
不接受稀疏参数。
示例 1 - 双输出形式的 ldl
双输出形式的 ldl
返回 L
和 D
,使得 A-(L*D*L')
较小,L
是置换单位下三角矩阵,D
是分块 2×2 对角矩阵。另外还要注意,由于 A
是正定矩阵,D
的对角元素全为正数:
[LA,DA] = ldl(A); fprintf(1, ... 'The factorization error ||A - LA*DA*LA''|| is %g\n', ... norm(A - LA*DA*LA')); neginds = find(diag(DA) < 0)
给定 a
b
,使用 LA
、DA
对 Ax=b
求解:
bA = sum(A,2); x = LA'\(DA\(LA\bA)); fprintf(... 'The absolute error norm ||x - ones(size(bA))|| is %g\n', ... norm(x - ones(size(bA))));
示例 2 - 三输出形式的 ldl
三输出形式也返回置换矩阵,这样 L
实际上就是单位下三角矩阵:
[Lm, Dm, Pm] = ldl(M); fprintf(1, ... 'The error norm ||Pm''*M*Pm - Lm*Dm*Lm''|| is %g\n', ... norm(Pm'*M*Pm - Lm*Dm*Lm')); fprintf(1, ... 'The difference between Lm and tril(Lm) is %g\n', ... norm(Lm - tril(Lm)));
给定 b
,使用 Lm
、Dm
和 Pm
对 Mx=b
求解:
bM = sum(M,2); x = Pm*(Lm'\(Dm\(Lm\(Pm'*bM)))); fprintf(... 'The absolute error norm ||x - ones(size(b))|| is %g\n', ... norm(x - ones(size(bM))));
示例 3 - D 的结构
D
是包含 1×1 和 2×2 个分块的分块对角矩阵。这使其成为对角矩阵的一个特殊情况。如果输入矩阵为正定矩阵,D
几乎总是对角矩阵(取决于矩阵的定性)。但如果矩阵是不定矩阵,D
可能是对角矩阵,也可能表示分块结构体。例如,通过上述的 A
,DA
就是对角矩阵。但如果您稍微移动 A
,最终会得到不定矩阵,然后您可以计算包含分块结构体的 D
。
figure; spy(DA); title('Structure of D from ldl(A)'); [Las, Das] = ldl(A - 4*eye(size(A))); figure; spy(Das); title('Structure of D from ldl(A - 4*eye(size(A)))');
示例 4 - 使用 'vector' 选项
与 lu
函数一样,ldl
接受决定此函数是返回置换向量还是返回置换矩阵的参数。ldl
默认情况下返回后者。如果选择 'vector'
,此函数执行速度较快并使用较少的内存。出于此原因考虑,建议指定 'vector'
选项。需要注意的另一点是,对于此类运算,建立索引通常比相乘快:
[Lm, Dm, pm] = ldl(M, 'vector'); fprintf(1, 'The error norm ||M(pm,pm) - Lm*Dm*Lm''|| is %g\n', ... norm(M(pm,pm) - Lm*Dm*Lm')); % Solve a system with this kind of factorization. clear x; x(pm,:) = Lm'\(Dm\(Lm\(bM(pm,:)))); fprintf('The absolute error norm ||x - ones(size(b))|| is %g\n', ... norm(x - ones(size(bM))));
示例 5 - 使用 'upper' 选项
和 chol
函数一样,ldl
接受确定要引用输入矩阵的哪个三角形以及 ldl
是返回下 (L
) 三角因子还是返回上 (L'
) 三角因子的参数。对于稠密矩阵,用上三角版本代替下三角版本并不能实现真正节省:
Ml = tril(M); [Lml, Dml, Pml] = ldl(Ml, 'lower'); % 'lower' is default behavior. fprintf(1, ... 'The difference between Lml and Lm is %g\n', norm(Lml - Lm)); [Umu, Dmu, pmu] = ldl(triu(M), 'upper', 'vector'); fprintf(1, ... 'The difference between Umu and Lm'' is %g\n', norm(Umu - Lm')); % Solve a system using this factorization. clear x; x(pm,:) = Umu\(Dmu\(Umu'\(bM(pmu,:)))); fprintf(... 'The absolute error norm ||x - ones(size(b))|| is %g\n', ... norm(x - ones(size(bM))));
如果同时指定 'upper'
和 'vector'
选项,'upper'
必须位于参数列表中 'vector'
的前面。
示例 6 - linsolve 和埃尔米特不定求解器
使用 linsolve
函数时,通过利用系统具有对称矩阵的知识可能会获得更好性能。上述示例中使用的矩阵由于较小而不能看出这一点,为此在本示例中生成一个较大的矩阵。此例的矩阵是对称正定矩阵,下面我们将看到通过有关矩阵的每条信息,都存在相应的加速。即,对称求解器比普通求解器快,而对称正定求解器比对称求解器快:
Abig = full(delsq(numgrid('L', 30))); bbig = sum(Abig, 2); LSopts.POSDEF = false; LSopts.SYM = false; tic; linsolve(Abig, bbig, LSopts); toc; LSopts.SYM = true; tic; linsolve(Abig, bbig, LSopts); toc; LSopts.POSDEF = true; tic; linsolve(Abig, bbig, LSopts); toc;
参考
[1] Ashcraft, Cleve, Roger G. Grimes, and John G. Lewis. “Accurate Symmetric Indefinite Linear Equation Solvers.” SIAM Journal on Matrix Analysis and Applications 20, no. 2 (January 1998): 513–561. https://doi.org/10.1137/S0895479896296921.
[2] Anderson, E., ed. LAPACK Users’ Guide. 3rd ed. Software, Environments, Tools. Philadelphia: Society for Industrial and Applied Mathematics, 1999. https://doi.org/10.1137/1.9780898719604.
[3] Duff, Iain S. “MA57---a Code for the Solution of Sparse Symmetric Definite and Indefinite Systems.” ACM Transactions on Mathematical Software 30, no. 2 (June 2004): 118–144. https://doi.org/10.1145/992200.992202.