LU分解シミュレーター 戻る
数値線形代数シミュレーター

LU分解シミュレーター — Ax=bの直接解法

3×3行列をL·Uに分解し、前進代入Ly=bと後退代入Ux=yでAx=bを解く過程を可視化。対角要素と右辺を変えて、行列式・解・残差・条件数の変化を学べます。

パラメータ設定
A[1,1](対角)
A[2,2](対角)
A[3,3](対角)
b[1](右辺)

非対角要素は固定: A[1,2]=3, A[1,3]=2, A[2,1]=1, A[2,3]=1, A[3,1]=2, A[3,2]=1。b[2]=8, b[3]=13 固定。

計算結果
行列式 det(A)
解 x_1
解 x_2
残差 ‖Ax-b‖_2
行列 A と右辺ベクトル b

左=係数行列A/右=右辺ベクトルb(青セル=非対角、緑セル=対角・可変、赤縁=ピボット行交換あり)

L·U 分解と前進・後退代入の解

緑=下三角L(対角=1)/青=上三角U/黄=中間ベクトルy=L⁻¹·b/橙=解x=U⁻¹·y

理論・主要公式

LU分解は、正方行列AをL(単位下三角)とU(上三角)の積に分解する手法です。部分ピボット選択を含めると P·A = L·U と書けます。

分解の式(Doolittle法、Lの対角が1):

$$P \cdot A = L \cdot U, \quad L_{ii} = 1$$

Ax=b は2段階の三角系で解けます。前進代入で y を、後退代入で x を求めます:

$$L\,y = P\,b \quad\Longrightarrow\quad U\,x = y$$

行列式と残差。det(A) はUの対角の積(行交換数の符号付き)、残差は数値誤差の指標:

$$\det(A) = (-1)^{s}\prod_i U_{ii}, \qquad r = \|A\,x - b\|_2$$

計算量は O(n³)。一度分解すれば同じ A に対し複数の右辺 b を O(n²) で解けるため、剛性行列を共有するFEMの多荷重ケースで威力を発揮します。

LU分解シミュレーターとは

🙋
連立一次方程式って、中学のときに加減法で解いたあれですよね?コンピュータも同じことやってるんですか?
🎓
基本は同じ「ガウス消去」だよ。ただ実務では、その消去の途中経過を「LとUの2つの行列」として残しておく流儀がある。これがLU分解だ。上のシミュレーターの2段目を見てごらん。緑のLが「消去のときに各行に掛けた倍率」、青のUが「上三角に整理し終わった姿」。この2つを組み合わせるとAに戻る。
🙋
わざわざ分解するメリットがあるんですか?1回解ければ十分な気もしますけど。
🎓
そこがミソだ。例えばFEMの線形静解析だと、同じ剛性行列Kに対して荷重ケースを10通り、100通り解きたいことがある。LUを一度作っておけば、各ケースは「前進代入と後退代入」だけで済むんだ。1回の分解はO(n³)で重いけど、その後の各解はO(n²)で安い。1000元の方程式なら、最初の分解1回ぶんのコストで1000ケースぐらい解けてしまう。
🙋
「部分ピボット選択」って何ですか?シミュレーターの下にちらっと書いてある。
🎓
消去のときに割り算をする「対角要素(ピボット)」が0だと破綻するし、0に近い小さい値だと誤差が爆発する。だから各段階で、いまの列の中から絶対値が一番大きい要素を選んで行を入れ替えてから消去する。これが部分ピボット選択だ。シミュレーターでA[1,1]を1にしてb[1]を-20にしてみて。行交換が走って解は変わらず計算できるはずだよ。
🙋
残差のカードがほぼ0なのは、計算が正しい証拠ってことですか?
🎓
そう。求めたxを元のAに掛け直して、bと比べた誤差ベクトルのノルムが残差だ。理論的には0、実機の倍精度浮動小数点でも10⁻¹⁰オーダーにしかならない。もしここが1とか10になったら、それは行列がほぼ特異(条件数が極端に大きい)か、実装にバグがあるかのどちらかだ。残差のチェックは「数値解の自己検算」として実務でも必ずやる癖をつけよう。

よくある質問

