Pythonで学ぶ数値計算のアルゴリズムと実装 第7回：行列分解の活用：LU分解とコレスキー分解

こんにちは、JS2IIUです。
前回の第6回では、連立一次方程式を解くための最も基本的な手法である「ガウス消去法」を学びました。ガウス消去法は非常に強力ですが、一つの大きな弱点があります。それは、係数行列 \(A\) は同じままで、右辺のベクトル \(b\) だけが次々と変わるような状況において、毎回 \(O(n^3)\) の重い計算を最初からやり直さなければならないという点です。

機械学習の最適化計算や、リアルタイムの物理シミュレーションにおいて、これは致命的な非効率を招きます。そこで登場するのが「行列分解（Matrix Decomposition）」という考え方です。行列 \(A\) をあらかじめ「扱いやすい二つの行列」に分解して保存しておけば、右辺が変わっても非常に高速に解を求めることができます。

今回は、行列分解の代表格である「LU分解」と、AI分野で頻出する対称正定値行列に特化した「コレスキー分解」について、そのアルゴリズムとPython/PyTorchでの実装を詳しく解説します。今回もよろしくお願いします。

1. なぜ行列を「分解」するのか？

連立一次方程式 \(Ax = b\) を解くプロセスを、以下の二つのフェーズに分けるのが行列分解の戦略です。

1. 分解フェーズ (\(O(n^3)\)): 行列 \(A\) を \(L\)（下三角行列）と \(U\)（上三角行列）の積に分解する。この計算は \(A\) が決まった時点で一度だけ行えば良い。
2. 代入フェーズ (\(O(n^2)\)): 分解された \(L\) と \(U\) を使い、与えられた \(b\) に対して解 \(x\) を求める。これは第6回で学んだ「後退代入」の親戚のような計算で、非常に高速。

ここで、下三角行列 \(L\) と上三角行列 \(U\) は、それぞれ以下のような形をしています。

$$L = \begin{bmatrix}
l_{11} & 0 & 0 & \cdots & 0 \\
l_{21} & l_{22} & 0 & \cdots & 0 \\
l_{31} & l_{32} & l_{33} & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
l_{n1} & l_{n2} & l_{n3} & \cdots & l_{nn}
\end{bmatrix}, \quad
U = \begin{bmatrix}
u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\
0 & u_{22} & u_{23} & \cdots & u_{2n} \\
0 & 0 & u_{33} & \cdots & u_{3n} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & u_{nn}
\end{bmatrix}$$

\(L\) は対角線より上側の要素がすべて 0、\(U\) は対角線より下側の要素がすべて 0 という構造を持っています。この特殊な構造により、連立方程式の解法が劇的に高速化されます。

例えば、1000次元の行列に対して、右辺のデータが100パターンある場合、ガウス消去法を100回繰り返すと莫大な時間がかかりますが、LU分解なら「重い計算1回 + 軽い計算100回」で済みます。これが、行列分解が数値計算の標準となっている最大の理由です。

2. LU分解のアルゴリズム：ガウス消去法の再構築

LU分解とは、行列 \(A\) を以下の二つの行列の積に変形することです。
$$A = LU$$
ここで、\(L\) (Lower) は対角成分がすべて 1 の下三角行列、\(U\) (Upper) は上三角行列です。

実は、LU分解のプロセスはガウス消去法の「前進消去」そのものです。前進消去で得られる結果が \(U\) であり、その過程で各行を何倍して引いたかという「倍率」を並べたものが \(L\) になります。

LU分解による解法の手順

\(Ax = b\) を \((LU)x = b\) と書き換えます。

1. まず \(Ly = b\) を解いて \(y\) を求める（前向き代入：Forward Substitution）。
2. 次に \(Ux = y\) を解いて \(x\) を求める（後ろ向き代入：Backward Substitution）。

具体的には、以下の手順で計算が進みます。

ステップ1: 前向き代入 \(Ly = b\)

下三角行列 \(L\) の構造を利用して、上から順に \(y\) の要素を求めていきます。

