lu

行列の LU 分解

構文

[L,U] = lu(A)

[L,U,P] = lu(A)

[L,U,P] = lu(A,outputForm)

[L,U,P,Q] = lu(S)

[L,U,P,Q,D] = lu(S)

[___] = lu(S,thresh)

[___] = lu(___,outputForm)

説明

[L,U] = lu(A) は非スパース行列またはスパース行列 A を、上三角行列 U と、置換された下三角行列 L に因数分解します。このとき A = L*U が成立します。

[L,U,P] = lu(A) はさらに、A = P'*L*U を満たす、置換行列 P を返します。この構文では、L は単位下三角行列で、U は上三角行列です。

[L,U,P] = lu(A,outputForm) は、outputForm で指定した形式で P を返します。outputForm を 'vector' として指定すると、P は A(P,:) = L*U を満たす置換ベクトルとして返されます。

[L,U,P,Q] = lu(S) は、スパース行列 S を、単位下三角行列 L、上三角行列 U、行置換行列 P、列置換行列 Q に因数分解します。このとき、P*S*Q = L*U が成立します。

[L,U,P,Q,D] = lu(S) はさらに、P*(D\S)*Q = L*U を満たす対角スケーリング行列 D も返します。一般的に、行のスケーリングは、よりスパースで安定した因数分解が得られます。

[___] = lu(S,thresh) は、前述の出力引数の任意の組み合わせを使用して、lu により適用されるピボットの手法のしきい値を指定します。指定した出力引数の数によって、thresh 入力の既定値と要件は異なります。詳細については、thresh 引数の説明を参照してください。

[___] = lu(___,outputForm) は、outputForm で指定した形式で P と Q を返します。P と Q を置換ベクトルとして返すには、outputForm を 'vector' として指定します。前述の構文にある任意の入力引数の組み合わせが使用できます。

例

行列の LU 分解

行列の LU 分解を計算し、生成される因子を確認します。LU 分解は、行列 $\mathit{A}$ を、上三角行列 $\mathit{U}$、下三角行列 $\mathit{L}$、置換行列 $\mathit{P}$ に分解する手法です。このとき、$\mathrm{PA}=\mathrm{LU}$ が成立します。これらの行列は、行列が行簡約階段形になるまでガウスの消去法を実行するために必要なステップを表します。行列 $\mathit{L}$ にはすべての乗算器が含まれます。また、置換行列 $\mathit{P}$ では行交換が考慮されます。

3 行 3 列の行列を作成し、LU 因子を計算します。

A = [10 -7 0
     -3  2 6
      5 -1 5];

[L,U] = lu(A)

L = 3×3

    1.0000         0         0
   -0.3000   -0.0400    1.0000
    0.5000    1.0000         0

U = 3×3

   10.0000   -7.0000         0
         0    2.5000    5.0000
         0         0    6.2000

これらの因子を乗算して、A を再作成します。lu は 2 入力の構文を使用して、置換行列 P を直接 L 因子に組み込みます。このとき、返される L は実際には P'*L となり、そのため A = L*U となります。

L*U

ans = 3×3

    10    -7     0
    -3     2     6
     5    -1     5

3 つの出力を指定して、置換行列を L の乗算器から分離することができます。

[L,U,P] = lu(A)

L = 3×3

    1.0000         0         0
    0.5000    1.0000         0
   -0.3000   -0.0400    1.0000

U = 3×3

   10.0000   -7.0000         0
         0    2.5000    5.0000
         0         0    6.2000

P = 3×3

     1     0     0
     0     0     1
     0     1     0

P'*L*U

ans = 3×3

    10    -7     0
    -3     2     6
     5    -1     5

LU 分解による線形システムの求解

LU 分解を実行し、因子を使用して問題を単純化することにより、線形システムを解きます。その結果を、バックスラッシュ演算子や decomposition オブジェクトを使用した他の方法と比較します。

5 行 5 列の魔方陣行列を作成し、線形システム $\mathrm{Ax}=\mathit{b}$ を解きます。ここで、b のすべての要素が 65 (魔方陣の和) と等しくなります。65 はこの行列の魔方陣の和である (すべての行および列ではその和が 65 になる) ため、x について予期される解は 1 から成るベクトルです。

