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) は、2 つの出力形式で並べ替えられた下三角行列 L のみを返します。ブロックの対角要素 D のような置換情報は失われます。既定の設定では、関数 ldlA の対角要素と下三角行列を参照し、上三角行列が下三角行列の複素共役転置であると推定します。そのため [L,D,P] = ldl(TRIL(A))[L,D,P] = ldl(A) は、両方とも正確に同じ要素を返します。この構文は、スパース A に対して有効ではありません。

[L,D] = ldl(A) には、A = L*D*L' となるようなブロック対角行列 D と並べ替えられた下三角行列 L が格納されます。ブロックの対角行列 D は、対角要素に 1 行 1 列と 2 行 2 列のブロックをもちます。この構文は、スパース A に対して有効ではありません。

[L,D,P] = ldl(A) は、P'*A*P = L*D*L' となる単位下三角行列 L、ブロック対角要素 D、および置換行列 P を返します。これは、[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 の対角要素と上三角行列のみを参照し、下三角行列が上三角行列の複素共役転置であると推定します。この構文は、P'*A*P = U'*D*U となる上三角行列 U を返します (A は Hermitian であり、上三角行列ではありません)。同様に [L,D,P] = ldl(A,'lower') は既定の動作を行います。

[U,D,p] = ldl(A,'upper','vector') は、[L,D,p] = ldl(A,'lower','vector') のように、ベクトルとして p の置換情報を返します。A は、非スパース行列でなければなりません。

[L,D,P,S] = ldl(A) は、P'*S*A*S*P = L*D*L' となる単位下三角行列 L、ブロック対角要素 D、置換行列 P、およびスケーリング行列 S を返します。この構文は、実スパース行列に対してのみ利用可能で、A の下三角行列のみを参照します。

[L,D,P,S] = LDL(A,THRESH) は、アルゴリズムのピボットの許容誤差として THRESH を使用します。THRESH は、区間 [0, 0.5] 内の double のスカラーでなければなりません。THRESH の既定値は 0.01 です。THRESH をより小さな値にすると、分解時間が短くなり、要素数が少なくなる可能性もありますが、安定しない場合もあります。この構文は、実スパース行列に対してのみ利用可能です。

[U,D,p,S] = LDL(A,THRESH,'upper','vector') は、上記のように、ピボット許容誤差を設定し、上三角行列 U と置換ベクトル p を返します。

以下の例は、関数 ldl のさまざまな形式の使用方法を説明します。1 から 3 までの出力形式と、vectorupper のオプションを使用します。トピックは以下になります。

これらの例を実行する前に、以下の正定値行列と Hermitian 不定値行列を作成する必要があります。

A = full(delsq(numgrid('L', 10)));
B = gallery('uniformdata',10,0);
M = [eye(10) B; B' zeros(10)]; 

ここで 構造体 M は、最適化と流体の問題で非常に一般的なもので、実際に M は不定です。関数 ldl はスパース引数を受け入れないため、正定値行列 A は非スパース行列であることに注意してください。

例 1 — ldl の 2 出力形式

関数 ldl の 2 出力形式は、A-(L*D*L') が小さくなるような LD を返します。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)

ab を与え、LADA を使用して 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 の 3 出力形式

3 出力形式も同様に置換行列を返すため、事実上 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 を与え、LmDmPm を使用して、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, 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.

拡張機能

参考

| |