matlab hermitian,Hermitian 不定矩阵的分块 LDL 分解

ldl

Hermitian 不定矩阵的分块 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 是 Hermitian 矩阵，而不仅仅是上三角矩阵)。[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 选项。涵盖的主题包括：

在运行上面任何示例之前，需要生成以下正定和不定 Hermitian 矩阵：

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 和 Hermitian 不定求解器

使用 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, C., R.G. Grimes, and J.G. Lewis.

“Accurate Symmetric Indefinite Linear Equations Solvers.” SIAM

J. Matrix Anal. Appl. Vol. 20. Number 2, 1998, pp. 513–561.

[2] Duff, I. S. "MA57 — A new code

for the solution of sparse symmetric definite and indefinite systems."

Technical Report RAL-TR-2002-024, Rutherford Appleton Laboratory,

2002.