どちらもLU分解の実装方式で、Lの対角を1とするのがDoolittle法、Uの対角を1とするのがCrout法です。本ツールはDoolittle法を使っており、Lの対角はすべて1、Uの対角がピボット値になります。数値的な結果は本質的に同じで、メモリ配置や演算順序の好みで使い分けられます。教科書ではDoolittleの方が広く使われています。
条件数は cond(A)=‖A‖·‖A⁻¹‖ で定義され、入力データの誤差がどれだけ解に増幅されるかの指標です。条件数が10ᵏ程度のとき、おおむね倍精度浮動小数点(約16桁の精度)の解はk桁失われると見積もれます。条件数が10¹⁶を超えると有効桁数が完全に消え、解は信用できません。実務では反復改良法(iterative refinement)やスケーリング、別の定式化への切り替えで対処します。
大規模FEMやCFDの行列は要素の99%以上が0の「疎行列」になります。素直にLU分解すると、消去の途中で0だった位置に非零要素が生まれる「fill-in」が起きてメモリが破綻します。実務ではAMD・METIS等の順序付けで行・列を並び替えてfill-inを最小化し、SuperLU・PARDISO・MUMPSといった疎LU分解ライブラリを使います。それでも巨大すぎる場合は、共役勾配法(CG)やGMRESなどの反復解法に切り替えます。
線形弾性解析の剛性行列Kは対称正定値になることが多く、その場合はLU分解より約半分の計算量とメモリで済むCholesky分解(K=L·Lᵀ)が第一選択です。接触解析や塑性解析、非対称な定式化(流体構造連成等)では一般のLU分解を使います。動解析の質量行列も対称正定値なので、固有値解析でも Cholesky が活躍します。いずれも「一度分解して複数の右辺を解く」場面で、直接法の真価が発揮されます。

実世界での応用

有限要素法(FEM)の線形静解析:構造解析ソルバー(NASTRAN・Abaqus・ANSYS・Marc等)の心臓部は、剛性行列に対する直接法です。対称正定値ならCholesky、一般ならLU分解を使い、同じ剛性行列に対して複数荷重ケース(自重・風荷重・地震荷重・温度荷重…)を解く場面で再利用性が効きます。スパース版(SuperLU・PARDISO・MUMPS)が業界標準で、数千万自由度のモデルでも実用時間内に解けます。

回路解析(SPICE):電子回路の節点方程式 G·v = i は、節点アドミタンス行列Gが疎な対称行列で、しばしばLU分解で解かれます。SPICE系シミュレータは時間積分の各ステップで線形化された方程式を解くため、行列パターンが変わらない区間ではLUの再利用が大きな高速化につながります。

最適化と機械学習の正規方程式:線形最小二乗 (XᵀX)·β = Xᵀy の正規方程式や、ガウス過程回帰の (K+σ²I)·α = y などは、対称正定値行列に対するCholesky分解(=特殊なLU)で解かれます。バッチ学習やベイズ推論の中核ルーチンとして、高速な行列分解の実装が性能を左右します。

制御工学のリアプノフ方程式:状態空間モデルの安定性評価で現れる連続時間リアプノフ方程式 AᵀP + PA = -Q は、Bartels-Stewart法でAをSchur分解した後、三角行列に対する代入で解きます。本質はLU分解と同じ「上三角・下三角に対する代入」の応用で、制御系設計ソフト(MATLAB Control Toolbox等)の標準実装です。

よくある誤解と注意点

最も多い誤解は、「LU分解は逆行列を計算する手段」と思い込むことです。実際には逆行列A⁻¹を陽に作る必要はほとんどなく、Ax=bを解くだけならLU分解で十分です。逆行列を作るとメモリと計算コストが約3倍になる上、数値誤差も増えます。MATLABの `inv(A)*b` ではなく `A\b`(バックスラッシュ演算子、内部でLUを使う)が推奨されるのも同じ理由です。本ツールも逆行列を一切作らず、L·yとU·xの2回の三角系解法だけで答えを得ています。

次に多いのが、部分ピボット選択を入れれば常に数値的に安定だと信じてしまうことです。部分ピボット選択は対角要素が小さくなりすぎる問題を防ぎますが、行列自体の条件数が大きい(ほぼ特異な)場合は依然として解の精度が落ちます。シミュレーターでA[1,1]=A[2,2]=A[3,3]=1のような対角が小さい組み合わせを試すと、行列式が小さくなり残差も増えがちです。実務では条件数の事前推定(rcond)や反復改良で精度を確保します。

最後に、計算量O(n³)の意味を見落として大規模問題に直接法を使う失敗が頻発します。3×3なら一瞬ですが、自由度100万のFEM行列を素直にLU分解しようとすると、メモリも時間も天文学的に膨らみます。疎行列向けのSuperLU・PARDISOでさえ、数百万自由度を超えると反復法(CG・GMRES・代数的マルチグリッド)の方が有利になることが多いです。「対称か」「正定値か」「疎か」「自由度はいくつか」を必ず確認してから解法を選びましょう。本ツールが扱うのは入門用の3×3ですが、その背後には数値線形代数の半世紀以上の蓄積があります。