A = magic(5);
b = 65*ones(5,1);
x = A\b

x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

一般的な正方行列では、バックスラッシュ演算子は LU 分解を使用して線形システムの解を計算します。LU 分解により、A は三角行列の積として表され、三角行列を含む線形システムは代入式により簡単に解けます。

バックスラッシュにより計算された解を再作成するには、A の LU 分解を計算します。その後、因子を使用して次の 2 つの三角線形システムを解きます。

y = L\(P*b);
x = U\y;

線形システムを解く前に行列因子を事前計算するこの方法は、多数の線形システムを解く場合のパフォーマンスを向上させることができます。これは、因数分解の発生が 1 回のみで、繰り返しが不要なためです。

[L,U,P] = lu(A)

L = 5×5

    1.0000         0         0         0         0
    0.7391    1.0000         0         0         0
    0.4783    0.7687    1.0000         0         0
    0.1739    0.2527    0.5164    1.0000         0
    0.4348    0.4839    0.7231    0.9231    1.0000

U = 5×5

   23.0000    5.0000    7.0000   14.0000   16.0000
         0   20.3043   -4.1739   -2.3478    3.1739
         0         0   24.8608   -2.8908   -1.0921
         0         0         0   19.6512   18.9793
         0         0         0         0  -22.2222

P = 5×5

     0     1     0     0     0
     1     0     0     0     0
     0     0     0     0     1
     0     0     1     0     0
     0     0     0     1     0

y = L\(P*b);
x = U\y

x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

decomposition オブジェクトも、特殊化された因数分解により線形システムを解く場合に便利です。この方法では、因子の使用方法の知識を求められることなく、行列因子の事前計算によるパフォーマンス上のメリットを数多く得られます。同じ結果を再作成するには、タイプが 'lu' の decomposition オブジェクトを使用します。

dA = decomposition(A,'lu');
x = dA\b

x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

スパース行列の LU 分解

スパース行列の LU 分解を計算して、恒等式 L*U = P*S*Q を検証します。

Buckminster-Fuller 幾何ドームの結合グラフの 60 行 60 列のスパース隣接行列を作成します。

S = bucky;

4 つの出力をもつスパース行列構文を使用して S の LU 分解を計算し、行と列の置換行列を返します。

[L,U,P,Q] = lu(S);

S の行と列を P*S*Q と置換した結果を、三角因子 L*U を乗算した場合と比較します。1 ノルムの違いは丸め誤差の範囲に入るため、L*U = P*S*Q となります。

e = P*S*Q - L*U;
norm(e,1)

ans = 2.2204e-16

置換ベクトルを使用してメモリを節約

行列の LU 分解を計算します。行置換を行列ではなくベクトルとして返すことで、メモリを節約します。

1000 行 1000 列のランダム行列を作成します。

A = rand(1000);

行列 P として保存されている置換情報を使用して LU 分解を計算します。結果を、ベクトル p として保存されている置換情報と比較します。行列が大きいほど、置換ベクトルを使用した場合にメモリ効率がより高くなります。

[L1,U1,P] = lu(A);
[L2,U2,p] = lu(A,'vector');
whos P p

  Name         Size                Bytes  Class     Attributes

  P         1000x1000            8000000  double              
  p            1x1000               8000  double              

また、置換ベクトルを使用すると、以降の操作の実行時間が短縮されます。たとえば、前述の LU 分解を使用して、線形システム $\mathrm{Ax}=\mathit{b}$ を解くことができます。置換ベクトルから得られる解と置換行列から得られる解は等価 (丸めを除く) ですが、通常は、置換ベクトルを使用して解を求める方が時間がやや短縮されます。

スパース行列の因数分解の非スパース要素を減らす

スパース行列の LU 分解の計算結果を、列置換がある場合と列置換がない場合で比較します。

行列 west0479 を読み込みます。これは 479 行 479 列の実数値のスパース行列です。

load west0479
A = west0479;

