ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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) は、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 の下三角行列のみを参照します。ldl 関数は、実数スパース A に対して MA57 を使用します。

[L,D,P,S] = LDL(A,THRESH) は、MA57 のピボットの許容誤差として 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;

詳細

すべて折りたたむ

アルゴリズム

ldl は、Harwell サブルーチン ライブラリ (HSL) で実スパース行列に対して MA57 を使用します。

参照

[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.

参考

| |

この情報は役に立ちましたか?