$$\begin{align}
y_1 &= \frac{b_1}{l_{11}} \\
y_2 &= \frac{b_2 – l_{21}y_1}{l_{22}} \\
y_3 &= \frac{b_3 – l_{31}y_1 – l_{32}y_2}{l_{33}} \\
&\vdots \\
y_i &= \frac{b_i – \sum_{j=1}^{i-1} l_{ij}y_j}{l_{ii}} \quad (i = 1, 2, \dots, n)
\end{align}$$

ステップ2: 後ろ向き代入 \(Ux = y\)

上三角行列 \(U\) の構造を利用して、下から順に \(x\) の要素を求めていきます。

$$\begin{align}
x_n &= \frac{y_n}{u_{nn}} \\
x_{n-1} &= \frac{y_{n-1} – u_{n-1,n}x_n}{u_{n-1,n-1}} \\
x_{n-2} &= \frac{y_{n-2} – u_{n-2,n-1}x_{n-1} – u_{n-2,n}x_n}{u_{n-2,n-2}} \\
&\vdots \\
x_i &= \frac{y_i – \sum_{j=i+1}^{n} u_{ij}x_j}{u_{ii}} \quad (i = n, n-1, \dots, 1)
\end{align}$$

どちらも三角行列を相手にするため、各ステップで既に求めた値だけを使って次の値を計算でき、全体で \(O(n^2)\) の計算量で完了します。

3. ピボット選択付きLU分解 (\(PA = LU\))

第6回で学んだ通り、数値計算の安定性を保つためには「ピボット選択（行の入れ替え）」が不可欠です。行を入れ替えた場合のLU分解は、数学的には以下のように表現されます。
$$PA = LU$$
ここで \(P\) は行の入れ替えを記録した「置換行列（Permutation Matrix）」です。実際のライブラリ（PyTorchの torch.linalg.lu など）では、この形式で分解が行われます。

PyTorchによるLU分解のスクラッチ実装

それでは、理解を深めるために、ピボット選択を考慮しない最も基本的なLU分解（Doolittle法）を実装してみましょう。

import torch

def lu_decomposition(A):
    """
    行列 A を LU 分解する (A = LU)
    注: 簡略化のためピボット選択は省略

    Args:
        A (torch.Tensor): n x n の正方行列

    Returns:
        L (torch.Tensor): 下三角行列
        U (torch.Tensor): 上三角行列
    """
    n = A.shape[0]
    L = torch.eye(n) # 対角成分が1の単位行列で初期化
    U = A.clone().float()

    for k in range(n):
        for i in range(k + 1, n):
            # ガウス消去法の倍率を L に格納
            factor = U[i, k] / U[k, k]
            L[i, k] = factor

            # i行目の更新 (前進消去)
            U[i, k:] -= factor * U[k, k:]

    return L, U

def solve_lu(L, U, b):
    """
    分解された L, U を用いて Ax = b を解く
    """
    n = L.shape[0]

    # 1. 前向き代入: Ly = b
    y = torch.zeros_like(b)
    for i in range(n):
        y[i] = b[i] - torch.dot(L[i, :i], y[:i])

    # 2. 後ろ向き代入: Ux = y
    x = torch.zeros_like(b)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - torch.dot(U[i, i+1:], x[i+1:])) / U[i, i]

    return x

# テスト実行
A_sample = torch.tensor([[2, 3, 1], [4, 1, 3], [2, -1, -1]], dtype=torch.float32)
b_sample = torch.tensor([1, 1, -3], dtype=torch.float32)

L, U = lu_decomposition(A_sample)
x_sol = solve_lu(L, U, b_sample)

print(f"A_sample:\n{A_sample}")
print(f"L:\n{L}")
print(f"U:\n{U}")
print(f"Solution x: {x_sol}")
# 確認: A @ x_sol が b に近いか
print(f"Check Ax: {torch.matmul(A_sample, x_sol)}")

import torch

def lu_decomposition(A):
    """
    行列 A を LU 分解する (A = LU)
    注: 簡略化のためピボット選択は省略

    Args:
        A (torch.Tensor): n x n の正方行列

    Returns:
        L (torch.Tensor): 下三角行列
        U (torch.Tensor): 上三角行列
    """
    n = A.shape[0]
    L = torch.eye(n) # 対角成分が1の単位行列で初期化
    U = A.clone().float()

    for k in range(n):
        for i in range(k + 1, n):
            # ガウス消去法の倍率を L に格納
            factor = U[i, k] / U[k, k]
            L[i, k] = factor

            # i行目の更新 (前進消去)
            U[i, k:] -= factor * U[k, k:]

    return L, U

