3つのデータ点 $(1, 2)$, $(2, 3)$, $(3, 5)$ を通る直線 $y = ax + b$ を求めたいとします。連立方程式を立てると
$$ \begin{pmatrix} 1 & 1 \\ 2 & 1 \\ 3 & 1 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} 2 \\ 3 \\ 5 \end{pmatrix} $$
3本の方程式に対して未知数は2個。一般にこの連立方程式は解を持ちません（方程式が過剰決定）。しかし「最もよい近似解」は求められます。この「最もよい」とは何を意味し、どう計算するのでしょうか。
答えは射影行列（projection matrix）にあります。射影行列は、あるベクトルを部分空間に「垂直に落とす」操作を表す行列です。最小二乗法は、右辺ベクトル $\bm{b}$ を列空間に射影する操作そのものであり、射影行列を理解すれば最小二乗法の本質が幾何学的に見えてきます。
射影行列と最小二乗法を理解すると、以下のような場面で活用できます。
- 線形回帰: 統計学・機械学習における最も基本的なモデル
- 信号処理: ノイズのある信号から最良の推定を得るフィルタリング
- 数値解析: 過剰決定系の近似解法
- コンピュータグラフィックス: 3次元点の2次元平面への射影
本記事の内容
- ベクトルの射影と射影行列の定義
- 射影行列の性質（冪等性・対称性）
- 最小二乗法の幾何学的導出
- 正規方程式とその導出
- 残差の性質と分散分析
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 固有値分解の理論と応用 — 固有値の概念
- 行列式の定義・計算方法・幾何学的意味 — 行列の正則性
ベクトルの射影
1次元の射影（直線への射影）
最も簡単な場合から始めましょう。ベクトル $\bm{b}$ を直線の方向 $\bm{a}$（$\bm{a} \neq \bm{0}$）に射影することを考えます。
射影ベクトル $\hat{\bm{b}}$ は $\bm{a}$ のスカラー倍 $\hat{\bm{b}} = x\bm{a}$ で、残差 $\bm{e} = \bm{b} – \hat{\bm{b}}$ が $\bm{a}$ に直交する条件から求まります。
$$ \bm{a}^T(\bm{b} – x\bm{a}) = 0 $$
$x$ について解くと
$$ x = \frac{\bm{a}^T\bm{b}}{\bm{a}^T\bm{a}} $$
したがって射影ベクトルは
$$ \hat{\bm{b}} = \frac{\bm{a}^T\bm{b}}{\bm{a}^T\bm{a}} \bm{a} = \frac{\bm{a}\bm{a}^T}{\bm{a}^T\bm{a}} \bm{b} $$
ここで $\bm{P} = \frac{\bm{a}\bm{a}^T}{\bm{a}^T\bm{a}}$ が直線 $\bm{a}$ 方向への射影行列です。
具体的に $\bm{a} = (1, 1)^T$ のとき
$$ \bm{P} = \frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} $$
$\bm{b} = (3, 1)^T$ を射影すると $\hat{\bm{b}} = \bm{Pb} = (2, 2)^T$ となり、残差 $\bm{e} = (1, -1)^T$ は確かに $\bm{a}$ に直交します（$\bm{a}^T\bm{e} = 1 \cdot 1 + 1 \cdot (-1) = 0$）。
部分空間への射影
次に、直線ではなく部分空間（平面や一般の線形部分空間）への射影を考えます。
$\bm{A}$ を $m \times n$ 行列（$m > n$、列が線形独立）とし、$\bm{A}$ の列が張る部分空間を $C(\bm{A})$（列空間）とします。ベクトル $\bm{b} \in \mathbb{R}^m$ を $C(\bm{A})$ に射影したベクトル $\hat{\bm{b}}$ は
$$ \hat{\bm{b}} = \bm{A}\hat{\bm{x}} $$
の形で書け、残差 $\bm{e} = \bm{b} – \bm{A}\hat{\bm{x}}$ が $C(\bm{A})$ に直交する条件
$$ \bm{A}^T(\bm{b} – \bm{A}\hat{\bm{x}}) = \bm{0} $$
から $\hat{\bm{x}}$ が決まります。これを整理すると
$$ \begin{equation} \bm{A}^T\bm{A}\hat{\bm{x}} = \bm{A}^T\bm{b} \end{equation} $$
これが正規方程式（normal equations）です。$\bm{A}$ の列が線形独立なら $\bm{A}^T\bm{A}$ は正則なので
$$ \hat{\bm{x}} = (\bm{A}^T\bm{A})^{-1}\bm{A}^T\bm{b} $$
射影ベクトルは
$$ \hat{\bm{b}} = \bm{A}\hat{\bm{x}} = \bm{A}(\bm{A}^T\bm{A})^{-1}\bm{A}^T\bm{b} $$
したがって列空間 $C(\bm{A})$ への射影行列は
$$ \begin{equation} \bm{P} = \bm{A}(\bm{A}^T\bm{A})^{-1}\bm{A}^T \end{equation} $$
1次元の場合（$\bm{A} = \bm{a}$ がベクトル）に戻すと $\bm{P} = \frac{\bm{a}\bm{a}^T}{\bm{a}^T\bm{a}}$ となり、先ほどの公式と一致します。
射影の定義がわかったところで、射影行列の代数的性質を調べましょう。
射影行列の性質
冪等性
射影行列の最も重要な性質は冪等性（idempotency）です。
$$ \begin{equation} \bm{P}^2 = \bm{P} \end{equation} $$
直感的には「一度射影したベクトルをもう一度射影しても変わらない」ことを意味します。部分空間上のベクトルを同じ部分空間に射影しても、そのベクトル自身が返ってくるのは自然です。
代数的に確認すると
$$ \bm{P}^2 = \bm{A}(\bm{A}^T\bm{A})^{-1}\bm{A}^T \cdot \bm{A}(\bm{A}^T\bm{A})^{-1}\bm{A}^T = \bm{A}(\bm{A}^T\bm{A})^{-1}\bm{A}^T = \bm{P} $$
途中で $\bm{A}^T\bm{A} \cdot (\bm{A}^T\bm{A})^{-1} = \bm{I}$ を使いました。
対称性
直交射影行列は対称です。
$$ \begin{equation} \bm{P}^T = \bm{P} \end{equation} $$
これは $(\bm{A}(\bm{A}^T\bm{A})^{-1}\bm{A}^T)^T = \bm{A}((\bm{A}^T\bm{A})^{-1})^T\bm{A}^T = \bm{A}((\bm{A}^T\bm{A})^T)^{-1}\bm{A}^T = \bm{A}(\bm{A}^T\bm{A})^{-1}\bm{A}^T$ から確認できます（$\bm{A}^T\bm{A}$ は対称なので $(\bm{A}^T\bm{A})^T = \bm{A}^T\bm{A}$）。
逆に、$\bm{P}^2 = \bm{P}$ かつ $\bm{P}^T = \bm{P}$ を満たす行列は何らかの部分空間への直交射影行列です。
固有値
冪等行列の固有値は0または1のみです。なぜなら $\bm{P}\bm{x} = \lambda\bm{x}$ に左から $\bm{P}$ を掛けると $\bm{P}^2\bm{x} = \lambda\bm{Px}$ つまり $\bm{Px} = \lambda^2\bm{x}$。$\lambda\bm{x} = \lambda^2\bm{x}$ なので $\lambda(\lambda – 1) = 0$、したがって $\lambda = 0$ または $\lambda = 1$ です。
固有値1に対応する固有ベクトルは射影先の部分空間に属するベクトル、固有値0に対応する固有ベクトルは直交補空間に属するベクトルです。
トレースとランク
射影行列のトレースは射影先の部分空間の次元に等しいです。
$$ \text{tr}(\bm{P}) = \text{rank}(\bm{P}) = \text{rank}(\bm{A}) $$
これは固有値が0と1のみで、1の個数がランク（= 射影先の次元）に等しいことから従います。
補射影行列
$\bm{P}$ が部分空間 $V$ への射影なら、$\bm{I} – \bm{P}$ は直交補空間 $V^\perp$ への射影です。
$$ (\bm{I} – \bm{P})^2 = \bm{I} – 2\bm{P} + \bm{P}^2 = \bm{I} – 2\bm{P} + \bm{P} = \bm{I} – \bm{P} $$
任意のベクトル $\bm{b}$ は $\bm{b} = \bm{Pb} + (\bm{I} – \bm{P})\bm{b}$ と直交分解され、$\bm{Pb} \in V$ と $(\bm{I} – \bm{P})\bm{b} \in V^\perp$ は直交します。
射影行列の性質を整理したところで、これが最小二乗法とどうつながるかを見ましょう。
最小二乗法の幾何学的理解
問題設定
$m \times n$ 行列 $\bm{A}$（$m > n$、$\text{rank}(\bm{A}) = n$）と $\bm{b} \in \mathbb{R}^m$ に対して、$\bm{Ax} = \bm{b}$ は一般に解を持ちません。最小二乗法（least squares method）は、残差の二乗和を最小化する $\hat{\bm{x}}$ を求めます。
$$ \hat{\bm{x}} = \arg\min_{\bm{x}} \|\bm{b} – \bm{Ax}\|^2 $$
幾何学的解釈
$\bm{Ax}$ は $\bm{x}$ を変化させると $\bm{A}$ の列空間 $C(\bm{A})$ の全てのベクトルを生成します。したがって「$\bm{b}$ に最も近い $C(\bm{A})$ 内のベクトル」を求めることが最小二乗問題です。
幾何学的に明らかなのは、$\bm{b}$ に最も近い $C(\bm{A})$ 内の点は $\bm{b}$ の $C(\bm{A})$ への直交射影 $\hat{\bm{b}} = \bm{Pb}$ です。三平方の定理（ピタゴラスの定理）から
$$ \|\bm{b} – \bm{Ax}\|^2 = \|\bm{b} – \hat{\bm{b}}\|^2 + \|\hat{\bm{b}} – \bm{Ax}\|^2 $$
第1項は $\bm{x}$ に依存しないので、全体を最小化するには第2項 $\|\hat{\bm{b}} – \bm{Ax}\|^2 = 0$、つまり $\bm{A}\hat{\bm{x}} = \hat{\bm{b}}$ とすればよいです。
この $\hat{\bm{x}}$ がまさに正規方程式 $\bm{A}^T\bm{A}\hat{\bm{x}} = \bm{A}^T\bm{b}$ の解です。
微分による導出
幾何学的な方法とは別に、目的関数を微分する方法でも同じ結果が得られます。
$$ J(\bm{x}) = \|\bm{b} – \bm{Ax}\|^2 = (\bm{b} – \bm{Ax})^T(\bm{b} – \bm{Ax}) = \bm{b}^T\bm{b} – 2\bm{x}^T\bm{A}^T\bm{b} + \bm{x}^T\bm{A}^T\bm{Ax} $$
$\bm{x}$ で微分して0とおくと
$$ \frac{\partial J}{\partial \bm{x}} = -2\bm{A}^T\bm{b} + 2\bm{A}^T\bm{A}\bm{x} = \bm{0} $$
整理すると正規方程式 $\bm{A}^T\bm{A}\bm{x} = \bm{A}^T\bm{b}$ が得られます。
ヘッセ行列は $\frac{\partial^2 J}{\partial \bm{x}^2} = 2\bm{A}^T\bm{A}$ です。$\bm{A}$ がフルランクのとき $\bm{A}^T\bm{A}$ は正定値なので、この臨界点は確かに最小値を与えます。
2つの異なる方法（幾何と微分）で同じ正規方程式に到達しました。次に、最小二乗法の重要な性質を見ていきましょう。
残差の性質
残差ベクトル
最小二乗解 $\hat{\bm{x}}$ に対する残差ベクトルは
$$ \bm{e} = \bm{b} – \bm{A}\hat{\bm{x}} = \bm{b} – \bm{Pb} = (\bm{I} – \bm{P})\bm{b} $$
残差は $\bm{b}$ の列空間の直交補空間への射影です。
直交性
残差の最も重要な性質は、$\bm{A}$ の列空間に直交することです。
$$ \bm{A}^T\bm{e} = \bm{A}^T(\bm{b} – \bm{A}\hat{\bm{x}}) = \bm{A}^T\bm{b} – \bm{A}^T\bm{A}\hat{\bm{x}} = \bm{0} $$
これは正規方程式の直接の帰結です。幾何学的には、射影の残差は射影先の部分空間に垂直であることを意味します。
二乗和の分解
$$ \|\bm{b}\|^2 = \|\hat{\bm{b}}\|^2 + \|\bm{e}\|^2 $$
これは $\hat{\bm{b}} \perp \bm{e}$ からのピタゴラスの定理です。
統計学では
- $\|\bm{b}\|^2$: 全変動（SST: Total Sum of Squares）
- $\|\hat{\bm{b}}\|^2$: 回帰変動（SSR: Regression Sum of Squares）
- $\|\bm{e}\|^2$: 残差変動（SSE: Error Sum of Squares）
として分散分析（ANOVA）に使われます。決定係数 $R^2 = 1 – \|\bm{e}\|^2 / \|\bm{b}\|^2 = \|\hat{\bm{b}}\|^2 / \|\bm{b}\|^2$ は、モデルがデータの変動をどれだけ説明しているかを表します。
ただし、ここで $\bm{b}$ は中心化済み（平均0）を仮定しています。切片項がある場合は適切に修正が必要です。
ハット行列
射影行列 $\bm{P} = \bm{A}(\bm{A}^T\bm{A})^{-1}\bm{A}^T$ は統計学ではハット行列（hat matrix）$\bm{H}$ とも呼ばれます。$\hat{\bm{b}} = \bm{Hb}$ と書けるため、$\bm{b}$ に「ハットを被せる」行列です。
ハット行列の対角成分 $h_{ii}$ はレバレッジ（leverage）と呼ばれ、$0 \leq h_{ii} \leq 1$ を満たします。$h_{ii}$ が大きいデータ点は、回帰直線に対して大きな影響力を持つ「てこ」のような点です。
$$ 0 \leq h_{ii} \leq 1, \quad \sum_{i=1}^{m} h_{ii} = n = \text{rank}(\bm{A}) $$
残差の性質が理解できたところで、Pythonで理論を実装・可視化しましょう。
Pythonでの実装と可視化
射影の幾何学的可視化
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(16, 5))
# (a) 2次元での直線への射影
ax = fig.add_subplot(131)
a = np.array([2, 1])
b = np.array([1, 3])
# 射影行列
P_line = np.outer(a, a) / np.dot(a, a)
b_hat = P_line @ b
e = b - b_hat
# 直線
t = np.linspace(-1, 3, 100)
ax.plot(t * a[0] / np.linalg.norm(a), t * a[1] / np.linalg.norm(a) * np.linalg.norm(a),
'b-', linewidth=1.5, alpha=0.5, label='Line $\\mathrm{span}(a)$')
ax.plot([0, a[0]*1.5], [0, a[1]*1.5], 'b-', linewidth=1.5, alpha=0.5)
ax.plot([0, -a[0]*0.5], [0, -a[1]*0.5], 'b-', linewidth=1.5, alpha=0.5)
# ベクトル
ax.annotate('', xy=b, xytext=[0,0],
arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax.annotate('', xy=b_hat, xytext=[0,0],
arrowprops=dict(arrowstyle='->', color='green', lw=2))
ax.plot([b[0], b_hat[0]], [b[1], b_hat[1]], 'k--', linewidth=1.5)
# 直角マーク
corner_size = 0.15
e_dir = e / np.linalg.norm(e)
a_dir = a / np.linalg.norm(a)
corner = b_hat + corner_size * e_dir
ax.plot([corner[0], corner[0] + corner_size * a_dir[0]],
[corner[1], corner[1] + corner_size * a_dir[1]], 'k-', linewidth=1)
ax.plot([b_hat[0] + corner_size * a_dir[0], corner[0] + corner_size * a_dir[0]],
[b_hat[1] + corner_size * a_dir[1], corner[1] + corner_size * a_dir[1]], 'k-', linewidth=1)
ax.text(b[0]+0.1, b[1]+0.1, '$b$', fontsize=14, color='red')
ax.text(b_hat[0]+0.1, b_hat[1]-0.3, '$\\hat{b}=Pb$', fontsize=12, color='green')
ax.text((b[0]+b_hat[0])/2+0.1, (b[1]+b_hat[1])/2, '$e$', fontsize=12)
ax.set_xlim(-0.5, 3)
ax.set_ylim(-0.5, 3.5)
ax.set_aspect('equal')
ax.set_title('Projection onto a Line', fontsize=13)
ax.grid(True, alpha=0.3)
ax.legend(fontsize=9, loc='lower right')
# (b) 3次元での平面への射影
ax = fig.add_subplot(132, projection='3d')
A = np.array([[1, 0], [0, 1], [1, 1]], dtype=float)
b_3d = np.array([1, 2, 4], dtype=float)
P_plane = A @ np.linalg.inv(A.T @ A) @ A.T
b_hat_3d = P_plane @ b_3d
e_3d = b_3d - b_hat_3d
# 列空間（平面）
s = np.linspace(-1, 2, 10)
t_grid = np.linspace(-1, 2, 10)
S, T = np.meshgrid(s, t_grid)
X_plane = A[0, 0] * S + A[0, 1] * T
Y_plane = A[1, 0] * S + A[1, 1] * T
Z_plane = A[2, 0] * S + A[2, 1] * T
ax.plot_surface(X_plane, Y_plane, Z_plane, alpha=0.2, color='blue')
# ベクトル
ax.quiver(0, 0, 0, b_3d[0], b_3d[1], b_3d[2], color='red', arrow_length_ratio=0.1, linewidth=2)
ax.quiver(0, 0, 0, b_hat_3d[0], b_hat_3d[1], b_hat_3d[2], color='green', arrow_length_ratio=0.1, linewidth=2)
ax.plot([b_3d[0], b_hat_3d[0]], [b_3d[1], b_hat_3d[1]], [b_3d[2], b_hat_3d[2]],
'k--', linewidth=1.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
ax.set_title('Projection onto Column Space', fontsize=13)
# (c) 射影行列の冪等性
ax = fig.add_subplot(133)
n_test = 50
A_test = np.random.randn(n_test, 3)
P_test = A_test @ np.linalg.inv(A_test.T @ A_test) @ A_test.T
P2_test = P_test @ P_test
# P^2 - P の要素をヒストグラム
diff = (P2_test - P_test).flatten()
ax.hist(diff, bins=50, color='steelblue', edgecolor='black', alpha=0.7)
ax.axvline(x=0, color='red', linewidth=2, linestyle='--')
ax.set_xlabel('$(P^2 - P)_{ij}$', fontsize=12)
ax.set_ylabel('Count', fontsize=12)
ax.set_title(f'Idempotency Check ($P^2 = P$)\nmax $|P^2-P|$ = {np.max(np.abs(diff)):.2e}', fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('projection_matrix.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、射影行列の幾何学的意味が確認できます。
-
左図（直線への射影）: ベクトル $\bm{b}$（赤）が直線（青）上の点 $\hat{\bm{b}} = \bm{Pb}$（緑）に射影され、残差 $\bm{e}$（黒破線）が直線に直交しています。直角マークが直交性を視覚的に示しています
-
中央図（平面への射影）: 3次元空間のベクトル $\bm{b}$（赤）が、行列 $\bm{A}$ の列空間（青い平面）上の点 $\hat{\bm{b}}$（緑）に射影されています。残差ベクトル（黒破線）は平面に垂直です。これが最小二乗法の幾何学的な姿です
-
右図（冪等性の検証）: $\bm{P}^2 – \bm{P}$ の全要素のヒストグラムが0に集中しており、$\bm{P}^2 = \bm{P}$（冪等性）が数値的に確認されています。最大誤差が $10^{-14}$ のオーダーで、機械精度の範囲です
最小二乗法の実装と可視化
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# データ生成
np.random.seed(42)
n_data = 30
x_data = np.linspace(0, 5, n_data)
y_true = 2.0 * x_data + 1.0
y_data = y_true + np.random.randn(n_data) * 1.5
# デザイン行列
A_ls = np.column_stack([x_data, np.ones(n_data)])
# 正規方程式で最小二乗解
ATA = A_ls.T @ A_ls
ATb = A_ls.T @ y_data
x_hat = np.linalg.solve(ATA, ATb)
# 射影行列
P_ls = A_ls @ np.linalg.inv(ATA) @ A_ls.T
y_hat = P_ls @ y_data
residuals = y_data - y_hat
print(f"最小二乗推定: a = {x_hat[0]:.4f}, b = {x_hat[1]:.4f}")
print(f"真のパラメータ: a = 2.0, b = 1.0")
# (a) データとフィッティング
ax = axes[0]
ax.scatter(x_data, y_data, c='steelblue', s=30, zorder=5, label='Data')
ax.plot(x_data, x_hat[0] * x_data + x_hat[1], 'r-', linewidth=2,
label=f'Fit: $y={x_hat[0]:.2f}x + {x_hat[1]:.2f}$')
ax.plot(x_data, y_true, 'g--', linewidth=1.5, alpha=0.7, label='True: $y=2x+1$')
# 残差の可視化
for i in range(0, n_data, 3):
ax.plot([x_data[i], x_data[i]], [y_data[i], y_hat[i]], 'k-', alpha=0.3, linewidth=0.8)
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$y$', fontsize=12)
ax.set_title('Least Squares Fitting', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (b) 残差の分析
ax = axes[1]
ax.scatter(y_hat, residuals, c='steelblue', s=30)
ax.axhline(y=0, color='red', linewidth=1.5, linestyle='--')
ax.set_xlabel('Fitted values $\\hat{y}$', fontsize=12)
ax.set_ylabel('Residuals $e$', fontsize=12)
ax.set_title('Residual Plot', fontsize=13)
ax.grid(True, alpha=0.3)
# A^T e = 0 の確認
print(f"\nA^T e = {A_ls.T @ residuals}") # ゼロに近いはず
# (c) レバレッジ（ハット行列の対角成分）
ax = axes[2]
leverage = np.diag(P_ls)
ax.bar(range(n_data), leverage, color='steelblue', alpha=0.7)
ax.axhline(y=2/n_data, color='red', linewidth=1.5, linestyle='--',
label=f'Average: $p/n = {2/n_data:.3f}$')
ax.axhline(y=2*2/n_data, color='orange', linewidth=1.5, linestyle='--',
label=f'Threshold: $2p/n = {4/n_data:.3f}$')
ax.set_xlabel('Data point index', fontsize=12)
ax.set_ylabel('Leverage $h_{ii}$', fontsize=12)
ax.set_title('Leverage Values (Hat Matrix Diagonal)', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
print(f"\ntr(P) = {np.trace(P_ls):.4f} (should be {A_ls.shape[1]})")
print(f"SSE = {np.sum(residuals**2):.4f}")
print(f"SSR = {np.sum(y_hat**2):.4f}")
plt.tight_layout()
plt.savefig('least_squares.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、最小二乗法の重要な性質が読み取れます。
-
左図（フィッティング）: 推定された回帰直線（赤）が真の直線（緑破線）に近く、データ点の「中央」を通っています。各データ点から回帰直線への垂直距離（黒線）の二乗和が最小化されています
-
中央図（残差プロット）: 残差がゼロの周りにランダムに散らばっており、系統的なパターンが見られません。これは線形モデルの仮定が適切であることを示しています。もし残差に曲線的なパターンが見られれば、非線形な関係が存在する示唆になります
-
右図（レバレッジ）: 端のデータ点（$x$ が小さい点や大きい点）のレバレッジが大きくなっています。これは直感に合います — 端のデータ点は回帰直線の傾きに大きな影響を与えるからです。赤線（平均値 $p/n$）とオレンジ線（閾値 $2p/n$）により、影響の大きいデータ点を特定できます
射影行列の固有値分布
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# 射影行列の固有値
ax = axes[0]
eigenvalues_P = np.linalg.eigvalsh(P_ls)
ax.bar(range(n_data), sorted(eigenvalues_P, reverse=True),
color='steelblue', alpha=0.7)
ax.set_xlabel('Index', fontsize=12)
ax.set_ylabel('Eigenvalue', fontsize=12)
ax.set_title('Eigenvalues of Projection Matrix $P$', fontsize=13)
ax.axhline(y=0.5, color='red', linewidth=1, linestyle='--', alpha=0.5)
ax.grid(True, alpha=0.3)
# (I-P) の固有値
ax = axes[1]
eigenvalues_IP = np.linalg.eigvalsh(np.eye(n_data) - P_ls)
ax.bar(range(n_data), sorted(eigenvalues_IP, reverse=True),
color='coral', alpha=0.7)
ax.set_xlabel('Index', fontsize=12)
ax.set_ylabel('Eigenvalue', fontsize=12)
ax.set_title('Eigenvalues of $I - P$ (Residual Maker)', fontsize=13)
ax.axhline(y=0.5, color='red', linewidth=1, linestyle='--', alpha=0.5)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('projection_eigenvalues.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、射影行列の固有値の性質が確認できます。
-
左図（$\bm{P}$ の固有値）: 固有値が0と1のみに集中しています。固有値1が2個（= ランク2 = パラメータ数）、固有値0が28個（= $n – p$）です。冪等行列の固有値は0と1のみであるという理論が数値的にも確認されています
-
右図（$\bm{I} – \bm{P}$ の固有値）: 補射影行列も同様に0と1の固有値のみを持ちますが、$\bm{P}$ と正反対に固有値1が28個、固有値0が2個です。$\bm{P}$ と $\bm{I} – \bm{P}$ が直交補空間への射影であることが固有値の「反転」に反映されています
正規方程式の数値的安定性
正規方程式 $\bm{A}^T\bm{A}\hat{\bm{x}} = \bm{A}^T\bm{b}$ には数値的な注意点があります。$\bm{A}^T\bm{A}$ の条件数は $\bm{A}$ の条件数の二乗になるため
$$ \kappa(\bm{A}^T\bm{A}) = \kappa(\bm{A})^2 $$
$\bm{A}$ の条件数が $10^6$ 程度でも、$\bm{A}^T\bm{A}$ の条件数は $10^{12}$ になり、数値的に不安定になります。
より安定な方法として、QR分解を使う方法があります。$\bm{A} = \bm{QR}$ と分解すると
$$ \bm{R}\hat{\bm{x}} = \bm{Q}^T\bm{b} $$
$\bm{R}$ は上三角行列なので後退代入で解けます。条件数が二乗にならないため、数値的に安定です。
さらに安定な方法としてSVDを使う方法があり、$\bm{A}$ のランクが不足する場合でも擬似逆行列を通じて最小ノルム解を得ることができます。
重み付き最小二乗法への拡張
標準的な最小二乗法では、全てのデータ点を等しく扱います。しかし、計測精度がデータ点ごとに異なる場合、精度の高いデータ点をより重視したいことがあります。これを実現するのが重み付き最小二乗法（weighted least squares, WLS）です。
重み行列 $\bm{W} = \text{diag}(w_1, w_2, \ldots, w_m)$（$w_i > 0$）を導入し、重み付き残差の二乗和を最小化します。
$$ \hat{\bm{x}} = \arg\min_{\bm{x}} (\bm{b} – \bm{Ax})^T \bm{W} (\bm{b} – \bm{Ax}) $$
微分してゼロとおくと、重み付き正規方程式が得られます。
$$ \bm{A}^T\bm{W}\bm{A}\hat{\bm{x}} = \bm{A}^T\bm{W}\bm{b} $$
対応する射影行列は $\bm{P}_W = \bm{A}(\bm{A}^T\bm{W}\bm{A})^{-1}\bm{A}^T\bm{W}$ です。この射影行列は $\bm{W}$ が定める内積に関して直交射影を行います。つまり、重み $w_i$ が大きいデータ点ほど、回帰直線がそのデータ点に近づくように引き寄せられます。
重み付き最小二乗法は、計測誤差の分散 $\sigma_i^2$ が既知の場合に $w_i = 1/\sigma_i^2$ と設定する一般化最小二乗法（GLS）の特殊な場合です。GLS推定量は、ガウス・マルコフの定理により、不偏線形推定量のクラスで最小分散を達成します。
まとめ
本記事では、射影行列と最小二乗法の理論を幾何学的な視点から解説しました。
- 射影行列 $\bm{P} = \bm{A}(\bm{A}^T\bm{A})^{-1}\bm{A}^T$ はベクトルを列空間に直交射影する行列で、冪等性 $\bm{P}^2 = \bm{P}$ と対称性 $\bm{P}^T = \bm{P}$ を持つ
- 最小二乗法は $\bm{b}$ を列空間に射影する操作であり、正規方程式 $\bm{A}^T\bm{A}\hat{\bm{x}} = \bm{A}^T\bm{b}$ から解が求まる
- 残差は列空間に直交し（$\bm{A}^T\bm{e} = \bm{0}$）、二乗和は $\|\bm{b}\|^2 = \|\hat{\bm{b}}\|^2 + \|\bm{e}\|^2$ と分解される
- ハット行列の対角成分（レバレッジ）は各データ点の影響度を測る
- 正規方程式は条件数が二乗になるため、数値的にはQR分解やSVDを使う方が安定
次のステップとして、以下の記事も参考にしてください。
- QR分解とグラム・シュミット法 — 数値的に安定な最小二乗法
- ムーア・ペンローズ擬逆行列の理論 — ランク不足の場合の一般化
- 特異値分解(SVD)の導出と応用 — 最も安定な最小二乗法