連立方程式 $\bm{A}\bm{x} = \bm{b}$ を解きたいとき、$\bm{A}$ が正方かつ正則なら $\bm{x} = \bm{A}^{-1}\bm{b}$ で解が求まります。しかし現実の問題では、$\bm{A}$ が長方形だったり、正方でもランク落ちしていたりして、通常の逆行列が存在しないことが頻繁にあります。
データ点が方程式の数より多い過決定系（$m > n$）では厳密解が存在せず、パラメータが方程式の数より多い劣決定系（$m < n$）では解が無数に存在します。このような場合に「最も合理的な解」を与えてくれるのがムーア・ペンローズ擬逆行列（Moore-Penrose pseudoinverse）$\bm{A}^+$ です。
擬逆行列を理解すると、以下のような応用が開けます。
- 最小二乗法: 過決定系の残差最小解の統一的な表現
- 最小ノルム解: 劣決定系における最小ノルム解の選択
- ロボティクス: 冗長マニピュレータの逆運動学
- 制御工学: 劣決定な制御入力の配分問題
- 機械学習: 正則化されていない線形回帰の解析解
本記事の内容
- ムーア・ペンローズ条件（4つの条件）の定義と直感的理解
- SVDによる擬逆行列の構成と証明
- 存在と一意性の証明
- 最小二乗解と最小ノルム解としての特徴付け
- Pythonによる実装と応用例
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 特異値分解（SVD） — 擬逆行列の構成に不可欠
- 行列の基本演算 — ランク、核空間、値域の概念
- QR分解とグラム・シュミット法 — 最小二乗法との関連
ムーア・ペンローズ擬逆行列とは — 直感的な理解
擬逆行列の直感を掴むために、まず「逆行列ができること」を整理しましょう。
正則行列 $\bm{A}$ の逆行列 $\bm{A}^{-1}$ は、$\bm{A}$ の操作を「完全に取り消す」ものです。$\bm{A}$ が $\bm{x}$ を $\bm{b}$ に写すなら、$\bm{A}^{-1}$ は $\bm{b}$ を $\bm{x}$ に戻します。
しかし、$\bm{A}$ が長方形のとき、この「完全な取り消し」は不可能です。たとえば $3 \times 2$ 行列 $\bm{A}$ は3次元空間の中の2次元平面にしか写せないので、平面の外にある $\bm{b}$ を正確に戻すことはできません。
擬逆行列 $\bm{A}^+$ は「できる限り逆行列に近いもの」です。
- $\bm{b}$ が $\bm{A}$ の値域にあれば、$\bm{A}^+\bm{b}$ は $\bm{A}\bm{x} = \bm{b}$ の解を返す
- $\bm{b}$ が値域にないなら、$\bm{A}^+\bm{b}$ は残差 $\|\bm{A}\bm{x} – \bm{b}\|$ を最小にする $\bm{x}$ を返す
- 解が複数あるなら、$\bm{A}^+\bm{b}$ はその中でノルム $\|\bm{x}\|$ が最小のものを返す
日常的なアナロジーで言えば、擬逆行列は「ベストエフォートの逆操作」です。完全に元に戻せなくても、最も近い状態まで戻してくれます。
歴史的には、E. H. ムーア（1920年）が「一般逆行列」の概念を最初に提案し、ロジャー・ペンローズ（1955年）が4つの条件による公理的な定義を確立しました。ペンローズは後にノーベル物理学賞を受賞する物理学者ですが、この数学的な貢献も深遠なものです。
直感を掴んだところで、ペンローズによる公理的な定義を見ていきましょう。
ムーア・ペンローズ条件
4つの条件
$m \times n$ 行列 $\bm{A}$ に対して、$n \times m$ 行列 $\bm{A}^+$ がムーア・ペンローズ擬逆行列であるとは、以下の4条件（ペンローズ条件）を全て満たすことです。
$$ \begin{align} (1) \quad & \bm{A}\bm{A}^+\bm{A} = \bm{A} \\ (2) \quad & \bm{A}^+\bm{A}\bm{A}^+ = \bm{A}^+ \\ (3) \quad & (\bm{A}\bm{A}^+)^T = \bm{A}\bm{A}^+ \\ (4) \quad & (\bm{A}^+\bm{A})^T = \bm{A}^+\bm{A} \end{align} $$
各条件の意味
条件(1): $\bm{A}\bm{A}^+\bm{A} = \bm{A}$ — これは $\bm{A}^+$ が $\bm{A}$ の「弱い逆行列」であることを意味します。$\bm{A}$ で写して $\bm{A}^+$ で戻して再び $\bm{A}$ で写すと、最初の $\bm{A}$ と同じ結果になります。$\bm{A}$ の値域上では $\bm{A}^+$ が正しく逆操作を行っていることを保証します。
条件(2): $\bm{A}^+\bm{A}\bm{A}^+ = \bm{A}^+$ — 条件(1)の「逆方向」です。$\bm{A}^+$ で写して $\bm{A}$ で写して再び $\bm{A}^+$ で写すと、最初の $\bm{A}^+$ と同じになります。
条件(3): $(\bm{A}\bm{A}^+)^T = \bm{A}\bm{A}^+$ — $\bm{A}\bm{A}^+$ は対称行列です。幾何学的には、$\bm{A}\bm{A}^+$ は $\bm{A}$ の値域（列空間）への直交射影であることを意味します。
条件(4): $(\bm{A}^+\bm{A})^T = \bm{A}^+\bm{A}$ — $\bm{A}^+\bm{A}$ も対称行列です。こちらは $\bm{A}$ の行空間への直交射影です。
条件(3)と(4)が「直交射影」を要求している点が重要です。一般の「一般逆行列」（条件(1)のみ満たすもの）は無数に存在しますが、直交射影の条件を追加することで擬逆行列は一意に定まります。
これらの条件が結果的に「最小二乗・最小ノルム解」を与えることは、後で証明します。次にSVDを使った擬逆行列の具体的な構成方法を見ましょう。
SVDによる擬逆行列の構成
特異値分解の復習
$m \times n$ 行列 $\bm{A}$ のランクを $r$ とします。特異値分解（SVD）は
$$ \bm{A} = \bm{U}\bm{\Sigma}\bm{V}^T $$
ここで $\bm{U}$（$m \times m$）と $\bm{V}$（$n \times n$）は直交行列、$\bm{\Sigma}$（$m \times n$）は対角成分に特異値 $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r > 0$ を持つ行列です。
擬逆行列の定義（SVD経由）
$\bm{A} = \bm{U}\bm{\Sigma}\bm{V}^T$ のとき、擬逆行列は
$$ \begin{equation} \bm{A}^+ = \bm{V}\bm{\Sigma}^+\bm{U}^T \end{equation} $$
ここで $\bm{\Sigma}^+$（$n \times m$）は $\bm{\Sigma}$ の非ゼロ対角成分の逆数を取り、転置したものです。
$$ \bm{\Sigma} = \begin{pmatrix} \sigma_1 & & & 0 \\ & \ddots & & \vdots \\ & & \sigma_r & 0 \\ 0 & \cdots & 0 & 0 \end{pmatrix}_{m \times n} \quad \Rightarrow \quad \bm{\Sigma}^+ = \begin{pmatrix} 1/\sigma_1 & & & 0 \\ & \ddots & & \vdots \\ & & 1/\sigma_r & 0 \\ 0 & \cdots & 0 & 0 \end{pmatrix}_{n \times m} $$
直感的には、SVDで $\bm{A}$ を「拡大・縮小」の各方向に分解し、各方向の拡大率の逆数を取ることで「逆操作」を作っています。拡大率がゼロの方向（$\sigma_k = 0$）には逆操作できないので、ゼロのまま残します。
ペンローズ条件の検証
$\bm{A}^+ = \bm{V}\bm{\Sigma}^+\bm{U}^T$ がペンローズ4条件を満たすことを確認します。
条件(1): $\bm{A}\bm{A}^+\bm{A} = \bm{U}\bm{\Sigma}\bm{V}^T \cdot \bm{V}\bm{\Sigma}^+\bm{U}^T \cdot \bm{U}\bm{\Sigma}\bm{V}^T$
直交行列の性質 $\bm{V}^T\bm{V} = \bm{I}$ と $\bm{U}^T\bm{U} = \bm{I}$ を使うと
$$ = \bm{U}\bm{\Sigma}\bm{\Sigma}^+\bm{\Sigma}\bm{V}^T $$
$\bm{\Sigma}\bm{\Sigma}^+\bm{\Sigma}$ を計算すると、各対角成分は $\sigma_i \cdot (1/\sigma_i) \cdot \sigma_i = \sigma_i$（$i \leq r$）またはゼロ（$i > r$）なので $\bm{\Sigma}\bm{\Sigma}^+\bm{\Sigma} = \bm{\Sigma}$ です。
よって $\bm{A}\bm{A}^+\bm{A} = \bm{U}\bm{\Sigma}\bm{V}^T = \bm{A}$ で条件(1)が成立します。
条件(2)も同様に $\bm{A}^+\bm{A}\bm{A}^+ = \bm{V}\bm{\Sigma}^+\bm{\Sigma}\bm{\Sigma}^+\bm{U}^T = \bm{V}\bm{\Sigma}^+\bm{U}^T = \bm{A}^+$ で成立します。
条件(3): $\bm{A}\bm{A}^+ = \bm{U}\bm{\Sigma}\bm{\Sigma}^+\bm{U}^T$
$\bm{\Sigma}\bm{\Sigma}^+$ は対角行列（対角成分が1または0）なので対称です。$\bm{U}\bm{D}\bm{U}^T$ の形は明らかに対称行列です。
条件(4): $\bm{A}^+\bm{A} = \bm{V}\bm{\Sigma}^+\bm{\Sigma}\bm{V}^T$
同様に対称です。$\square$
一意性の証明
ペンローズ4条件を満たす行列は一意であることを示します。
$\bm{X}$ と $\bm{Y}$ が共にペンローズ4条件を満たすとします。
$\bm{A}\bm{X}$ は対称射影行列なので $\bm{A}\bm{X} = (\bm{A}\bm{X})^T$、同様に $\bm{A}\bm{Y} = (\bm{A}\bm{Y})^T$。
条件(1)より $\bm{A}\bm{X}\bm{A} = \bm{A}$ なので $\bm{A}\bm{X}\bm{A}\bm{Y} = \bm{A}\bm{Y}$。両辺の転置を取ると $\bm{Y}^T\bm{A}^T(\bm{A}\bm{X})^T = (\bm{A}\bm{Y})^T$、すなわち $\bm{Y}^T\bm{A}^T\bm{A}\bm{X} = \bm{A}\bm{Y}$ です（条件(3)を使用）。
同様に $\bm{Y}\bm{A}\bm{X} = \bm{Y}\bm{A}\bm{X}\bm{A}\bm{X} = \bm{Y}\bm{A}\bm{X}$ を経由して（条件(4)を利用）、$\bm{X} = \bm{Y}$ が示されます。 $\square$
SVDによる構成と一意性が確認できたので、次に擬逆行列が最小二乗法と最小ノルム解にどう結びつくかを見ましょう。
最小二乗解と最小ノルム解
過決定系: 最小二乗解
$m > n$ で $\text{rank}(\bm{A}) = n$ の場合、$\bm{A}\bm{x} = \bm{b}$ は一般に厳密解を持ちません。この場合の最小二乗解は
$$ \hat{\bm{x}} = \arg\min_{\bm{x}} \|\bm{A}\bm{x} – \bm{b}\|^2 $$
であり、$\hat{\bm{x}} = \bm{A}^+\bm{b} = (\bm{A}^T\bm{A})^{-1}\bm{A}^T\bm{b}$ で与えられます。
証明: 残差 $\bm{r} = \bm{b} – \bm{A}\bm{x}$ を $\bm{A}$ の列空間（値域）$\mathcal{R}(\bm{A})$ とその直交補空間に分解します。
$\bm{A}\bm{A}^+$ は $\mathcal{R}(\bm{A})$ への直交射影（条件(3)より）なので
$$ \bm{b} = \bm{A}\bm{A}^+\bm{b} + (\bm{I} – \bm{A}\bm{A}^+)\bm{b} $$
$\hat{\bm{x}} = \bm{A}^+\bm{b}$ とおくと、残差は
$$ \bm{r} = \bm{b} – \bm{A}\hat{\bm{x}} = \bm{b} – \bm{A}\bm{A}^+\bm{b} = (\bm{I} – \bm{A}\bm{A}^+)\bm{b} $$
これは $\mathcal{R}(\bm{A})$ に直交する成分のみを含むので、$\|\bm{r}\|$ はこれ以上小さくできません。$\square$
劣決定系: 最小ノルム解
$m < n$ で $\text{rank}(\bm{A}) = m$ の場合、$\bm{A}\bm{x} = \bm{b}$ は無数の解を持ちます。その中でノルムが最小の解は
$$ \hat{\bm{x}} = \arg\min_{\bm{x}} \|\bm{x}\|^2 \quad \text{subject to} \quad \bm{A}\bm{x} = \bm{b} $$
であり、$\hat{\bm{x}} = \bm{A}^+\bm{b} = \bm{A}^T(\bm{A}\bm{A}^T)^{-1}\bm{b}$ で与えられます。
証明: $\bm{A}\bm{x} = \bm{b}$ の一般解は $\bm{x} = \bm{A}^+\bm{b} + (\bm{I} – \bm{A}^+\bm{A})\bm{z}$（$\bm{z}$ は任意のベクトル）です。ここで $\bm{A}^+\bm{A}$ は $\bm{A}$ の行空間への直交射影（条件(4)より）なので、$\bm{A}^+\bm{b}$ は行空間に属し、$(\bm{I} – \bm{A}^+\bm{A})\bm{z}$ は核空間に属します。
行空間と核空間は直交するので、ピタゴラスの定理より
$$ \|\bm{x}\|^2 = \|\bm{A}^+\bm{b}\|^2 + \|(\bm{I} – \bm{A}^+\bm{A})\bm{z}\|^2 \geq \|\bm{A}^+\bm{b}\|^2 $$
等号は $(\bm{I} – \bm{A}^+\bm{A})\bm{z} = \bm{0}$ のとき、すなわち $\bm{x} = \bm{A}^+\bm{b}$ のときに成立します。$\square$
一般の場合: 最小二乗最小ノルム解
ランク落ちの一般の場合は、$\bm{A}^+\bm{b}$ は
$$ \hat{\bm{x}} = \bm{A}^+\bm{b} = \arg\min_{\bm{x}} \|\bm{x}\| \quad \text{subject to} \quad \|\bm{A}\bm{x} – \bm{b}\| = \min $$
つまり、残差ノルムを最小にする解の中で、さらにベクトルのノルムが最小のものを選びます。これが擬逆行列の最も一般的な特徴付けです。
これらの性質を確認するために、Pythonで実装してみましょう。
Pythonでの実装
SVDによる擬逆行列の計算
まず、SVDを使って擬逆行列をスクラッチで計算し、NumPyの結果と比較します。
import numpy as np
def pseudoinverse_svd(A, tol=1e-10):
"""SVDによるムーア・ペンローズ擬逆行列の計算"""
U, s, Vt = np.linalg.svd(A, full_matrices=True)
# 特異値の逆数（ゼロに近いものは無視）
r = np.sum(s > tol)
S_plus = np.zeros((A.shape[1], A.shape[0]))
for i in range(r):
S_plus[i, i] = 1.0 / s[i]
return Vt.T @ S_plus @ U.T
# テスト1: フルランクの長方形行列（過決定系）
np.random.seed(42)
A1 = np.random.randn(5, 3)
A1_pinv_my = pseudoinverse_svd(A1)
A1_pinv_np = np.linalg.pinv(A1)
print("=== フルランク長方形行列 (5x3) ===")
print(f"自作 vs NumPy: ||差|| = {np.linalg.norm(A1_pinv_my - A1_pinv_np):.2e}")
# ペンローズ4条件の検証
print("\nペンローズ条件の検証:")
print(f" (1) ||AA+A - A|| = {np.linalg.norm(A1 @ A1_pinv_my @ A1 - A1):.2e}")
print(f" (2) ||A+AA+ - A+|| = {np.linalg.norm(A1_pinv_my @ A1 @ A1_pinv_my - A1_pinv_my):.2e}")
P1 = A1 @ A1_pinv_my
P2 = A1_pinv_my @ A1
print(f" (3) ||(AA+)^T - AA+|| = {np.linalg.norm(P1.T - P1):.2e}")
print(f" (4) ||(A+A)^T - A+A|| = {np.linalg.norm(P2.T - P2):.2e}")
# テスト2: ランク落ち行列
print("\n=== ランク落ち行列 (4x4, rank=2) ===")
B = np.random.randn(4, 2)
A2 = B @ B.T # ランク2の4x4行列
A2_pinv = pseudoinverse_svd(A2)
print(f"rank(A) = {np.linalg.matrix_rank(A2)}")
print(f"ペンローズ条件(1): ||AA+A - A|| = {np.linalg.norm(A2 @ A2_pinv @ A2 - A2):.2e}")
print(f"ペンローズ条件(3): ||(AA+)^T - AA+|| = {np.linalg.norm((A2 @ A2_pinv).T - A2 @ A2_pinv):.2e}")
上のコードでは、SVDによる擬逆行列の計算をスクラッチで実装し、NumPyの np.linalg.pinv と結果が一致することを確認しています。さらにペンローズ4条件が数値的に成立していることを検証しています。フルランクの長方形行列とランク落ちの正方行列の両方で、全条件の誤差が機械イプシロンのオーダーに収まっていることが重要です。
最小二乗解と最小ノルム解の可視化
擬逆行列が「最小二乗最小ノルム解」を返すことを視覚的に確認します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) 過決定系: 最小二乗解
ax = axes[0]
A_over = np.array([[1, 0], [0, 1], [1, 1]], dtype=float)
b_over = np.array([1, 2, 4], dtype=float)
# 擬逆行列による最小二乗解
x_ls = np.linalg.pinv(A_over) @ b_over
# 残差ノルムの等高線
x1_range = np.linspace(-1, 4, 200)
x2_range = np.linspace(-1, 5, 200)
X1, X2 = np.meshgrid(x1_range, x2_range)
residual = np.zeros_like(X1)
for i in range(len(x1_range)):
for j in range(len(x2_range)):
xvec = np.array([X1[j, i], X2[j, i]])
residual[j, i] = np.linalg.norm(A_over @ xvec - b_over)
cs = ax.contour(X1, X2, residual, levels=20, cmap="Blues", alpha=0.7)
ax.plot(x_ls[0], x_ls[1], "r*", markersize=15, label=f"$A^+b$ = ({x_ls[0]:.2f}, {x_ls[1]:.2f})")
ax.set_xlabel("$x_1$", fontsize=12)
ax.set_ylabel("$x_2$", fontsize=12)
ax.set_title("Overdetermined: Least Squares", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.colorbar(cs, ax=ax, label="$\\|Ax - b\\|$")
# (b) 劣決定系: 最小ノルム解
ax = axes[1]
A_under = np.array([[1, 2, 1]], dtype=float)
b_under = np.array([3.0])
# 擬逆行列による最小ノルム解
x_mn = np.linalg.pinv(A_under) @ b_under
# 解の集合（Ax = b を満たす x の2Dスライス）
t1 = np.linspace(-2, 4, 100)
t2 = np.linspace(-2, 4, 100)
T1, T2 = np.meshgrid(t1, t2)
# x = (t1, t2, 3 - t1 - 2*t2)
X3 = 3 - T1 - 2 * T2
norms = np.sqrt(T1**2 + T2**2 + X3**2)
cs = ax.contourf(T1, T2, norms, levels=20, cmap="Reds", alpha=0.7)
# A^+b の最初の2成分
ax.plot(x_mn[0], x_mn[1], "b*", markersize=15,
label=f"$A^+b$ = ({x_mn[0]:.2f}, {x_mn[1]:.2f}, {x_mn[2]:.2f})")
ax.set_xlabel("$x_1$", fontsize=12)
ax.set_ylabel("$x_2$", fontsize=12)
ax.set_title("Underdetermined: Minimum Norm", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.colorbar(cs, ax=ax, label="$\\|x\\|$")
# (c) ランク落ち: 最小二乗最小ノルム解
ax = axes[2]
n_examples = 20
ranks = range(1, 6)
errors_ls = []
norms_x = []
for r in ranks:
errs = []
nrms = []
for _ in range(n_examples):
m, n = 10, 8
U = np.random.randn(m, r)
V = np.random.randn(r, n)
A = U @ V # ランクrの行列
b = np.random.randn(m)
x = np.linalg.pinv(A) @ b
errs.append(np.linalg.norm(A @ x - b))
nrms.append(np.linalg.norm(x))
errors_ls.append(np.mean(errs))
norms_x.append(np.mean(nrms))
ax2 = ax.twinx()
bars1 = ax.bar(np.array(list(ranks)) - 0.2, errors_ls, 0.35,
color="salmon", alpha=0.8, label="$\\|Ax - b\\|$ (residual)")
bars2 = ax2.bar(np.array(list(ranks)) + 0.2, norms_x, 0.35,
color="steelblue", alpha=0.8, label="$\\|x\\|$ (solution norm)")
ax.set_xlabel("rank(A)", fontsize=12)
ax.set_ylabel("Residual norm", fontsize=12, color="red")
ax2.set_ylabel("Solution norm", fontsize=12, color="blue")
ax.set_title("Rank-Deficient: LS + Min Norm", fontsize=13)
ax.set_xticks(list(ranks))
lines1, labels1 = ax.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax.legend(lines1 + lines2, labels1 + labels2, fontsize=9, loc="upper right")
ax.grid(True, alpha=0.3, axis="y")
plt.tight_layout()
plt.savefig("pseudoinverse_properties.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、擬逆行列の3つの性質が視覚的に確認できます。
-
左図（過決定系の最小二乗解）: 残差ノルム $\|\bm{A}\bm{x} – \bm{b}\|$ の等高線上で、赤い星が最小点に位置しています。これは $\bm{A}^+\bm{b}$ が残差を最小にする最小二乗解であることを示しています。
-
中央図（劣決定系の最小ノルム解）: 解の集合（$\bm{A}\bm{x} = \bm{b}$ を満たす $\bm{x}$ の空間）上でのノルム $\|\bm{x}\|$ が色で表示されています。青い星がノルム最小の点に位置しており、$\bm{A}^+\bm{b}$ が最小ノルム解を返していることが確認できます。
-
右図（ランク落ちの場合）: ランクが上がるにつれて残差ノルム（赤）が減少し、解のノルム（青）も変化しています。ランクが低いほど情報が失われるため残差が大きくなりますが、擬逆行列は常に最小二乗最小ノルム解を返しています。
応用: 劣決定系の逆運動学
ロボティクスでの応用例として、冗長マニピュレータの逆運動学を擬逆行列で解きます。
import numpy as np
import matplotlib.pyplot as plt
def forward_kinematics(theta, lengths):
"""順運動学: 関節角度から先端位置を計算"""
x, y = 0.0, 0.0
positions = [(x, y)]
angle_sum = 0.0
for t, l in zip(theta, lengths):
angle_sum += t
x += l * np.cos(angle_sum)
y += l * np.sin(angle_sum)
positions.append((x, y))
return np.array([x, y]), positions
def jacobian(theta, lengths):
"""ヤコビアン行列の計算"""
n = len(theta)
J = np.zeros((2, n))
for j in range(n):
angle_sum = np.sum(theta[:j+1])
for k in range(j, n):
angle_sum_k = np.sum(theta[:k+1])
# dpos/dtheta_j
angle_j = np.sum(theta[:j+1])
remaining_x = sum(l * (-np.sin(np.sum(theta[:k+1])))
for k, l in enumerate(lengths) if k >= j)
remaining_y = sum(l * np.cos(np.sum(theta[:k+1]))
for k, l in enumerate(lengths) if k >= j)
J[0, j] = remaining_x
J[1, j] = remaining_y
return J
# 3リンクマニピュレータ（2D位置制御 → 劣決定系）
n_joints = 3
lengths = [1.0, 0.8, 0.6]
theta = np.array([0.3, 0.5, -0.2])
# 目標軌道（円）
n_steps = 100
t = np.linspace(0, 2 * np.pi, n_steps)
target_path = np.column_stack([1.2 + 0.3 * np.cos(t), 0.5 + 0.3 * np.sin(t)])
# 擬逆行列による逆運動学
theta_history = [theta.copy()]
actual_path = []
dt = 0.1
for i in range(n_steps):
pos, _ = forward_kinematics(theta, lengths)
actual_path.append(pos)
# 目標速度
target = target_path[i]
dx = (target - pos) * 3.0 # 比例制御
# ヤコビアン
J = jacobian(theta, lengths)
# 擬逆行列による関節速度
dtheta = np.linalg.pinv(J) @ dx
theta = theta + dtheta * dt
theta_history.append(theta.copy())
actual_path = np.array(actual_path)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# (a) 軌道追従
ax = axes[0]
ax.plot(target_path[:, 0], target_path[:, 1], "r--", linewidth=2, label="Target")
ax.plot(actual_path[:, 0], actual_path[:, 1], "b-", linewidth=2, label="Actual")
# 最終姿勢を描画
_, positions = forward_kinematics(theta, lengths)
positions = np.array(positions)
ax.plot(positions[:, 0], positions[:, 1], "ko-", linewidth=3, markersize=8)
ax.plot(positions[0, 0], positions[0, 1], "ks", markersize=12, label="Base")
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_title("3-Link Manipulator: Pseudoinverse IK", fontsize=13)
ax.legend(fontsize=10)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
# (b) 関節角度の推移
ax = axes[1]
theta_hist = np.array(theta_history)
for j in range(n_joints):
ax.plot(theta_hist[:, j], linewidth=2, label=f"Joint {j+1}")
ax.set_xlabel("Step", fontsize=12)
ax.set_ylabel("Joint angle (rad)", fontsize=12)
ax.set_title("Joint Angle History", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("pseudoinverse_robotics.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、擬逆行列を用いた逆運動学の動作が確認できます。
-
左図（軌道追従）: 3リンクマニピュレータ（関節3つ、制御変数2つの劣決定系）が、擬逆行列による逆運動学で目標の円軌道（赤い破線）を追従しています。実際の軌道（青い実線）が目標に近い位置を辿っていることがわかります。擬逆行列は $\bm{J}^+$ によって最小ノルムの関節速度 $\Delta\bm{\theta}$ を選ぶため、不必要に大きな関節運動を避けています。
-
右図（関節角度の推移）: 3つの関節角度が滑らかに変化しています。擬逆行列が最小ノルム解を選ぶため、関節角度の変化が最小限に抑えられ、滑らかな運動が実現されています。
擬逆行列の計算コストと注意点
計算コスト
擬逆行列の計算コストは本質的にSVDのコストに等しく、$m \times n$ 行列に対して $O(\min(mn^2, m^2n))$ です。
実用上は以下の選択指針があります。
- フルランクの過決定系 ($m > n$, $\text{rank} = n$): $\bm{A}^+ = (\bm{A}^T\bm{A})^{-1}\bm{A}^T$ が効率的。QR分解経由でさらに安定
- フルランクの劣決定系 ($m < n$, $\text{rank} = m$): $\bm{A}^+ = \bm{A}^T(\bm{A}\bm{A}^T)^{-1}$ が効率的
- ランク落ち: SVD経由の $\bm{A}^+ = \bm{V}\bm{\Sigma}^+\bm{U}^T$ が必要
数値的な注意点
実際の計算では、ゼロに近い特異値の扱いが重要です。理論上はゼロでなくても、丸め誤差やノイズの影響で「ほぼゼロ」の特異値が現れることがあります。
- 閾値の設定: 小さな特異値を無視するための閾値（tolerance）を適切に設定する必要があります。NumPyの
pinvではデフォルトで $\text{tol} = \max(m, n) \cdot \varepsilon_{\text{mach}} \cdot \sigma_1$ を使用しています - 正則化: ティホノフ正則化 $\bm{A}^+_\alpha = (\bm{A}^T\bm{A} + \alpha\bm{I})^{-1}\bm{A}^T$ で小さな特異値の影響を抑制できます。この正則化はリッジ回帰と本質的に同じであり、$\alpha \to 0$ の極限で擬逆行列に一致します
まとめ
本記事では、ムーア・ペンローズ擬逆行列の理論について解説しました。
- ペンローズ4条件は擬逆行列を公理的に定義し、その一意性を保証する
- SVDによる構成 $\bm{A}^+ = \bm{V}\bm{\Sigma}^+\bm{U}^T$ はあらゆる行列に適用でき、特異値の逆数を取るという直感的な操作に基づく
- 過決定系では最小二乗解を、劣決定系では最小ノルム解を、ランク落ちの場合は最小二乗最小ノルム解を与える
- $\bm{A}\bm{A}^+$ は列空間への直交射影、$\bm{A}^+\bm{A}$ は行空間への直交射影である
- ロボティクスの逆運動学など、劣決定系の問題で実用的に重要な役割を果たす
次のステップとして、以下の記事も参考にしてください。
- QR分解とグラム・シュミット法 — 最小二乗法の安定な計算
- 行列の微分公式を完全まとめ — 最適化への応用
- 勾配・発散・回転の定義と直感的理解 — ベクトル解析の基本