def solve_lu(L, U, b):
    """
    分解された L, U を用いて Ax = b を解く
    """
    n = L.shape[0]

    # 1. 前向き代入: Ly = b
    y = torch.zeros_like(b)
    for i in range(n):
        y[i] = b[i] - torch.dot(L[i, :i], y[:i])

    # 2. 後ろ向き代入: Ux = y
    x = torch.zeros_like(b)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - torch.dot(U[i, i+1:], x[i+1:])) / U[i, i]

    return x

# テスト実行
A_sample = torch.tensor([[2, 3, 1], [4, 1, 3], [2, -1, -1]], dtype=torch.float32)
b_sample = torch.tensor([1, 1, -3], dtype=torch.float32)

L, U = lu_decomposition(A_sample)
x_sol = solve_lu(L, U, b_sample)

print(f"A_sample:\n{A_sample}")
print(f"L:\n{L}")
print(f"U:\n{U}")
print(f"Solution x: {x_sol}")
# 確認: A @ x_sol が b に近いか
print(f"Check Ax: {torch.matmul(A_sample, x_sol)}")

A_sample:
tensor([[ 2.,  3.,  1.],
        [ 4.,  1.,  3.],
        [ 2., -1., -1.]])
L:
tensor([[1.0000, 0.0000, 0.0000],
        [2.0000, 1.0000, 0.0000],
        [1.0000, 0.8000, 1.0000]])
U:
tensor([[ 2.0000,  3.0000,  1.0000],
        [ 0.0000, -5.0000,  1.0000],
        [ 0.0000,  0.0000, -2.8000]])
Solution x: tensor([-0.7143,  0.4286,  1.1429])
Check Ax: tensor([ 1.0000,  1.0000, -3.0000])

A_sample:
tensor([[ 2.,  3.,  1.],
        [ 4.,  1.,  3.],
        [ 2., -1., -1.]])
L:
tensor([[1.0000, 0.0000, 0.0000],
        [2.0000, 1.0000, 0.0000],
        [1.0000, 0.8000, 1.0000]])
U:
tensor([[ 2.0000,  3.0000,  1.0000],
        [ 0.0000, -5.0000,  1.0000],
        [ 0.0000,  0.0000, -2.8000]])
Solution x: tensor([-0.7143,  0.4286,  1.1429])
Check Ax: tensor([ 1.0000,  1.0000, -3.0000])

コードの解説

• lu_decomposition: ガウス消去法の過程で発生する「倍率（factor）」を \(L\) に記録し、消去後の行列を \(U\) としています。
• solve_lu: 二つの代入ステップを連続して実行します。この関数は \(O(n^2)\) なので、一度 \(L, U\) を求めてしまえば、異なる \(b\) に対して非常に高速に動作します。

4. コレスキー分解：対称正定値行列のための高速戦略

機械学習、特にガウス過程や最小二乗法、あるいは統計学の共分散行列を扱う際、行列はしばしば「対称正定値（Symmetric Positive Definite）」という性質を持ちます。これは、行列が \(A = A^T\) であり、かつすべての固有値が正であることを意味します。

この性質を持つ行列に対しては、LU分解をさらに特殊化・効率化した「コレスキー分解（Cholesky Decomposition）」が使えます。
コレスキー分解では、\(A\) を次のように一つの下三角行列 \(L\) とその転置 \(L^T\) の積に分解します：

$$A = LL^T$$

\(A\) が \(n \times n\) の対称正定値行列のとき、\(L\) は次のような下三角行列です：

$$L = \begin{bmatrix}
l_{11} & 0 & 0 & \cdots & 0 \\
l_{21} & l_{22} & 0 & \cdots & 0 \\
l_{31} & l_{32} & l_{33} & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
l_{n1} & l_{n2} & l_{n3} & \cdots & l_{nn}
\end{bmatrix}$$

各成分 \(l_{ij}\) は以下のように計算されます：

• 対角成分：
  $$l_{jj} = \sqrt{a_{jj} – \sum_{k=1}^{j-1} l_{jk}^2}$$
