ガウス過程回帰で予測分布を計算するとき、あるいは多変量正規分布からサンプルを生成するとき、必ず登場する操作があります。それは「正定値対称行列の連立方程式を解く」ことです。一般のLU分解でも解けますが、対称性を活用すれば計算量を半分にでき、しかも数値的にも安定な方法が存在します。
この方法がコレスキー分解（Cholesky decomposition）です。正定値対称行列 $\bm{A}$ を下三角行列 $\bm{L}$ とその転置 $\bm{L}^T$ の積に分解します。
$$ \bm{A} = \bm{L}\bm{L}^T $$
コレスキー分解を理解すると、以下のような応用が開けます。
- ガウス過程: カーネル行列の分解によるベイズ予測の高速化
- 多変量正規分布: 共分散行列の分解によるサンプル生成
- 最適化: ニュートン法におけるヘッセ行列の連立方程式の効率的な求解
- カルマンフィルタ: 誤差共分散行列の更新
本記事の内容
- コレスキー分解の直感的な理解と数学的定義
- 分解が一意に存在するための条件（正定値性）の証明
- アルゴリズムの導出と計算量の解析
- Pythonによる実装と応用例
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 行列の基本演算 — 行列の積と転置の基礎
- 固有値と固有ベクトル — 正定値性の理解に必要
- LU分解 — 三角行列分解の基本
コレスキー分解とは — 直感的な理解
コレスキー分解を直感的に理解するために、まず「正定値対称行列」がどんなものかを考えましょう。
2次元の正定値対称行列 $\bm{A}$ は、2次形式 $\bm{x}^T\bm{A}\bm{x}$ を通じて楕円を定義します。この楕円は「ゆがんだ円」と見なせます。コレスキー分解 $\bm{A} = \bm{L}\bm{L}^T$ は、この楕円を「$\bm{L}$ による線形変換で単位円が写る先」として表現します。言い換えれば、$\bm{L}$ は単位円を楕円に変形する変換の「平方根」のようなものです。
もう一つのアナロジーとして、スカラーの場合を考えてみましょう。正の実数 $a > 0$ に対して $a = l^2$（$l > 0$）と書けますが、これはコレスキー分解の1次元版です。行列に拡張すると「下三角行列の平方根」になります。下三角に限定するのは、分解を一意にするためです。
歴史的には、フランスの陸軍士官アンドレ＝ルイ・コレスキー（1875-1918）が測量計算のためにこの方法を開発しました。彼は第一次世界大戦で戦死し、この方法は死後に同僚によって発表されました。実用的な計算手法として開発された背景があり、今日でも数値計算の現場で最もよく使われる行列分解の一つです。
コレスキー分解の直感を掴んだところで、次に正定値行列の定義とコレスキー分解の存在定理を見ていきましょう。
正定値対称行列の復習
コレスキー分解の前提条件である「正定値対称行列」を正確に定義しておきます。
$n \times n$ の実行列 $\bm{A}$ が正定値対称行列（symmetric positive definite, SPD）であるとは、以下の2条件を満たすことです。
- 対称性: $\bm{A} = \bm{A}^T$
- 正定値性: 全てのゼロでないベクトル $\bm{x} \in \mathbb{R}^n$ に対して $\bm{x}^T\bm{A}\bm{x} > 0$
正定値性は「行列が定義する2次形式が常に正の値をとる」ことを意味します。幾何学的には、2次形式の等高面が原点を中心とした閉じた楕円体になります。
正定値性の同値条件
正定値対称行列には多くの同値な特徴付けがあります。$\bm{A}$ が対称行列のとき、以下は全て同値です。
- 全ての固有値が正: $\lambda_i > 0$（$i = 1, \ldots, n$）
- 全ての主小行列式（leading principal minor）が正: $\det(\bm{A}_k) > 0$（$k = 1, \ldots, n$）
- 正の対角成分を持つ下三角行列 $\bm{L}$ が存在して $\bm{A} = \bm{L}\bm{L}^T$（コレスキー分解の存在）
- 正則行列 $\bm{B}$ が存在して $\bm{A} = \bm{B}^T\bm{B}$
最後の同値条件は、正定値行列が「ある行列の列ベクトル同士の内積（グラム行列）」として現れることを示しています。機械学習でのカーネル行列が正定値であるのは、まさにこの性質によるものです。
これらの性質を踏まえて、コレスキー分解の数学的な定義と存在証明に入りましょう。
コレスキー分解の数学的定義と存在定理
定義
正定値対称行列 $\bm{A}$ のコレスキー分解とは
$$ \begin{equation} \bm{A} = \bm{L}\bm{L}^T \end{equation} $$
を満たす、正の対角成分を持つ下三角行列 $\bm{L}$ を求めることです。$\bm{L}$ をコレスキー因子と呼びます。
この定義が「なぜ下三角行列なのか」を考えてみましょう。LU分解 $\bm{A} = \bm{L}\bm{U}$ では、$\bm{L}$（下三角）と $\bm{U}$（上三角）の2つの行列が必要でした。しかし $\bm{A}$ が対称なら $\bm{L}\bm{U} = \bm{A} = \bm{A}^T = \bm{U}^T\bm{L}^T$ であり、分解の一意性から $\bm{U} = \bm{D}\bm{L}^T$（$\bm{D}$ は対角行列）が導かれます。正定値性により $\bm{D}$ の対角成分は全て正なので、$\bm{D}^{1/2}$ を $\bm{L}$ に吸収すれば $\bm{A} = \bm{L}\bm{L}^T$ が得られるのです。
存在と一意性の証明
定理: $\bm{A}$ が $n \times n$ の正定値対称行列ならば、正の対角成分を持つ下三角行列 $\bm{L}$ が一意に存在して $\bm{A} = \bm{L}\bm{L}^T$ が成り立つ。
証明: $n$ に関する数学的帰納法で示します。
基底ケース ($n = 1$): $\bm{A} = (a)$ で $a > 0$（正定値性より）。$\bm{L} = (\sqrt{a})$ とすれば $\bm{L}\bm{L}^T = (a) = \bm{A}$。
帰納ステップ: $n-1$ 次の正定値対称行列に対してコレスキー分解が存在すると仮定します。$n$ 次の正定値対称行列 $\bm{A}$ をブロック分割します。
$$ \bm{A} = \begin{pmatrix} \bm{A}_{11} & \bm{a}_{12} \\ \bm{a}_{12}^T & a_{22} \end{pmatrix} $$
ここで $\bm{A}_{11}$ は $(n-1) \times (n-1)$、$\bm{a}_{12}$ は $(n-1) \times 1$、$a_{22}$ はスカラーです。
$\bm{A}$ が正定値であることから $\bm{A}_{11}$ も正定値です（$\bm{A}_{11}$ の任意の主小行列式が $\bm{A}$ の主小行列式に等しいため）。帰納法の仮定より、$\bm{A}_{11} = \bm{L}_{11}\bm{L}_{11}^T$ が存在します。
コレスキー因子を次のように仮定します。
$$ \bm{L} = \begin{pmatrix} \bm{L}_{11} & \bm{0} \\ \bm{l}_{21}^T & l_{22} \end{pmatrix} $$
$\bm{L}\bm{L}^T = \bm{A}$ を展開すると
$$ \begin{pmatrix} \bm{L}_{11}\bm{L}_{11}^T & \bm{L}_{11}\bm{l}_{21} \\ \bm{l}_{21}^T\bm{L}_{11}^T & \bm{l}_{21}^T\bm{l}_{21} + l_{22}^2 \end{pmatrix} = \begin{pmatrix} \bm{A}_{11} & \bm{a}_{12} \\ \bm{a}_{12}^T & a_{22} \end{pmatrix} $$
ブロックを比較すると
$(1, 2)$ ブロックから $\bm{L}_{11}\bm{l}_{21} = \bm{a}_{12}$ が得られます。$\bm{L}_{11}$ は正則（対角成分が正）なので、前進代入で $\bm{l}_{21}$ が一意に求まります。
$$ \bm{l}_{21} = \bm{L}_{11}^{-1}\bm{a}_{12} $$
$(2, 2)$ ブロックから $l_{22}^2 = a_{22} – \bm{l}_{21}^T\bm{l}_{21}$ が得られます。
ここで $a_{22} – \bm{l}_{21}^T\bm{l}_{21} > 0$ を示す必要があります。シュア補行列を考えると
$$ a_{22} – \bm{a}_{12}^T\bm{A}_{11}^{-1}\bm{a}_{12} = a_{22} – \bm{l}_{21}^T\bm{l}_{21} $$
$\bm{A}$ が正定値ならシュア補行列も正であることから（$\det \bm{A} = \det \bm{A}_{11} \cdot (a_{22} – \bm{a}_{12}^T\bm{A}_{11}^{-1}\bm{a}_{12}) > 0$ と $\det \bm{A}_{11} > 0$ より）、$l_{22} = \sqrt{a_{22} – \bm{l}_{21}^T\bm{l}_{21}} > 0$ が一意に定まります。
以上により帰納法が完了し、コレスキー分解の存在と一意性が証明されました。$\square$
この証明自体がアルゴリズムのレシピになっていることに注目してください。次にこれを計算手順として整理しましょう。
コレスキー分解のアルゴリズム
算出公式の導出
$\bm{A} = \bm{L}\bm{L}^T$ の成分を直接比較して、$\bm{L}$ の各成分を求める公式を導出します。
$\bm{A}$ の $(i, j)$ 成分を $a_{ij}$、$\bm{L}$ の $(i, j)$ 成分を $l_{ij}$ とします。$(\bm{L}\bm{L}^T)_{ij} = \sum_{k=1}^{n} l_{ik} l_{jk}$ ですが、$\bm{L}$ は下三角なので $k > \min(i, j)$ では $l_{ik}$ または $l_{jk}$ がゼロになります。
対角成分 ($i = j$):
$$ a_{ii} = \sum_{k=1}^{i} l_{ik}^2 = l_{ii}^2 + \sum_{k=1}^{i-1} l_{ik}^2 $$
$l_{ii}$ について解くと
$$ \begin{equation} l_{ii} = \sqrt{a_{ii} – \sum_{k=1}^{i-1} l_{ik}^2} \end{equation} $$
非対角成分 ($i > j$):
$$ a_{ij} = \sum_{k=1}^{j} l_{ik} l_{jk} = l_{ij} l_{jj} + \sum_{k=1}^{j-1} l_{ik} l_{jk} $$
$l_{ij}$ について解くと
$$ \begin{equation} l_{ij} = \frac{1}{l_{jj}}\left(a_{ij} – \sum_{k=1}^{j-1} l_{ik} l_{jk}\right) \end{equation} $$
これらの公式から、$j = 1, 2, \ldots, n$ の順に $\bm{L}$ の第 $j$ 列を上から下へ計算していけばよいことがわかります。
計算量の解析
コレスキー分解の計算量を見積もりましょう。第 $j$ 列の計算で
- 対角成分 $l_{jj}$: $j-1$ 回の乗算 + 1回の平方根
- 非対角成分 $l_{ij}$（$i = j+1, \ldots, n$）: 各 $j-1$ 回の乗算 × $(n-j)$ 個
全体の乗算回数は
$$ \sum_{j=1}^{n} \left[(j-1) + (n-j)(j-1)\right] = \sum_{j=1}^{n} (j-1)(n-j+1) \approx \frac{n^3}{6} $$
LU分解の計算量 $n^3/3$ と比較してちょうど半分です。これは対称性を活用している恩恵です。
LU分解との比較
| 特性 | コレスキー分解 | LU分解 |
|---|---|---|
| 適用条件 | 正定値対称行列 | 一般の正則行列 |
| 計算量 | $\frac{n^3}{6}$ | $\frac{n^3}{3}$ |
| メモリ | $\frac{n(n+1)}{2}$（下三角のみ） | $n^2$（$\bm{L}$ と $\bm{U}$） |
| ピボット選択 | 不要 | 必要（部分ピボット） |
| 数値安定性 | 無条件に安定 | ピボット選択で安定化 |
コレスキー分解はピボット選択が不要です。これは正定値性により、分解過程で現れる対角要素が常に正であることが保証されるためです。ゼロ除算や大きな要素の発生がないので、追加の安定化処理なしに数値的に安定な分解が得られます。
理論面の準備が整ったので、次にPythonで実装して動作を確認しましょう。
Pythonでの実装
スクラッチ実装
まず、上で導出した公式に基づいてコレスキー分解をスクラッチで実装します。
import numpy as np
import matplotlib.pyplot as plt
def cholesky_decomposition(A):
"""コレスキー分解のスクラッチ実装"""
n = A.shape[0]
L = np.zeros_like(A, dtype=float)
for j in range(n):
# 対角成分
s = A[j, j] - np.sum(L[j, :j]**2)
if s <= 0:
raise ValueError(f"行列が正定値でありません (ステップ {j}, s={s:.2e})")
L[j, j] = np.sqrt(s)
# 非対角成分（第j列のj+1行目以降）
for i in range(j + 1, n):
L[i, j] = (A[i, j] - np.sum(L[i, :j] * L[j, :j])) / L[j, j]
return L
# テスト: 3x3の正定値対称行列
A = np.array([[4, 2, 1],
[2, 5, 3],
[1, 3, 6]], dtype=float)
L = cholesky_decomposition(A)
print("A =")
print(A)
print("\nL (コレスキー因子) =")
print(np.round(L, 6))
print("\nL L^T (復元) =")
print(np.round(L @ L.T, 10))
print("\n復元誤差 ||A - LL^T|| =", np.linalg.norm(A - L @ L.T))
# NumPyの結果と比較
L_np = np.linalg.cholesky(A)
print("\nNumPyとの差 ||L - L_np|| =", np.linalg.norm(L - L_np))
上のコードでは、3行3列の正定値対称行列に対してコレスキー分解を実行しています。スクラッチ実装の結果がNumPyの np.linalg.cholesky と一致すること、および $\bm{L}\bm{L}^T$ が元の行列 $\bm{A}$ を正確に復元することを確認しています。分解過程で対角要素が負になった場合は正定値でないことを示すエラーが発生する安全機構も組み込んでいます。
コレスキー分解による連立方程式の解法
コレスキー分解の最も基本的な応用として、正定値対称行列の連立方程式 $\bm{A}\bm{x} = \bm{b}$ を解きましょう。
import numpy as np
def solve_cholesky(A, b):
"""コレスキー分解による連立方程式の解法"""
L = np.linalg.cholesky(A)
# 前進代入: Ly = b
y = np.linalg.solve(L, b)
# 後退代入: L^T x = y
x = np.linalg.solve(L.T, y)
return x
# テスト
n = 5
np.random.seed(42)
B = np.random.randn(n, n)
A = B.T @ B + 0.1 * np.eye(n) # 正定値対称行列
b = np.random.randn(n)
x_chol = solve_cholesky(A, b)
x_direct = np.linalg.solve(A, b)
print("コレスキー法の解:", np.round(x_chol, 8))
print("直接解法の解: ", np.round(x_direct, 8))
print("差のノルム:", np.linalg.norm(x_chol - x_direct))
コレスキー分解を使った連立方程式の解法は、$\bm{A}\bm{x} = \bm{b}$ を2段階に分けます。まず $\bm{L}\bm{y} = \bm{b}$ を前進代入で解き、次に $\bm{L}^T\bm{x} = \bm{y}$ を後退代入で解きます。どちらも三角行列の方程式なので $O(n^2)$ で解けます。
応用1: 多変量正規分布からのサンプル生成
コレスキー分解の重要な応用として、多変量正規分布 $\mathcal{N}(\bm{\mu}, \bm{\Sigma})$ からのサンプル生成があります。
標準正規分布のサンプル $\bm{z} \sim \mathcal{N}(\bm{0}, \bm{I})$ を $\bm{x} = \bm{\mu} + \bm{L}\bm{z}$ と変換すると、$\bm{x}$ の共分散行列は $E[(\bm{x} – \bm{\mu})(\bm{x} – \bm{\mu})^T] = \bm{L}E[\bm{z}\bm{z}^T]\bm{L}^T = \bm{L}\bm{L}^T = \bm{\Sigma}$ となります。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# 2次元の共分散行列
mu = np.array([2, 3])
Sigma = np.array([[2.0, 1.2],
[1.2, 1.5]])
# コレスキー分解
L = np.linalg.cholesky(Sigma)
# サンプル生成
n_samples = 2000
z = np.random.randn(2, n_samples) # 標準正規
x = mu[:, None] + L @ z # 多変量正規
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) 標準正規 → 多変量正規への変換
ax = axes[0]
ax.scatter(z[0, :500], z[1, :500], s=5, alpha=0.4, color="blue", label="Standard normal z")
ax.scatter(x[0, :500], x[1, :500], s=5, alpha=0.4, color="red", label="Transformed x = μ + Lz")
ax.set_xlabel("$x_1$", fontsize=12)
ax.set_ylabel("$x_2$", fontsize=12)
ax.set_title("Cholesky Transform", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
ax.set_xlim(-5, 7)
ax.set_ylim(-4, 8)
# (b) サンプルの分布（2D密度）
ax = axes[1]
ax.scatter(x[0], x[1], s=3, alpha=0.3, color="steelblue")
# 理論的な楕円（1σ, 2σ, 3σ）
theta = np.linspace(0, 2*np.pi, 100)
circle = np.array([np.cos(theta), np.sin(theta)])
for k, color in zip([1, 2, 3], ["red", "orange", "yellow"]):
ellipse = mu[:, None] + k * L @ circle
ax.plot(ellipse[0], ellipse[1], color=color, linewidth=2,
label=f"{k}σ ellipse")
ax.set_xlabel("$x_1$", fontsize=12)
ax.set_ylabel("$x_2$", fontsize=12)
ax.set_title("Multivariate Normal Samples", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
# (c) サンプルの統計量の確認
ax = axes[2]
sample_mean = np.mean(x, axis=1)
sample_cov = np.cov(x)
text = (f"Theory:\n"
f" μ = [{mu[0]:.1f}, {mu[1]:.1f}]\n"
f" Σ = [[{Sigma[0,0]:.1f}, {Sigma[0,1]:.1f}],\n"
f" [{Sigma[1,0]:.1f}, {Sigma[1,1]:.1f}]]\n\n"
f"Sample (n={n_samples}):\n"
f" μ̂ = [{sample_mean[0]:.3f}, {sample_mean[1]:.3f}]\n"
f" Σ̂ = [[{sample_cov[0,0]:.3f}, {sample_cov[0,1]:.3f}],\n"
f" [{sample_cov[1,0]:.3f}, {sample_cov[1,1]:.3f}]]\n\n"
f"Cholesky factor:\n"
f" L = [[{L[0,0]:.4f}, 0],\n"
f" [{L[1,0]:.4f}, {L[1,1]:.4f}]]")
ax.text(0.05, 0.95, text, transform=ax.transAxes, fontsize=10,
verticalalignment="top", fontfamily="monospace",
bbox=dict(boxstyle="round", facecolor="lightyellow", alpha=0.8))
ax.axis("off")
ax.set_title("Statistics Comparison", fontsize=13)
plt.tight_layout()
plt.savefig("cholesky_mvn.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、コレスキー分解による多変量正規分布のサンプル生成の仕組みが読み取れます。
-
左図（変換の可視化）: 青い点が標準正規分布 $\mathcal{N}(\bm{0}, \bm{I})$ のサンプル、赤い点がコレスキー変換 $\bm{x} = \bm{\mu} + \bm{L}\bm{z}$ 後のサンプルです。等方的な円形の分布が、共分散構造を持つ楕円形の分布に変換されていることが確認できます。
-
中央図（サンプルと信頼楕円）: 生成されたサンプル（青い点）の上に、理論的な1σ, 2σ, 3σの信頼楕円を重ねています。サンプルが楕円の範囲内に理論通りの割合で分布していることがわかります。楕円の傾きは共分散行列の非対角成分（相関）を反映しています。
-
右図（統計量の比較）: サンプルの平均と共分散が理論値に近いことを数値で確認しています。2000サンプルでの推定値が理論値によく一致しており、コレスキー変換が正しく機能していることがわかります。
応用2: 計算速度の比較
正定値対称行列の連立方程式に対して、コレスキー分解がLU分解やnumpyの一般解法と比べてどの程度高速かを測定します。
import numpy as np
import matplotlib.pyplot as plt
import time
def benchmark_solvers(n, n_trials=5):
"""n次元の正定値対称行列で各手法の計算時間を比較"""
np.random.seed(42)
B = np.random.randn(n, n)
A = B.T @ B + 0.01 * np.eye(n)
b = np.random.randn(n)
times = {}
# numpy.linalg.solve（一般）
t_list = []
for _ in range(n_trials):
t0 = time.perf_counter()
np.linalg.solve(A, b)
t_list.append(time.perf_counter() - t0)
times["linalg.solve"] = np.median(t_list)
# コレスキー分解
t_list = []
for _ in range(n_trials):
t0 = time.perf_counter()
L = np.linalg.cholesky(A)
y = np.linalg.solve(L, b)
np.linalg.solve(L.T, y)
t_list.append(time.perf_counter() - t0)
times["Cholesky"] = np.median(t_list)
# scipy.linalg.cho_solve
from scipy.linalg import cho_factor, cho_solve
t_list = []
for _ in range(n_trials):
t0 = time.perf_counter()
c, low = cho_factor(A)
cho_solve((c, low), b)
t_list.append(time.perf_counter() - t0)
times["scipy cho_solve"] = np.median(t_list)
return times
sizes = [50, 100, 200, 500, 1000, 2000]
results = {method: [] for method in ["linalg.solve", "Cholesky", "scipy cho_solve"]}
for n in sizes:
times = benchmark_solvers(n)
for method in results:
results[method].append(times[method])
fig, ax = plt.subplots(figsize=(10, 6))
for method, marker, color in [("linalg.solve", "o", "red"),
("Cholesky", "s", "blue"),
("scipy cho_solve", "^", "green")]:
ax.loglog(sizes, results[method], f"{color[0]}{marker}-", linewidth=2,
markersize=8, label=method, color=color)
# O(n^3) 参照線
ref = np.array(sizes, dtype=float)
ax.loglog(ref, 1e-7 * ref**3, "k:", alpha=0.3, label="$O(n^3)$ reference")
ax.set_xlabel("Matrix size n", fontsize=13)
ax.set_ylabel("Computation time (sec)", fontsize=13)
ax.set_title("Solver Performance: SPD Linear System", fontsize=14)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which="both")
plt.tight_layout()
plt.savefig("cholesky_benchmark.png", dpi=150, bbox_inches="tight")
plt.show()
このベンチマーク結果から、コレスキー分解の計算効率の優位性が確認できます。
-
全ての行列サイズでコレスキー分解が最速: scipy の
cho_solveは内部でFortranの最適化されたLAPACKルーチンを呼ぶため、特に大規模行列で一般的なlinalg.solveより大きなアドバンテージがあります。 -
スケーリングは全て $O(n^3)$: いずれの方法も計算量のオーダーは $O(n^3)$ ですが、定数係数がコレスキー分解の方が約2倍小さいため、同じ $n$ に対して約2倍高速になっています。
-
$n$ が大きいほど差が顕著: $n = 2000$ 程度になると計算時間の絶対差が大きくなり、正定値性を活用することの実用的な意義が明確になります。
不完全コレスキー分解と正定値でない場合の対処
不完全コレスキー分解
大規模な疎行列では、完全なコレスキー分解の因子 $\bm{L}$ がフィルイン（ゼロだった位置に非ゼロ要素が出現する現象）により密行列になってしまうことがあります。不完全コレスキー分解（incomplete Cholesky, IC）は、特定のスパースパターンを維持するよう制約を加えた近似的なコレスキー分解です。前処理付き共役勾配法（PCG）の前処理行列として広く使われています。
最もシンプルな変種である IC(0) は、元の行列 $\bm{A}$ と同じスパースパターンだけを因子 $\bm{L}$ に保持します。つまり $\bm{A}$ でゼロの位置は $\bm{L}$ でもゼロのままに強制します。こうすることでメモリ使用量が元の行列と同程度に抑えられ、大規模な有限要素法の連立方程式などでも実用的な前処理が可能になります。より高精度が必要な場合は、一定レベルのフィルインを許容する IC($p$) も利用されます。
正定値でない場合の対処
実際のアプリケーションでは、丸め誤差により理論的には正定値であるはずの行列が数値的に正定値でなくなることがあります。この場合の対処法として
- ジッタの追加: $\bm{A} + \epsilon \bm{I}$（$\epsilon$ は小さな正の数）として正定値性を回復する。ガウス過程でよく使われる手法です。
- 修正コレスキー分解: 分解過程で負の対角要素が出現した場合に、それを正の値に置き換える方法（Gill-Murray修正など）。最適化における非正定値ヘッシアンの処理に使われます。
なお、コレスキー分解が途中で失敗する（対角要素の平方根の中身が負になる）こと自体を「行列が正定値でないことの検出」に利用するテクニックもあります。固有値を全て計算するよりもはるかに低コストで正定値性を判定できるため、行列が正定値かどうかを事前チェックする実用的な手段として広く使われています。
このように、コレスキー分解は正定値対称行列という構造を最大限に活かした分解法であり、理論的な美しさと実用的な効率性を兼ね備えています。ガウス過程回帰ではカーネル行列の分解、カルマンフィルタでは共分散行列の更新、ベイズ推定では事後分布の共分散行列の操作と、科学技術計算のあらゆる場面で中心的な役割を果たします。正定値対称行列が現れたら、まずコレスキー分解を検討するのが数値計算における基本的な指針です。さらに、コレスキー分解を一度計算しておけば、右辺ベクトル $\bm{b}$ が異なる複数の連立方程式を追加コスト $O(n^2)$ で解けるという利点もあります。これはガウス過程回帰で超パラメータの勾配を計算する際など、同一のカーネル行列に対して異なる右辺を繰り返し解く場面で特に有効です。
まとめ
本記事では、コレスキー分解の理論と実装について解説しました。
- コレスキー分解は正定値対称行列 $\bm{A}$ を $\bm{L}\bm{L}^T$ に分解する手法であり、対称性を活用してLU分解の半分の計算量で実行できる
- 存在と一意性は数学的帰納法で証明され、分解が可能であるための必要十分条件は行列の正定値性である
- 数値安定性に優れ、ピボット選択が不要であるため実装が簡潔になる
- 多変量正規分布のサンプル生成では、$\bm{x} = \bm{\mu} + \bm{L}\bm{z}$ の変換により効率的にサンプルが得られる
- 連立方程式の解法では、前進代入と後退代入の2段階で解が求まり、一般の解法より約2倍高速である
次のステップとして、以下の記事も参考にしてください。
- QR分解とグラム・シュミット法 — 直交分解の別の手法
- ムーア・ペンローズ擬逆行列の理論 — SVDによる一般的な逆行列
- ガウス過程回帰 — コレスキー分解の重要な応用先