3 つの出力を指定して lu を呼び出し、A の LU 分解を計算します。L 因子と U 因子の spy プロットを生成します。

[L,U,P] = lu(A);
subplot(1,2,1)
spy(L)
title('L factor')
subplot(1,2,2)
spy(U)
title('U factor')

次に、4 つの出力を指定して lu を使用し、A の LU 分解を計算します。これにより、A の列が置換され、各因子について非ゼロの数が少なくなります。この結果生成される因子は、列置換を使用しない場合と比べて、スパースの傾向が強くなっています。

[L,U,P,Q] = lu(A);
subplot(1,2,1)
spy(L)
title('L factor')
subplot(1,2,2)
spy(U)
title('U factor')

入力引数

A — 入力行列

行列

入力行列。A は非スパースまたはスパースで、正方形または長方形のサイズにすることができます。

データ型: single | double
複素数のサポート: あり

S — スパース入力行列

スパース行列

スパース入力行列。S は正方形または長方形のサイズにすることができます。

データ型: double
複素数のサポート: あり

thresh — スパース行列に対するピボットしきい値

スカラー | 2 要素ベクトル

スパース行列に対するピボットしきい値。スカラーまたは 2 要素ベクトルとして指定します。有効な値は、[0 1] 間の値です。thresh を指定する方法は、lu の呼び出しで指定される出力の数によって異なります。

• 出力の数が 3 つ以下の場合、thresh はスカラーでなければならず、既定値は 1.0 です。

• 出力の数が 4 つ以上の場合、thresh はスカラーまたは 2 要素ベクトルにできます。既定値は [0.1 0.001] です。thresh をスカラーとして指定した場合、ベクトルの最初の値のみが置き換えられます。

上位レベルでは、この入力により精度と合計実行時間がトレードオフされます。thresh の値を小さくするほど、スパース LU 分解を導く傾向がありますが、解は非対称になる可能性があります。より大きな値は、より正確な解を導くことができ (常にではありませんが)、通常、作業とメモリ使用量の合計が増加します。

lu は、ピボットの手法を、まず出力引数の数、次に因数分解する行列のプロパティに基づいて選択します。いずれの場合でも、しきい値を 1.0 に設定すると部分ピボットが行われ、0 に設定すると、生成される行列のスパース性のみに基づいてピボットが選択されます。L のすべての値は、絶対値が 1/min(thresh) 以下です。

• 出力引数が 3 つ以下 — アルゴリズムは以下の方程式を満たす場合に対角ピボットを選択します。

  A(j,j) >= thresh * max(abs(A(j:m,j)))

  それ以外の場合、絶対値が最大である要素を含む行を選択します。

• 対称なピボットの手法 — S が、ほぼ対称構造で対角の大半が非ゼロの正方スパース行列である場合、lu は対称なピボットの手法を使用します。この手法では、アルゴリズムは以下の不等式を満たす場合に対角ピボット j を選択します。

  A(i,j) >= thresh(2) * max(abs(A(j:m,j)))

  対角要素がこのテストに失敗した場合、lu は以下の不等式を満たす最もスパースな行 i を選択します。

  A(i,j) >= thresh(1) * max(abs(A(j:m,j)))

• 非対称なピボットの手法 — S が対称なピボットの手法の要件を満たさない場合、lu は非対称な手法を使用します。この場合、lu は以下の不等式を満たす最もスパースな行 i を選択します。

  A(i,j) >= thresh(1) * max(abs(A(j:m,j)))

  thresh(1) の値が 1.0 の場合、通常の部分ピボットを行います。L この要素は、1/thresh(1) の絶対値またはそれより小さい値をもちます。非対称の手法では、入力ベクトル thresh の 2 番目の要素は使用されません。

メモ

非常にまれではありますが、誤った因数分解により P*S*Q ≠ L*U となることがあります。このような場合は、thresh を最大 1.0 (通常の部分ピボット) まで大きくして、再度試みてください。

outputForm — 置換出力の形状

'matrix' (既定値) | 'vector'