• 非対角成分（\(i > j\)）：
  $$l_{ij} = \frac{1}{l_{jj}} \left(a_{ij} – \sum_{k=1}^{j-1} l_{ik}l_{jk}\right)$$

LU分解では異なる二つの行列 \(L\) と \(U\) を求めましたが、コレスキー分解では一つの下三角行列 \(L\) とその転置 \(L^T\) だけで済みます。

コレスキー分解のメリット

1. 計算量: LU分解の約半分の手数（\(\frac{1}{3}n^3\) 回の演算）で完了します。
2. メモリ: \(L\) だけを保存すれば良いため、メモリ効率が高いです。
3. 数値的安定性: ピボット選択を行わなくても、数値的に非常に安定していることが数学的に保証されています。

コレスキー分解のPyTorch実装

def cholesky_decomposition(A):
    """
    対称正定値行列 A をコレスキー分解する (A = LL^T)
    """
    n = A.shape[0]
    L = torch.zeros_like(A)

    for j in range(n):
        # 対角成分の計算
        sum_sq = torch.sum(L[j, :j]**2)
        val = A[j, j] - sum_sq
        if val <= 0:
            raise ValueError("行列が正定値ではありません")
        L[j, j] = torch.sqrt(val)

        # 非対角成分の計算
        for i in range(j + 1, n):
            L[i, j] = (A[i, j] - torch.dot(L[i, :j], L[j, :j])) / L[j, j]

    return L

# 対称正定値行列の例 (A = M^T @ M は必ず対称正定値になる)
M = torch.tensor([[3, 1], [1, 2]], dtype=torch.float32)
A_spd = torch.matmul(M.T, M)

L_chol = cholesky_decomposition(A_spd)

print(f"A_spd:\n{A_spd}")
print(f"Cholesky L:\n{L_chol}")
print(f"Verification LL^T:\n{torch.matmul(L_chol, L_chol.T)}")

def cholesky_decomposition(A):
    """
    対称正定値行列 A をコレスキー分解する (A = LL^T)
    """
    n = A.shape[0]
    L = torch.zeros_like(A)

    for j in range(n):
        # 対角成分の計算
        sum_sq = torch.sum(L[j, :j]**2)
        val = A[j, j] - sum_sq
        if val <= 0:
            raise ValueError("行列が正定値ではありません")
        L[j, j] = torch.sqrt(val)

        # 非対角成分の計算
        for i in range(j + 1, n):
            L[i, j] = (A[i, j] - torch.dot(L[i, :j], L[j, :j])) / L[j, j]

    return L

# 対称正定値行列の例 (A = M^T @ M は必ず対称正定値になる)
M = torch.tensor([[3, 1], [1, 2]], dtype=torch.float32)
A_spd = torch.matmul(M.T, M)

L_chol = cholesky_decomposition(A_spd)

print(f"A_spd:\n{A_spd}")
print(f"Cholesky L:\n{L_chol}")
print(f"Verification LL^T:\n{torch.matmul(L_chol, L_chol.T)}")

A_spd:
tensor([[10.,  5.],
        [ 5.,  5.]])
Cholesky L:
tensor([[3.1623, 0.0000],
        [1.5811, 1.5811]])
Verification LL^T:
tensor([[10.,  5.],
        [ 5.,  5.]])

A_spd:
tensor([[10.,  5.],
        [ 5.,  5.]])
Cholesky L:
tensor([[3.1623, 0.0000],
        [1.5811, 1.5811]])
Verification LL^T:
tensor([[10.,  5.],
        [ 5.,  5.]])

5. 逆行列の構成アルゴリズム

エンジニアがよく使う inv(A) は、内部でどのように計算されているのでしょうか。実は、逆行列を直接求める公式（余因子行列など）を使うことはまずありません。

逆行列 \(A^{-1}\) を求めるということは、連立方程式\(AX = I\)（\(I\) は単位行列）を解くことと同じです。
つまり、単位行列の各列（\(e_1, e_2, \dots, e_n\)）を右辺として、LU分解済みの \(L, U\) を使って \(n\) 回の代入計算を行うことで、逆行列の各列を求めるのが最も効率的です。

実践的なアドバイス：安易に逆行列を求めない

