Main Content

本页翻译不是最新的。点击此处可查看最新英文版本。

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 函数(包括单输出、双输出和三输出形式)以及如何使用 vectorupper 选项。涵盖的主题包括:

在运行上面任何示例之前,需要生成以下正定和不定埃尔米特矩阵:

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 返回 LD,使得 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,使用 LADAAx=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,使用 LmDmPmMx=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 可能是对角矩阵,也可能表示分块结构体。例如,通过上述的 ADA 就是对角矩阵。但如果您稍微移动 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.

扩展功能

另请参阅

| |