5. 連立一次方程式： Gauss消去法とLU分解#
第3章では A \ b で解を求め、第4章ではその誤差を評価した。この章では、その中で使われる基本的な考え方を学ぶ。出発点は、高校数学でも行う「一方の式を何倍かして、もう一方の式から引く」操作である。
5.1. この章の目標#
Gauss消去法で連立一次方程式を三角形の形に変形できる。
上三角行列を後退代入で解ける。
LU分解がGauss消去の記録であることを説明できる。
行の入れ替え（ピボット選択）が必要な理由を説明できる。
n元のGauss消去を実装し、演算回数が約 \(\frac23n^3\) となる理由を説明できる。
JuliaでLU分解を使い、残差を確認できる。
産業連関分析のように、同じ行列で右辺だけを変える問題で分解を再利用できる。
5.2. まずは消去法を思い出す#
次の連立方程式を考える。
第2式から第1式の2倍を引くと
が得られる。これを第1式へ戻すと \(x_1=2\) である。行列では
のように、係数と右辺をまとめて書ける。縦線の右が右辺であり、同じ行操作を両方に施す。係数行列の対角線より下を0にしたことになる。この操作を大きな行列へ順番に行うのがGauss消去法である。
色と枠の読み方： 青い枠の \(2\) は割るときの基準（ピボット）、橙色の \(4\) は今回0にしたい成分、右の枠の \(0\) は消去した結果である。倍率は \(m=4/2=2\) となる。
0にしたいのは一つの成分だが、引き算は右辺を含む行全体に行う。次のコードでは、この三つの引き算がどの行に対応するか確認しよう。
using LinearAlgebra
A = [2.0 1.0; 4.0 3.0]
b = [5.0, 11.0]
m = A[2,1] / A[1,1] # 第1行を何倍するか
A2 = copy(A)
b2 = copy(b)
A2[2,:] .-= m .* A2[1,:]
b2[2] -= m * b2[1]
println("消去後のA = \n", A2)
println("消去後のb = ", b2)
消去後のA =
[2.0 1.0; 0.0 1.0]
消去後のb = [5.0, 1.0]
5.3. 三角行列なら順番に解ける#
消去後の行列は上三角行列である。下の式から上へ順番に未知数を求める操作を後退代入という。
なら、まず \(x_2=c_2/u_{22}\)、次に \(x_1=(c_1-u_{12}x_2)/u_{11}\) と求める。下三角行列では、上から下へ進む前進代入を使う。
function backward_substitution(U, b)
n = length(b)
x = zeros(Float64, n)
for i in n:-1:1
known = sum((U[i,j]*x[j] for j in i+1:n); init=0.0)
x[i] = (b[i] - known) / U[i,i]
end
return x
end
U = UpperTriangular(A2)
x = backward_substitution(U, b2)
println("後退代入で得た解 = ", x)
println("元の方程式の残差 = ", norm(b-A*x))
後退代入で得た解 = [2.0, 1.0]
元の方程式の残差 = 0.0
5.4. LU分解は消去操作の記録#
今の消去では「第2行から第1行の2倍を引く」と、上三角行列 \(U\) が得られた。逆に第2行へ第1行の2倍を足せば元の行列に戻る。この関係は
と書ける。 \(L\) の下側の青い枠の2は、消去で使った倍率 \(m=4/2\) の記録である。行列の積を1行だけ計算すると、元の第2行が \(2(2,1)+(0,1)=(4,3)\) と復元される。
これがLU分解の出発点である。大きな行列でも、消去で使った倍率を下三角行列 \(L\) に、消去後の係数を上三角行列 \(U\) に保存する。\(L\) の対角成分は1とする。正則な行列について、行交換なしに各段階を0で割らずに進められれば、この形の分解は存在し一意になる。正則でも途中に0のピボットが現れることがあり、その場合は後述の行交換が必要である。
\(Ax=b\) を解くとき、\(y=Ux\) と置けば、
の順で解ける。今の例では \(y_1=5\)、\(2y_1+y_2=11\) から \(y=(5,1)^T\) となる。これは先ほど消去で作った右辺と同じである。続く後退代入で \(x=(2,1)^T\) が得られる。
分解は係数行列の準備、代入は各右辺を解く作業である。この区別により、同じ行列で右辺だけが変わる計算を再利用できる。3次行列での消去係数と基本行列の関係は、後の補足で一度だけ詳しく確認する。
5.5. ピボットと行交換#
ここでも用いるのは行基本変形であり、列基本変形ではない。「第 \(k\) 列の消去」とは、各 \(i>k\) に対して
という行操作を行い、結果として第 \(k\) 列の対角成分より下を0にする、という意味である。第 \(k\) 段階の消去係数を
と計算する。ここで \(a_{ij}^{(k-1)}\) は、\(k-1\) 段階までの行操作を終えた行列の \((i,j)\) 成分を表す。実際、第 \(k\) 列の \((i,k)\) 成分は
となる。
この分母 \(a_{kk}^{(k-1)}\) が第 \(k\) 段階の ピボット であり、その位置 \((k,k)\) を ピボット位置 という。ピボットを使って、同じ列のピボットより下にある成分を0にしていく。
例えば
では、最初のピボット位置にある \(a_{11}=0\) で割れないため、そのままでは第1列を消去できない。第1行と第2行を交換すれば、左上の成分は2になり、これをピボットとして消去を進められる。
行交換を右辺も含めて描くと、次のようになる（\(b=(1,5)^T\) の例）。
青い2が左上のピボット位置へ移る。右辺の1と5も同じ行と一緒に移る。
ピボットが0でなくても、絶対値が非常に小さいと、消去係数 \(m_{ik}\) が大きくなり、途中の値や丸め誤差が増幅されることがある。そこで第 \(k\) 段階では、まだ消去していない第 \(k\) 行から第 \(n\) 行の中から
を選び、第 \(k\) 行と第 \(p\) 行を交換してから、\(a_{kk}^{(k-1)}\) をピボットとして使う。この方法を部分ピボット選択（partial pivoting） という。「部分」は、未処理の部分行列全体ではなく、消去中の第 \(k\) 列から候補を選ぶことを表す。候補を列の中から選ぶことと、列そのものを入れ替えることは別である。 実際に入れ替えるのは行である。
行交換をまとめて置換行列 \(P\) で表すと、行交換後の行列 \(PA\) にGauss消去を行うため、分解は
となる。\(P\) を左からかける \(PA\) は行交換を表す。列交換を表すのは \(AP\) のような右からの乗算であり、ここでは行っていない。連立方程式の右辺にも同じ行交換を適用すれば、
であり、方程式の解は変わらない。まず \(Ly=Pb\) を前進代入で解き、次に \(Ux=y\) を後退代入で解けばよい。
Juliaの lu(A) は通常、部分ピボット選択を用いる。返り値 F の F.p は行の並べ替えを表し、
A[F.p, :] ≈ F.L*F.U
が成り立つ。次の例では、最初のピボットが0の行列で実際に行が交換されることを確認する。
A_zero = [0.0 1.0; 2.0 3.0]
b_zero = [1.0, 5.0]
F_zero = lu(A_zero)
println("選ばれた行の順番 = ", F_zero.p)
println("行交換後の行列 = \n", A_zero[F_zero.p, :])
println("L*U = \n", F_zero.L*F_zero.U)
x_zero = F_zero \ b_zero
println("解 = ", x_zero)
println("残差 = ", norm(b_zero-A_zero*x_zero))
選ばれた行の順番 = [2, 1]
行交換後の行列 =
[2.0 3.0; 0.0 1.0]
L*U =
[2.0 3.0; 0.0 1.0]
解 = [1.0, 1.0]
残差 = 0.0
次に、3元の連立方程式の係数行列を lu で分解する。ここではJuliaの部分ピボット付き分解と \(PA=LU\) の確認に焦点を絞る。したがって、A ではなく行を F.p の順に並べた A[F.p, :] と F.L*F.U を比較する。同じ行列を行交換なしで消去し、基本行列から \(A=LU\) を導く計算は後の補足にまとめる。行交換の有無によって \(L,U\) は異なるが、それぞれ \(PA=LU\) と \(A=LU\) を満たす。
A = [2.0 1.0 1.0; 4.0 -6.0 0.0; -2.0 7.0 2.0]
F = lu(A)
println("L = \n", F.L)
println("U = \n", F.U)
println("行の並べ替え = ", F.p)
println("分解の確認 = ", norm(A[F.p,:] - F.L*F.U))
b = [5.0, -2.0, 9.0]
x = F \ b
println("解 = ", x)
println("残差 = ", norm(b-A*x))
L =
[1.0 0.0 0.0; 0.5 1.0 0.0; -0.5 1.0 1.0]
U =
[4.0 -6.0 0.0; 0.0 4.0 1.0; 0.0 0.0 1.0]
行の並べ替え = [2, 1, 3]
分解の確認 = 0.0
解 = [1.0, 1.0, 2.0]
残差 = 0.0
5.5.1. 0でなくても、小さなピボットには注意する#
\(A=\begin{bmatrix}\varepsilon&1\\1&1\end{bmatrix}\)、\(b=(1,2)^T\)、\(\varepsilon=10^{-20}\) とする。連立方程式は
である。行交換をしないと、最初のピボットは \(\varepsilon\) であり、第2行の \(x_1\) を消すための消去係数は
となる。行操作 \(R_2\leftarrow R_2-mR_1\) を行うと、第2式は
となる。\(m\varepsilon=1\) なので \(x_1\) の項が消え、後退代入により
を得る。その後、第1式 \(\varepsilon x_1+x_2=1\) に戻って
と求める。下の x2_no_pivot と x_no_pivot は、この2式をそのままコードにしたものである。
ただし、倍精度浮動小数点数では \(2-10^{20}\) と \(1-10^{20}\) の小さな差を保持できず、どちらも実質的に \(-10^{20}\) となる。そのため計算上は \(x_2=1\)、続いて \(x_1=(1-1)/\varepsilon=0\) となり、x_no_pivot = [0.0, 1.0] が得られる。一方、厳密解は
であり、\(\varepsilon=10^{-20}\) ではほぼ \((1,1)^T\) である。
この行列の条件数は約2.6である。第4章の言葉では、ここでの難しさは問題の悪条件性ではなく、計算手順にある。部分ピボット選択はこの例を改善し、通常の問題で丸めの増幅を抑える。ただし、あらゆる行列で小さな前進誤差を保証するわけではない。
eps_pivot = 1.0e-20
A_small = [eps_pivot 1.0; 1.0 1.0]
b_small = [1.0, 2.0]
multiplier = 1/eps_pivot # 行交換なしの消去係数 m
# R2 ← R2 - m*R1 により (1-m)x₂ = 2-m
x2_no_pivot = (2-multiplier)/(1-multiplier)
# 第1式 εx₁+x₂=1 に戻って x₁ を求める
x_no_pivot = [(1-x2_no_pivot)/eps_pivot, x2_no_pivot]
x_pivot = A_small \ b_small
println("条件数 = ", cond(A_small))
println("行交換なし = ", x_no_pivot, ", 残差 = ", norm(b_small-A_small*x_no_pivot))
println("部分ピボット選択 = ", x_pivot, ", 残差 = ", norm(b_small-A_small*x_pivot))
条件数 = 2.6180339887498953
行交換なし = [0.0, 1.0], 残差 = 1.0
部分ピボット選択 = [1.0, 1.0], 残差 = 0.0
5.6. n元の消去をコードにする#
ここまでは2元・3元の例を手で追った。同じ手順を \(n\) 元へ広げると、次の二重の繰り返しになる。
第 \(k\) 列（\(k=1,\ldots,n-1\)）について、ピボット \(u_{kk}\) を確認する。
その下の各行 \(i=k+1,\ldots,n\) で消去係数 \(m_{ik}=u_{ik}/u_{kk}\) を計算し、\(L\) の \((i,k)\) 成分に記録する。
第 \(i\) 行から第 \(k\) 行の \(m_{ik}\) 倍を引く。
本文の3元の行列を、行交換せずに消去する流れを先に見よう。各段階で、青い枠がピボット、橙色がその列で消す成分を示す。
第1列を消してから第2列へ進む。最後の行操作では第2行の第1成分が0なので、すでに作った第1列の0を保てる。記録する倍率は \(m_{21}=2\)、\(m_{31}=-1\)、\(m_{32}=-1\) である。
次の lu_nopivot は、理解のために行交換を行わないGauss消去をそのまま書いたものである。実際の計算では、部分ピボット選択を行う lu を使う。前進代入 forward_substitution は、後退代入を上から下へ進める形に書き換えたものである。
実行前に、本文の3元の例で \(L\) の対角線より下に何が並ぶか予想してみよう（後の補足で手計算した \(m_{21},m_{31},m_{32}\) である）。
function lu_nopivot(A)
n = size(A, 1)
size(A, 2) == n || error("正方行列を入力してください")
U = Matrix{Float64}(A)
L = Matrix{Float64}(I, n, n)
for k in 1:n-1 # 第k列を消去する
iszero(U[k,k]) && error("ピボットが0です（第$(k)段階）。行交換が必要です")
for i in k+1:n # ピボットより下の各行
m = U[i,k]/U[k,k] # 消去係数
L[i,k] = m # 倍率をLに記録する
U[i,k:n] .-= m .* U[k,k:n] # R_i ← R_i - m R_k
end
end
return L, U
end
function forward_substitution(L, b)
n = length(b)
y = zeros(Float64, n)
for i in 1:n
known = sum((L[i,j]*y[j] for j in 1:i-1); init=0.0)
y[i] = (b[i] - known)/L[i,i]
end
return y
end
A3 = [2.0 1.0 1.0; 4.0 -6.0 0.0; -2.0 7.0 2.0]
b3 = [5.0, -2.0, 9.0]
L3, U3 = lu_nopivot(A3)
println("L = ", L3)
println("U = ", U3)
println("||A-LU|| = ", norm(A3-L3*U3))
y3 = forward_substitution(L3, b3) # Ly = b
x3 = backward_substitution(U3, y3) # Ux = y
println("y = ", y3, ", x = ", x3, ", 残差 = ", norm(b3-A3*x3))
try
lu_nopivot([0.0 1.0; 2.0 3.0])
catch err
println("行交換なしの失敗: ", err.msg)
end
L = [1.0 0.0 0.0; 2.0 1.0 0.0; -1.0 -1.0 1.0]
U = [2.0 1.0 1.0; 0.0 -8.0 -2.0; 0.0 0.0 1.0]
||A-LU|| = 0.0
y = [5.0, -12.0, 2.0], x = [1.0, 1.0, 2.0], 残差 = 0.0
行交換なしの失敗: ピボットが0です（第1段階）。行交換が必要です
3元の例では、\(L\) の対角線より下に消去係数 \(2,-1,-1\) が並び、解 \((1,1,2)^T\) が得られた。これはJuliaの lu の結果と \(L,U\) は異なるが（あちらは行交換をしている）、解は同じである。最初のピボットが0の行列では、この実装は止まる。lu が行交換を行う理由がここにある。
5.6.1. 演算回数を数える#
第4章と同じように、何回の四則演算が必要かを数えよう。第 \(k\) 段階では、ピボットより下の \(n-k\) 行それぞれについて、消去係数の除算が1回、第 \(k\) 列以降の約 \(n-k\) 個の成分に乗算と減算が1回ずつ必要である。したがって、おおよそ
となる。\(\sum_{j=1}^{n-1}j^2=(n-1)n(2n-1)/6\approx n^3/3\) を使った。一方、前進代入・後退代入は各行で約 \(2i\) 回の演算を行うので、合わせて約 \(2n^2\) 回である。
作業 |
演算回数の目安 |
\(n\) を倍にすると |
|---|---|---|
LU分解 |
\(\frac{2}{3}n^3\) |
約8倍 |
前進代入＋後退代入（右辺1本） |
\(2n^2\) |
約4倍 |
次のコードは、上の和を実際に計算して \(\frac23n^3\) と比べ、さらにJuliaの lu と代入の実行時間を測る。
for n in (100, 200, 400)
elimination = sum(2*(n-k)^2 + (n-k) for k in 1:n-1)
println("n=$n: 消去の演算回数=", elimination, ", 2n^3/3=", round(Int, 2n^3/3))
end
using Random
Random.seed!(1)
function best_time(f; samples=3)
f() # 初回のコンパイルを先に済ませる
return minimum([@elapsed f() for _ in 1:samples])
end
for n in (500, 1000, 2000)
A_rand = randn(n, n)
b_rand = randn(n)
F_rand = lu(A_rand)
t_factor = best_time(() -> lu(A_rand))
t_solve = best_time(() -> F_rand \ b_rand)
println("n=$n: 分解 ", round(t_factor*1e3, sigdigits=3), " ms, 代入 ",
round(t_solve*1e3, sigdigits=3), " ms")
end
n=100: 消去の演算回数=661650, 2n^3/3=666667
n=200: 消去の演算回数=5313300, 2n^3/3=5333333
n=400: 消去の演算回数=42586600, 2n^3/3=42666667
n=500: 分解 1.39 ms, 代入 0.0412 ms
n=1000: 分解 6.15 ms, 代入 0.149 ms
n=2000: 分解 30.9 ms, 代入 0.732 ms
数えた演算回数は \(\frac23n^3\) とよく一致し、\(n\) を倍にすると約8倍になる。一方、実行時間は計算機によって異なり、8倍ちょうどにはならないことが多い。lu の内部では、行列計算ライブラリ（BLAS/LAPACK）が複数のコアやキャッシュを活用するためである（第16章）。演算回数は規模の見通しを与え、実測はその計算機での実際を与える。両者を区別して読もう。
それでも「分解は代入の数十倍以上かかる」という関係は変わらない。同じ行列で右辺だけを変えて何度も解くなら、分解を1回だけ行って再利用する利益は大きい。
逆行列を計算しないのはなぜか
\(x=A^{-1}b\) を inv(A)*b と計算すると、逆行列を作るのに約 \(2n^3\) 回の演算が必要で、LU分解の約3倍である。丸め誤差の面でも、一般に A \ b の方が有利である。「逆行列を掛ける」は式の上の表現であり、計算では「分解して代入する」と読み替える。ただし、次の応用例のように、逆行列の各成分そのものに意味があり、それを知りたい場合は別である。
5.7. 応用例：産業連関分析#
経済の各産業は、自分の生産のために他の産業の製品を使う。例えば製造業が1単位生産するには、農林水産業から0.05単位、製造業自身から0.40単位、サービス業から0.20単位を投入する、というように考える。この比率を並べた行列 \(A\) を投入係数行列という。第 \(j\) 列は「産業 \(j\) が1単位生産するのに必要な各産業からの投入」である。
各産業の生産額を \(x\)、消費者や輸出などの最終需要を \(d\) とすると、生産額は「他産業への中間投入」と「最終需要」の和に等しいので、
となる。これは連立一次方程式であり、Leontief（レオンチェフ）モデルと呼ばれる。実際の産業連関表は数百部門に及び、政府や自治体が経済波及効果を見積もる際に使われている。
ここでは3部門の架空の例を使う。最終需要をいくつも変えて比較する（シナリオ分析）ことが多いので、\(I-A\) を一度だけ分解し、各シナリオで代入だけを行う。
Ain = [0.10 0.05 0.02; # 行: 投入する側（農林水産, 製造, サービス）
0.20 0.40 0.15; # 列: 生産する側
0.10 0.20 0.25]
d_base = [50.0, 300.0, 400.0] # 最終需要
M_leontief = I - Ain
F_leontief = lu(M_leontief) # 一度だけ分解する
x_base = F_leontief \ d_base
println("投入係数の列和 = ", sum(Ain, dims=1))
println("必要な生産額 = ", round.(x_base, digits=1))
println("残差 = ", norm(M_leontief*x_base - d_base))
for (label, Δd) in (("製造業への最終需要 +30", [0.0, 30.0, 0.0]),
("サービス業への最終需要 +30", [0.0, 0.0, 30.0]))
Δx = F_leontief \ Δd # 分解を再利用して代入だけを行う
println(label, ": 生産額の増加 = ", round.(Δx, digits=2),
", 合計 = ", round(sum(Δx), digits=2))
end
投入係数の列和 = [0.4 0.65 0.42]
必要な生産額 = [112.2, 722.6, 741.0]
残差 = 1.6077746776921858e-13
製造業への最終需要 +30: 生産額の増加 = [3.39, 54.9, 15.09], 合計 = 73.38
サービス業への最終需要 +30: 生産額の増加 = [1.59, 11.34, 43.24], 合計 = 56.17
製造業への最終需要を30増やすと、製造業自身の生産は30より多い約54.9増え、農林水産業・サービス業の生産も増える。製造業が原材料を他産業から買い、その産業がさらに別の産業から買う、という連鎖が方程式に含まれているためである。生産額の増加の合計（約73.4）を最初の需要増30で割った値は、経済波及効果の大きさの目安になる。
\((I-A)^{-1}\) はレオンチェフ逆行列と呼ばれ、その第 \(j\) 列は「産業 \(j\) の最終需要が1増えたときの各産業の生産増」を表す。このように逆行列の成分そのものに意味がある場合には、単位行列の各列を右辺として F_leontief \ Matrix(I, 3, 3) のように計算する。これも \(n\) 本の右辺に対する代入であり、分解は1回で済む。
モデルとしての注意もある。投入係数が一定、価格が変わらない、生産能力に上限がない、という仮定を置いている（第1章のモデル誤差）。また、投入係数の列和が1未満でなければ \(x\ge0\) となる解が得られるとは限らない。数値計算で解が求まったことと、その解が経済的に意味を持つことは、別に確かめる必要がある。
5.8. \(\clubsuit\) 基本行列で消去を表す#
ここでは、本文のJulia例と同じ行列に対し、行交換を行わない消去を基本行列の積として確認する。このため、Juliaが部分ピボット選択で返した \(L,U\) と、ここで得る \(L,U\) は異なる。一般の分解を導くための補足であり、基本の解法と行交換は前の節だけで追える。
行基本変形は、行列を左からかける操作として表せる。第 \(i\) 行から第 \(k\) 行の \(m\) 倍を引く操作
に対応する基本行列を
と書く。ここで \(e_i\) は第 \(i\) 成分だけが1の単位ベクトルである。\(E_{ik}(m)A\) を計算すると、\(A\) の第 \(i\) 行だけが \(R_i-mR_k\) に置き換わる。
具体例として
を上三角行列へ変形する。第1列の消去係数は
なので、対応する基本行列は
これらを順に左からかけると
次に \(m_{32}=8/(-8)=-1\) として
を左からかければ
となる。ここで重要なのは、\(U\) が突然現れたのではなく、消去を表す基本行列を順番に左からかけた結果だという点である。
上式の消去操作を逆向きに戻すと
したがって
※ここでは、基本行列の逆操作が、上の基本行列の係数を反転させた操作であることを使った。基本行列の逆行列では、引いた操作を足して戻すため符号が反転する。
結果として、Gauss消去で使った係数 \(m_{21},m_{31},m_{32}\) が \(L\) の対角線より下に並ぶ。
この計算で消去係数の分母に使った \(a_{11}=2\) と \(a_{22}^{(1)}=-8\) が、前節で定義したピボットである。どちらも0でないため、この補足では行を入れ替えずに消去できた。
A_basic = [2.0 1.0 1.0; 4.0 -6.0 0.0; -2.0 7.0 2.0]
I3 = Matrix{Float64}(I,3,3)
m21 = A_basic[2,1]/A_basic[1,1]
m31 = A_basic[3,1]/A_basic[1,1]
E21 = copy(I3); E21[2,1] = -m21
E31 = copy(I3); E31[3,1] = -m31
A_step1 = E31*E21*A_basic
m32 = A_step1[3,2]/A_step1[2,2]
E32 = copy(I3); E32[3,2] = -m32
U_basic = E32*A_step1
# 各消去操作を逆向きに戻す基本行列
E21_inv = copy(I3); E21_inv[2,1] = m21
E31_inv = copy(I3); E31_inv[3,1] = m31
E32_inv = copy(I3); E32_inv[3,2] = m32
L_basic = E21_inv*E31_inv*E32_inv
println("消去係数 = ", (m21,m31,m32))
println("U = \n", U_basic)
println("L = \n", L_basic)
println("||A-LU|| = ", norm(A_basic-L_basic*U_basic))
消去係数 = (2.0, -1.0, -1.0)
U =
[2.0 1.0 1.0; 0.0 -8.0 -2.0; 0.0 0.0 1.0]
L =
[1.0 0.0 0.0; 2.0 1.0 0.0; -1.0 -1.0 1.0]
||A-LU|| = 0.0
5.9. まとめ#
本章のまとめ
Gauss消去法は、係数と右辺へ同じ行操作を施して三角形の方程式を作る。LU分解はその操作を再利用するための記録で、行交換がある場合は \(PA=LU\) と書く。行交換は丸めの増幅を抑える工夫であり、元の問題の悪条件性を解消するものではない。
密な \(n\times n\) 行列では、分解が \(O(n^3)\)、1本の右辺に対する前進・後退代入が \(O(n^2)\) である。右辺が3本なら、分解を3回行う代わりに、1回の分解と3回の代入で済む。
1回の求解：
x = A \ b同じ行列で繰り返す：
F = lu(A)の後にF \ b1,F \ b2検証：元の式の相対残差
norm(b-A*x)/norm(b)と、必要に応じて条件数（ここではb != 0）。
逆行列は、その成分自体に意味がある場合（レオンチェフ逆行列など）を除き、明示的に作らない。
未知数が数十万を超える疎行列では、分解のメモリが足りなくなることがある。そのような場合に使う反復法は、第10章「連立一次方程式の反復法」で扱う。産業連関の「波及を1段ずつ足す」計算が、その入口になる。
次の第6章では、方程式が未知数より多く、すべてを同時には満たせない場合を考える。
5.10. 演習#
本文の式・表・出力を使って説明してよい。
演習
消去と右辺：本文の2元連立方程式で、第2式から第1式の2倍を引き、消去後の右辺が1になることを確かめよ。その後、下の式から順に解け。
行交換の意味：\(A=\begin{bmatrix}0&1\\2&3\end{bmatrix}\)、\(b=(1,5)^T\) で、最初に行交換する理由を説明せよ。右辺にも同じ交換をした式を書き、\(x=(1,1)^T\) を代入して確かめよ。
分解を使い回す：同じ係数行列で右辺だけが3回変わる場合、LU分解のどの部分を再利用できるか。問題サイズを倍にすると、分解と1回の代入の仕事量はそれぞれ約何倍になるか。
実装を確かめる：
lu_nopivotを本文の2元の例 \(A=\begin{bmatrix}2&1\\4&3\end{bmatrix}\) に使い、手計算の \(L,U\) と一致するか確かめよ。次に、小さなピボットの例（\(\varepsilon=10^{-20}\)）に使うと何が起こるかを実行前に予想して書き、実行結果と比べよ。予想と違った点があれば、第4章の「条件数」と「安定性」のどちらの問題かを説明せよ。探究：波及効果を比べる：産業連関の例で、どの産業の最終需要を1増やすと生産額の合計が最も増えるか。まずレオンチェフ逆行列の列和から予想し、計算で確かめよ。次に、製造業の自部門投入係数（本文の0.40）を0.60、0.80と大きくしたとき、生産額と
cond(I - Ain)がどう変わるかを調べ、列和が1に近づくと何が起こるかを考察せよ。結論だけでなく、どの値を変えて何を比べたかを記録すること。
5.11. 参考文献と関連する節#
Driscoll and Braun, Fundamentals of Numerical Computation, LU factorization.