数値計算の定石として、「逆行列そのものが必要な場面は非常に少ない」という言葉があります。
多くの場合は \(A^{-1}b\)（逆行列とベクトルの積）が計算できれば十分であり、その場合は solve(A, b) を使うほうが精度も高く、計算も速いです。逆行列を明示的に求めるのは、どうしてもその行列自体を使い回す必要があるときだけにしましょう。

6. PyTorchの標準機能との比較

ここまで自作アルゴリズムを見てきましたが、実務では PyTorch の高度に最適化された関数を使用します。これらは内部で LAPACK という伝統ある高速ライブラリを呼んでおり、マルチコアCPUやGPUをフルに活用します。

# LU分解 (ピボット付き)
# P は置換行列を圧縮したインデックス形式で返されることが多い
A_large = torch.randn(100, 100)
LU, pivots = torch.linalg.lu_factor(A_large)

# コレスキー分解
# A_spd は対称正定値である必要がある
L_torch = torch.linalg.cholesky(A_spd)

# LU分解 (ピボット付き)
# P は置換行列を圧縮したインデックス形式で返されることが多い
A_large = torch.randn(100, 100)
LU, pivots = torch.linalg.lu_factor(A_large)

# コレスキー分解
# A_spd は対称正定値である必要がある
L_torch = torch.linalg.cholesky(A_spd)

まとめ

第7回では、数値計算の実践において極めて重要な「行列分解」について学びました。

• LU分解はガウス消去法を保存可能な形にしたもので、右辺が複数ある問題で圧倒的な効率を誇る。
• ピボット選択 (\(PA=LU\)) は実用上の必須テクニックであり、数値的安定性を守る。
• コレスキー分解 (\(A=LL^T\)) は対称正定値行列に特化しており、LU分解の2倍速く、安定している。
• 逆行列は LU 分解を応用して求められるが、直接 inv を求めるよりも solve を使うほうが推奨されることが多い。

AIアルゴリズムの論文を読んでいると、「Cholesky decomposition」や「LU factorization」という言葉が頻繁に登場します。それらが出てきたとき、その背後でこうした「三角行列への変換」と「高速な代入計算」が行われていることを思い出してください。

次回は、行列の要素の多くが 0 である場合に特化した「帯行列」や「スパース行列」の高速計算について解説します。大規模データ処理には欠かせない知識です。お楽しみに。

最後まで読んでいただきありがとうございます。

連載記事、各回へのリンク

Pythonで学ぶ数値計算のアルゴリズムと実装 第1回：多項式計算とテイラー展開 | アマチュア無線局JS2IIU
https://js2iiu.com/2025/12/31/python_numerical_01/

Pythonで学ぶ数値計算のアルゴリズムと実装 第2回：多項式補間とチェビシェフ補間の戦略 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_02/

Pythonで学ぶ数値計算のアルゴリズムと実装 第3回：積分の基礎：中点則、台形則からシンプソン則まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_03/

Pythonで学ぶ数値計算のアルゴリズムと実装 第4回：高度な数値積分：二重指数型（DE）公式と広義積分 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_04/

Pythonで学ぶ数値計算のアルゴリズムと実装 第5回：計算の「正しさ」を評価する：ノルムと条件数 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numeric_05/

Pythonで学ぶ数値計算のアルゴリズムと実装 第6回：連立一次方程式の直截解法：ガウス消去法とピボット選択 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_06/

Pythonで学ぶ数値計算のアルゴリズムと実装 第8回：スパース性を活かす：帯行列の高速計算アルゴリズム | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerical_08/

Pythonで学ぶ数値計算のアルゴリズムと実装 第9回：最小二乗法とQR分解：ハウスホルダー変換による安定化 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerial_09/

Pythonで学ぶ数値計算のアルゴリズムと実装 第10回：非線形方程式の探索：2分法からニュートン法まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_10/

Pythonで学ぶ数値計算のアルゴリズムと実装 第11回：行列の本質を捉える：固有値問題の数値解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_11/

Pythonで学ぶ数値計算のアルゴリズムと実装 第12回：現象の変化を追う：ルンゲ・クッタ法による常微分方程式の解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_12/

Pythonで学ぶ数値計算のアルゴリズムと実装 第13回：空間の熱と振動を解く：偏微分方程式の差分近似と反復解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_13/