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のような置換情報は失われます。既定の設定では、関数ldlはAの対角要素と下三角行列を参照し、上三角行列が下三角行列の複素共役転置であると推定します。そのため[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 までの出力形式と、vectorと上のオプションを使用します。トピックは以下になります。

これらの例を実行する前に、以下の正定値行列と 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')が小さくなるようなLとDを返します。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 の 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を与え、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.