置換出力の形状。'matrix' または 'vector' として指定します。このフラグは、lu が行置換 P と列置換 Q を置換行列として返すか置換ベクトルとして返すかを制御します。

行列の場合、P と Q は以下の恒等式を満たします。

• 出力が 3 つ — P は P*A = L*U を満たします。

• 出力が 4 つ — P と Q は P*S*Q = L*U を満たします。

• 出力が 5 つ — P、Q、D は P*(D\S)*Q = L*U を満たします。

ベクトルの場合、出力 P と Q は以下の恒等式を満たします。

• 出力が 3 つ — P は A(P,:) = L*U を満たします。

• 出力が 4 つ — P と Q は S(P,Q) = L*U を満たします。

• 出力が 5 つ — P、Q、D は D(:,P)\S(:,Q) = L*U を満たします。

例: [L,U,P] = lu(A,'vector')

出力引数

L — 下三角因子

行列

下三角因子。行列として返されます。L の形式は、行置換 P が個別の出力で返されるかどうかによって異なります。

• 3 番目の出力 P を指定した場合、L は単位下三角行列 (つまり、主対角が 1 の下三角行列) として返されます。

• 3 番目の出力 P を指定しない場合、L は単位下三角行列の行置換として返されます。具体的には、これは出力が 3 つの場合に返される出力 P と L の積 P'*L です。

U — 上三角因子

行列

上三角因子。上三角行列として返されます。

P — 行置換

ベクトル | 行列

行置換。置換行列、または置換ベクトル ('vector' オプションが指定されている場合) として返されます。この出力を使用すると、計算の数値安定性が向上します。

この出力が満たす恒等式の詳細については、outputForm を参照してください。

Q — 列置換

ベクトル | 行列

列置換。置換行列、または置換ベクトル ('vector' オプションが指定されている場合) として返されます。この出力を使用すると、スパース行列の因子の非スパース要素 (非ゼロの数) が減ります。

この出力が満たす恒等式の詳細については、outputForm を参照してください。

D — 行のスケーリング

対角行列

行のスケーリング。対角行列として返されます。D は、P*(D\S)*Q = L*U を満たすように S をスケーリングするために使用されます。常にではありませんが、一般的に、行のスケーリングを使用すると、よりスパースで安定した因数分解が得られます。

アルゴリズム

LU 分解は、ガウスの消去法の変形を使用して計算されます。正確な解を計算できるかどうかは、元の行列の条件数の値 cond(A) によって異なります。行列の条件数が大きい (ほぼ特異である) 場合、計算された因数分解は正確でないことがあります。

LU 分解は、inv を使用して逆行列を求めたり、det を使用して行列式の値を求めるための重要な部分です。線形方程式の求解や、演算子 \ や / を使用した行列の除算に対しても基本的な役割を果たします。そのため、必然的に、これらの依存関数にも lu の数値制限が存在することになります。

参照

[1] Gilbert, John R., and Tim Peierls. “Sparse Partial Pivoting in Time Proportional to Arithmetic Operations.” SIAM Journal on Scientific and Statistical Computing 9, no. 5 (September 1988): 862–874. https://doi.org/10.1137/0909058.

[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] Davis, Timothy A. "Algorithm 832: UMFPACK V4.3 – an unsymmetric-pattern multifrontal method." ACM Transactions on Mathematical Software 30, no. 2 (June 2004): 196–199. https://doi.org/10.1145/992200.992206.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意事項および制限事項:

• コード生成では、この関数のスパース行列入力はサポートされません。

GPU コード生成

GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

使用上の注意事項および制限事項:

• コード生成では、この関数のスパース行列入力はサポートされません。

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

この関数はスレッドベースの環境を完全にサポートしています。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

使用上の注意事項および制限事項:

• A は非スパースでなければなりません。

詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

分散配列

Parallel Computing Toolbox™ を使用して、クラスターの結合メモリ上で大きなアレイを分割します。

使用上の注意事項および制限事項:

• A は非スパースでなければなりません。

詳細については、分散配列を使用した MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

参考

cond | det | inv | ilu | qr | rref | mldivide | decomposition

トピック

• 分解