前進代入

この章の目標

下三角行列を係数行列とする連立一次方程式を前進代入で解く手順を理解する。計算量 $O(n^2)$、単位下三角行列での簡略化、LU分解・コレスキー分解における役割を把握する。

前提知識

  • 行列とベクトルの基本演算
  • 下三角行列・上三角行列の定義
目次

1. 定義と原理

前進代入（forward substitution）は、下三角行列 $L$ を係数行列とする連立一次方程式

$$L\mathbf{y} = \mathbf{b}$$

を解くアルゴリズムである。下三角行列では $l_{ij} = 0 \; (i < j)$ であるから、第1行は1変数のみの方程式 $l_{11} y_1 = b_1$ となり、$y_1$ が直ちに求まる。得られた $y_1$ を第2行に代入すれば $y_2$ が求まり、これを最後の行まで繰り返すことで全ての変数を決定できる。

展開して書くと

$$\begin{pmatrix} l_{11} & 0 & \cdots & 0 \\ l_{21} & l_{22} & \cdots & 0 \\ \vdots & & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix}$$

である。$l_{ii} \neq 0$（$L$ が非特異）であれば解が一意に存在する。

前進代入：第1行の y₁ から始め、上から下へ順に未知数を求める 2 0 0 1 3 0 −1 2 4 y₁ y₂ y₃ = 4 7 10 1 2 3 上から下へ 2·y₁ = 4 → y₁ = 2 1·y₁ + 3·y₂ = 7 −y₁ + 2·y₂ + 4·y₃ = 10 対角の上はすべて 0 だから、各行は新しい未知数を1つだけ含む。 ① 第1行 l₁₁y₁ = b₁ から y₁ が即座に求まる（緑＝確定）。 ② 求めた y₁ を次の行に代入して y₂ を求める。 ③ 一般に yᵢ = (bᵢ − 既知項) / lᵢᵢ。上から下へ繰り返す。 後退代入（上三角を下から上へ解く）のちょうど鏡像である。
図1. 前進代入。下三角行列 $L\mathbf{y}=\mathbf{b}$ では対角の上が $0$ なので、第1行の $y_1$ から始め、求めた値を下の行へ代入しながら上から下へ順に解く。各 $y_i = (b_i - \sum_{j<i} l_{ij}y_j)/l_{ii}$ で、後退代入の鏡像になっている。

2. 公式

$$y_1 = \dfrac{b_1}{l_{11}}$$
$$y_i = \dfrac{1}{l_{ii}} \left(b_i - \displaystyle\sum_{j=1}^{i-1} l_{ij} y_j\right), \quad i = 2, 3, \ldots, n$$

各 $y_i$ は、右辺 $b_i$ から既に求まっている $y_1, \ldots, y_{i-1}$ の寄与を差し引いた後、対角要素 $l_{ii}$ で割ることで得られる。後退代入とは逆方向に進む点が対称的である。

3. アルゴリズム

function forwardSubstitution(L, b, n):
    y[1] = b[1] / L[1][1]
    for i = 2 to n:
        sum = 0
        for j = 1 to i-1:
            sum = sum + L[i][j] * y[j]
        y[i] = (b[i] - sum) / L[i][i]
    return y

4. 単位下三角行列の場合

LU 分解で得られる $L$ はしばしば単位下三角行列（unit lower triangular matrix）、すなわち $l_{ii} = 1$ である。この場合、除算が不要になり

$$y_i = b_i - \displaystyle\sum_{j=1}^{i-1} l_{ij} y_j$$

で直接 $y_i$ が求まる。除算 $n$ 回分の計算が節約されるとともに、除算に伴う丸め誤差も生じず、乗算・加減算による丸め誤差だけが残る。

5. 計算量

演算一般の $L$単位下三角 $L$
除算$n$$0$
乗算$\dfrac{n(n-1)}{2}$$\dfrac{n(n-1)}{2}$
加減算$\dfrac{n(n-1)}{2}$$\dfrac{n(n-1)}{2}$
合計$n^2$ ($O(n^2)$)$n^2 - n$ ($O(n^2)$)

6. LU分解・コレスキー分解での役割

LU 分解

$PA = LU$ の LU 分解が計算済みの場合、$A\mathbf{x} = \mathbf{b}$ は以下の2ステップで解ける。

  1. $L\mathbf{y} = P\mathbf{b}$ を前進代入で解く。
  2. $U\mathbf{x} = \mathbf{y}$ を後退代入で解く。

コレスキー分解

正定値対称行列 $A = LL^T$ のコレスキー分解の場合も同様に

  1. $L\mathbf{y} = \mathbf{b}$ を前進代入で解く。
  2. $L^T \mathbf{x} = \mathbf{y}$ を後退代入で解く。

コレスキー分解は $L$ が一般の下三角行列（$l_{ii} \neq 1$）であるため、前進代入で除算が必要になる。

7. 数値安定性

前進代入は後退代入と同様に後退安定であり、丸め誤差の増大は穏やかである。LU 分解で得られる単位下三角行列 $L$ の場合、部分ピボット選択により $|l_{ij}| \leq 1$ が保証されるため、特に安定した計算が期待できる。

8. 計算例

例1: 3x3 下三角系

$\begin{pmatrix} 2 & 0 & 0 \\ 1 & 3 & 0 \\ -1 & 2 & 4 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} 4 \\ 7 \\ 10 \end{pmatrix}$

ステップ計算結果
$y_1$$4 / 2$$y_1 = 2$
$y_2$$(7 - 1 \times 2) / 3$$y_2 = 5/3$
$y_3$$(10 - (-1) \times 2 - 2 \times 5/3) / 4$$y_3 = 23/12$

解は $\mathbf{y} = (2,\; 5/3,\; 23/12)^T$ である。

例2: 単位下三角行列の場合

$\begin{pmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ -0.25 & 0.75 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} 6 \\ 1 \\ 5 \end{pmatrix}$

ステップ計算結果
$y_1$$6$$y_1 = 6$
$y_2$$1 - 0.5 \times 6$$y_2 = -2$
$y_3$$5 - (-0.25) \times 6 - 0.75 \times (-2)$$y_3 = 8$

除算が不要であることに注意されたい。

9. よくある質問

Q1. 前進代入とは？

下三角行列 $L\mathbf{y} = \mathbf{b}$ を最初の変数 $y_1$ から順方向に求めていくアルゴリズムである。LU分解やコレスキー分解で使われる基本アルゴリズムである。

Q2. 計算量は？

$n \times n$ の下三角行列に対して $O(n^2)$ である。後退代入と同じオーダーである。

Q3. 単位下三角行列の場合は？

$l_{ii} = 1$ であるため除算が不要になり、$y_i = b_i - \displaystyle\sum_{j=1}^{i-1} l_{ij} y_j$ で直接求まる。

10. 参考資料