NumPyによる特異値分解（SVD）の実装方法：理論から応用まで徹底解説
1. はじめに：特異値分解（SVD）とは何か、なぜ重要なのか
データサイエンス、機械学習、信号処理、画像処理といった様々な分野において、行列はデータの構造や関係性を表現するための基本的なツールです。行列が大きくなると、その性質を理解したり、効率的に処理したりすることが難しくなります。ここで強力なツールとなるのが「行列分解（Matrix Decomposition）」あるいは「行列因子分解（Matrix Factorization）」と呼ばれる手法群です。これは、与えられた行列を、より扱いやすく、構造が明確な複数の行列の積として表現する技術です。
行列分解の中でも特に強力で汎用的なのが特異値分解（Singular Value Decomposition, SVD）です。SVDは、正方行列だけでなく、任意の次元を持つ行列に対しても適用できるという特徴があります。これにより、どんな線形変換も、回転、スケーリング、そして再び回転という三つの基本的な操作の組み合わせとして理解できることが示されます。
SVDは多岐にわたる応用を持っています。例えば、データの次元削減（Principal Component Analysis; PCAとの関連）、ノイズ除去、画像圧縮、情報検索における潜在的意味解析（Latent Semantic Analysis; LSA）、レコメンデーションシステム（協調フィルタリング）、そして線形方程式の解法における疑似逆行列の計算などです。これらの応用は、SVDが行列の最も重要な情報を抽出し、冗長性を取り除く能力に基づいています。
Pythonにおいて、数値計算ライブラリのデファクトスタンダードであるNumPyは、SVDを効率的に計算するための関数を提供しています。本記事では、NumPyを用いてSVDを実装する方法を、その理論的な背景から具体的なコード例、そして様々な応用例まで、詳細かつ包括的に解説します。約5000語をかけて、SVDの理解を深め、実際にPythonで使いこなせるようになることを目指します。
2. 特異値分解（SVD）の理論：数学的な基礎
SVDの実装に入る前に、まずはその数学的な定義と性質を理解しておくことが重要です。SVDは、任意の $m \times n$ の実数行列 $A$ を、以下の形式に分解するものです。
$A = U \Sigma V^T$
ここで、各行列は以下の性質を持ちます。
- $U$：$m \times m$ のユニタリ行列（実数の場合は直交行列）。つまり、$U^T U = U U^T = I_m$ となります。$U$ の列ベクトルは、行列 $A$ の左特異ベクトルと呼ばれます。これらは $m$ 次元空間の正規直交基底を形成します。
- $\Sigma$：$m \times n$ の対角行列（厳密には「対角成分以外がゼロ」）。対角成分には、行列 $A$ の特異値 $\sigma_i$ が非負の実数として降順に並べられます。
$\Sigma = \begin{pmatrix}
\sigma_1 & 0 & \cdots & 0 & \cdots & 0 \
0 & \sigma_2 & \cdots & 0 & \cdots & 0 \
\vdots & \vdots & \ddots & \vdots & \ddots & \vdots \
0 & 0 & \cdots & \sigma_r & \cdots & 0 \
\vdots & \vdots & \ddots & \vdots & \ddots & \vdots \
0 & 0 & \cdots & 0 & \cdots & 0
\end{pmatrix}$
ここで、$r = \min(m, n)$ です。特異値 $\sigma_i$ は、一般的に $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_r \ge 0$ の順に並べられます。非ゼロの特異値の数は、行列 $A$ のランクに等しくなります。 - $V^T$：$n \times n$ のユニタリ行列 $V$ の転置行列（実数の場合は直交行列の転置）。つまり、$V^T V = V V^T = I_n$ となります。$V$ の列ベクトル（$V^T$ の行ベクトル）は、行列 $A$ の右特異ベクトルと呼ばれます。これらは $n$ 次元空間の正規直交基底を形成します。
重要な点として、SVDは任意の実数行列（あるいは複素数行列）に対して存在し、その分解は（特異値の順序と、複数の特異値が等しい場合の特異ベクトルの符号・選択を除いて）一意的です。
2.1. 特異値、特異ベクトルの意味
特異値 $\sigma_i$ は、行列 $A$ がベクトルを「どれだけ引き伸ばすか」の度合いを示します。具体的には、右特異ベクトル $v_i$（$V$ の $i$ 番目の列ベクトル）は、線形変換 $x \mapsto Ax$ によって $u_i$（$U$ の $i$ 番目の列ベクトル）の方向に写像され、その際に長さが $\sigma_i$ 倍に引き伸ばされるという関係があります。すなわち、$A v_i = \sigma_i u_i$ が成り立ちます。これは、SVDの定義式 $A = U \Sigma V^T$ を変形することで確認できます。右から $V$ をかけると $AV = U \Sigma$ となり、これを列ベクトルごとに見ると $A [v_1 | v_2 | \cdots | v_n] = [u_1 | u_2 | \cdots | u_m] \Sigma$ となり、最終的に $A v_i = \sigma_i u_i$ が得られます（ただし、$i > \min(m, n)$ の場合は $\sigma_i = 0$）。
左特異ベクトル $u_i$ は出力空間の基底を、右特異ベクトル $v_i$ は入力空間の基底を形成します。SVDは、これらの特別な正規直交基底の下で、線形変換 $A$ が単なる各軸方向へのスケーリングとして表現されることを示しています。変換 $A$ は、まず入力空間の基底 $v_1, \ldots, v_n$ に対して回転（$V^T$ の操作）を行い、次に各軸方向へのスケーリング（$\Sigma$ の操作）を行い、最後に新しい基底 $u_1, \ldots, u_m$ の方向へ回転（$U$ の操作）を行う、と解釈できます。
2.2. SVDと固有値分解の違い
行列分解のもう一つの重要な手法に固有値分解（Eigenvalue Decomposition）があります。固有値分解は、正方行列 $A$ を $A = P \Lambda P^{-1}$ の形に分解するものです。ここで、$P$ は固有ベクトルを並べた行列、$\Lambda$ は固有値を対角成分に並べた対角行列です。
SVDと固有値分解は関連がありますが、重要な違いがあります。
1. 適用範囲: 固有値分解は正方行列にのみ適用可能ですが、SVDは任意の $m \times n$ 行列に適用可能です。
2. 基底: 固有値分解における固有ベクトルは一般に直交しません（特に $A$ が対称行列でない場合）。一方、SVDの左特異ベクトル $U$ と右特異ベクトル $V$ は常に直交行列の列ベクトルです。
3. 存在: 固有値分解は、全ての正方行列に対して存在するわけではありません。特に、対角化可能でない行列は固有値分解できません。一方、SVDは全ての行列に対して存在します。
4. 値の性質: 固有値は複素数になることも、負になることもありますが、特異値は常に非負の実数です。
5. 関連性: 行列 $A$ のSVD $A = U \Sigma V^T$ を考えると、行列 $A^T A$ の固有値分解は $V (\Sigma^T \Sigma) V^T$ となります。ここで、$V$ は $A^T A$ の固有ベクトルからなる行列、$\Sigma^T \Sigma$ は対角成分に $A$ の特異値の2乗 $\sigma_i^2$ を並べた対角行列です。同様に、$A A^T$ の固有値分解は $U (\Sigma \Sigma^T) U^T$ となります。したがって、特異値の2乗は $A^T A$（または $A A^T$）の固有値であり、右特異ベクトル $V$ は $A^T A$ の固有ベクトル、左特異ベクトル $U$ は $A A^T$ の固有ベクトルに対応します。特に、$A$ が対称行列の場合は $A = A^T$ なので、$A^T A = A^2$ となり、SVDと固有値分解は密接に関連します。
SVDの最大の強みは、その普遍性と安定性にあります。どんな行列に対しても適用でき、数値的に安定して計算できます。
3. NumPyによるSVDの実装：numpy.linalg.svd を使う
PythonでSVDを計算するには、主にNumPyライブラリの linalg サブモジュールにある svd 関数を使用します。
“`python
import numpy as np
例として、3×4 の行列を作成
A = np.array([[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12]])
print(“元の行列 A:”)
print(A)
“`
numpy.linalg.svd 関数の基本的な使い方は非常にシンプルです。行列を引数として渡すだけです。
“`python
SVDを実行
U, s, Vh = np.linalg.svd(A)
print(“\n左特異ベクトル行列 U:”)
print(U)
print(“\n特異値配列 s:”)
print(s)
print(“\n右特異ベクトル行列 Vh (Vの転置):”)
print(Vh)
“`
3.1. numpy.linalg.svd の戻り値
numpy.linalg.svd(A) は、デフォルトでは以下の3つのオブジェクトを返します。
U: 行列 $A$ の左特異ベクトルからなる行列 $U$ です。これは $m \times m$ の直交行列です。NumPyの出力では、列ベクトルが左特異ベクトルに対応します。s: 特異値の配列です。これは対角行列 $\Sigma$ の対角成分 $\sigma_1, \sigma_2, \ldots, \sigma_{\min(m, n)}$ を降順に並べた1次元配列です。$\Sigma$ 自体が必要な場合は、この配列を使って別途作成する必要があります。Vh: 行列 $A$ の右特異ベクトルからなる行列 $V$ の転置 $V^T$（あるいは $V^H$、共役転置）です。これは $n \times n$ の直交行列です。NumPyの出力では、行ベクトルが右特異ベクトル $v_i^T$ に対応します。
上記の例では、行列 $A$ は $3 \times 4$ なので、$m=3, n=4$ です。
* U の形状は $(3, 3)$ となります。
* s の形状は $(3$,$) となります（\min(3, 4)=3$）。
* Vh の形状は $(4, 4)$ となります。
3.2. 取得した要素からの元の行列の復元
SVDの定義 $A = U \Sigma V^T$ に従って、NumPyから得られた U, s, Vh を使って元の行列 $A$ を復元できることを確認してみましょう。注意が必要なのは、s は対角行列 $\Sigma$ そのものではなく、対角成分の配列であるという点です。$\Sigma$ 行列を作成するには、np.diag() 関数を使います。また、$m \times n$ 行列 $A$ のSVDにおいて、$\Sigma$ は $m \times n$ 行列ですが、s の長さは $\min(m, n)$ です。NumPyの linalg.svd は、デフォルトでReduced SVDと呼ばれる形式を計算します。この場合、$U$ は $m \times k$、$s$ は長さ $k$、$V^T$ は $k \times n$ となります。ここで $k = \min(m, n)$ です。NumPyから得られる U は $m \times m$ ですが、s の長さは $\min(m, n)$ であり、Vh は $n \times n$ です。元の行列を復元するためには、適切な次元の $\Sigma$ 行列を作成する必要があります。
NumPyの svd 関数が返す U, s, Vh を使って元の行列を復元するには、s を対角行列に変換し、その形状を $m \times n$ に合わせる必要があります。Reduced SVDの形式に合わせて復元するのが一般的です。Reduced SVDでは、$A = U_k \Sigma_k V_k^T$ となります。ここで $U_k$ は $U$ の最初の $\min(m, n)$ 列、$V_k$ は $V$ の最初の $\min(m, n)$ 列、$\Sigma_k$ は $\Sigma$ の左上 $\min(m, n) \times \min(m, n)$ 部分です。NumPyの svd は、full_matrices=False の場合にこのReduced SVDの $U_k$ と $V_k^T$ を返します。しかし、full_matrices=True (デフォルト)の場合でも、返される s の長さは $\min(m, n)$ であり、この s を使って $\min(m, n) \times \min(m, n)$ の対角行列を作成し、$U$ の最初の $\min(m, n)$ 列と $V^T$ の最初の $\min(m, n)$ 行を使って復元するのが一般的です。NumPyの戻り値 U の最初の $\min(m, n)$ 列、s 全体、Vh の最初の $\min(m, n)$ 行を使って復元します。
“`python
特異値配列 s から対角行列 Sigma を作成
s は長さ min(m, n) の配列
U は m x m, Vh は n x n (デフォルト full_matrices=True の場合)
A は m x n なので、復元には U の最初の min(m, n) 列と Vh の最初の min(m, n) 行を使うか、
s から m x n の対角行列を作成する必要がある。
一般的には、Reduced SVD の形式で復元する。
Reduced SVD: A = U[:, :k] @ Sigma[:k, :k] @ Vh[:k, :]
k = min(m, n)
k = min(A.shape)
Sigma = np.zeros((A.shape[0], A.shape[1]))
対角部分に特異値をセット
Sigma[:k, :k] = np.diag(s)
元の行列 A を復元
A_restored = U @ Sigma @ Vh
print(“\n復元された行列 A_restored:”)
print(A_restored)
元の行列と復元された行列がほぼ同じであることを確認
print(“\n元の行列と復元された行列の差の最大絶対値:”)
print(np.max(np.abs(A – A_restored))) # 浮動小数点演算誤差程度の値になるはず
“`
出力を見ると、元の行列と復元された行列が非常に近い値になっていることが確認できます。差の最大絶対値が非常に小さい（通常は $10^{-14}$ のオーダー）のは、浮動小数点演算による誤差です。
3.3. full_matrices 引数
numpy.linalg.svd には full_matrices という重要な引数があります。デフォルトは True です。
-
full_matrices=True(デフォルト): Full SVD を計算します。Uは $m \times m$ の形状になります。Vhは $n \times n$ の形状になります。sは長さ $\min(m, n)$ の配列になります。
このモードでは、UとV^Tは完全な直交行列を提供しますが、特異値がゼロに対応する列/行も含まれるため、計算量やメモリ使用量が多くなる可能性があります。
-
full_matrices=False: Reduced SVD を計算します。Uは $m \times k$ の形状になります。ここで $k = \min(m, n)$ です。Vhは $k \times n$ の形状になります。sは長さ $k$ の配列になります。
このモードでは、UとV^Tは非ゼロ特異値に対応する部分のみを含みます。多くの応用（特に次元削減）では、非ゼロ特異値に対応する特異ベクトルのみが必要なため、このモードが効率的です。特異値が全て非ゼロの場合、Reduced SVDとFull SVDは同じ情報を含みますが、行列の形状が異なります。
Reduced SVD形式での復元は以下のようになります。full_matrices=False の場合の戻り値 U_reduced, s_reduced, Vh_reduced を使います。
“`python
Reduced SVD を実行
U_reduced, s_reduced, Vh_reduced = np.linalg.svd(A, full_matrices=False)
print(“\nReduced SVD – 左特異ベクトル行列 U_reduced:”)
print(U_reduced.shape)
print(U_reduced)
print(“\nReduced SVD – 特異値配列 s_reduced:”)
print(s_reduced.shape)
print(s_reduced)
print(“\nReduced SVD – 右特異ベクトル行列 Vh_reduced (Vの転置):”)
print(Vh_reduced.shape)
print(Vh_reduced)
Reduced SVD 結果から元の行列を復元
A = U_k @ Sigma_k @ Vh_k
ここで Sigma_k は s_reduced を対角成分とする k x k 対角行列
Sigma_k = np.diag(s_reduced)
A_restored_reduced = U_reduced @ Sigma_k @ Vh_reduced
print(“\nReduced SVD から復元された行列 A_restored_reduced:”)
print(A_restored_reduced)
元の行列と復元された行列がほぼ同じであることを確認
print(“\n元の行列とReduced SVD から復元された行列の差の最大絶対値:”)
print(np.max(np.abs(A – A_restored_reduced)))
``full_matrices=Falseの場合、U_reducedの形状は $(m, \min(m, n))$、Vh_reducedの形状は $(\min(m, n), n)$ となります。s_reducedは長さ $\min(m, n)$ の配列です。復元には、s_reducedを対角成分とする $\min(m, n) \times \min(m, n)$ の対角行列を作成し、U_reducedとVh_reduced` との積を計算します。これも元の行列とほぼ一致することが確認できます。
多くの場合、特に大規模な行列で非ゼロ特異値の数が少ない場合は、full_matrices=False を使う方がメモリ効率も計算効率も優れています。
3.4. compute_uv 引数
compute_uv 引数は、特異ベクトル行列 $U$ と $V^T$ を計算するかどうかを指定します。デフォルトは True です。
compute_uv=True(デフォルト): $U, s, V^T$ の全てを計算して返します。compute_uv=False: 特異値配列 $s$ のみ計算して返します。この場合、戻り値は1つの配列 (特異値配列) のみとなります。特異ベクトルが必要ない場合は、計算時間を節約できます。
“`python
特異値のみ計算
s_only = np.linalg.svd(A, compute_uv=False)
print(“\n特異値のみ (compute_uv=False):”)
print(s_only)
“`
これは、例えば行列のランクを調べたいだけの場合などに有用です。
4. SVDの応用例：NumPyを使った実装
SVDは非常に多岐にわたる応用を持っています。ここでは、NumPyのSVD機能を使ったいくつかの代表的な応用例を紹介します。
4.1. 低ランク近似と次元削減
SVDの最も一般的な応用の一つは、元の行列をより低いランクの行列で近似することです。これは、特異値が降順に並んでいることを利用します。大きい特異値に対応する特異ベクトルは、元の行列のデータにおける主要なパターンや情報を含んでいると考えられます。一方、小さい特異値はノイズや重要でない情報を表していることが多いです。
元の行列 $A$ のSVDが $A = U \Sigma V^T$ であるとき、上位 $k$ 個の最大の特異値 $\sigma_1, \ldots, \sigma_k$ と、それらに対応する左特異ベクトル $u_1, \ldots, u_k$ および右特異ベクトル $v_1, \ldots, v_k$ を使うことで、$A$ をランク $k$ の行列 $A_k$ で近似することができます。
$A_k = U_k \Sigma_k V_k^T$
ここで、$U_k$ は $U$ の最初の $k$ 列、$V_k^T$ は $V^T$ の最初の $k$ 行、$\Sigma_k$ は $\Sigma$ の左上 $k \times k$ 部分（対角成分が $\sigma_1, \ldots, \sigma_k$ の対角行列）です。この近似 $A_k$ は、フロベニウスノルムの意味で、ランク $k$ の全ての行列の中で元の行列 $A$ に最も近いことが知られています（Eckart-Young-Mirsky theorem）。
NumPyのSVD結果（Reduced SVD形式が都合が良い）を使って、この低ランク近似を実装してみましょう。
“`python
Reduced SVD を計算
U, s, Vh = np.linalg.svd(A, full_matrices=False)
元の行列 A: 3×4
U: 3×3 (full_matrices=True) or 3×3 (full_matrices=False) -> Reduced SVD では 3×3
s: 長さ 3
Vh: 4×4 (full_matrices=True) or 3×4 (full_matrices=False) -> Reduced SVD では 3×4
Reduced SVDの結果の形状を確認
print(“\nReduced SVD 結果の形状:”)
print(“U:”, U.shape)
print(“s:”, s.shape)
print(“Vh:”, Vh.shape)
上位 k 個の特異値/特異ベクトルを使って近似行列 A_k を作成
例えば、k=1 で近似してみる
k = 1
U_k: U の最初の k 列
U_k = U[:, :k]
s_k: s の最初の k 個の要素
s_k = s[:k]
Vh_k: Vh の最初の k 行
Vh_k = Vh[:k, :]
print(f”\nランク {k} 近似のための行列形状:”)
print(“U_k:”, U_k.shape)
print(“s_k:”, s_k.shape)
print(“Vh_k:”, Vh_k.shape)
近似行列 A_k を計算: A_k = U_k @ Sigma_k @ Vh_k
Sigma_k は s_k を対角成分とする k x k 行列
Sigma_k = np.diag(s_k)
A_k = U_k @ Sigma_k @ Vh_k
print(f”\n元の行列 A:\n{A}”)
print(f”\nランク {k} 近似行列 A_{k}:\n{A_k}”)
元の行列と近似行列の差
print(f”\n元の行列とランク {k} 近似行列の差:\n{A – A_k}”)
``s
この例では、元の行列 $A$ はランクが2であることが特異値を見ればわかります（が[30.8577…, 1.6999…, 0.]` のように、3番目の特異値がほぼゼロに近いことから）。したがって、ランクを1に削減すると情報がかなり失われ、近似行列は元の行列から大きく乖離します。しかし、ランクを2に設定すれば、ほぼ元の行列が復元されるはずです。
この低ランク近似は、特にデータ行列が大きい場合に、計算量やストレージを削減するのに役立ちます。また、小さい特異値に対応する成分を捨てることで、ノイズ除去の効果も期待できます。
次元削減という文脈では、データ行列 $A$（各行がデータ点、各列が特徴量）に対してSVDを適用し、上位 $k$ 個の特異ベクトルを用いてデータを新しい $k$ 次元空間に射影することを指す場合があります。特に、共分散行列に対してSVD（あるいは固有値分解）を行う主成分分析（PCA）は、これと密接に関連しています。PCAは、データ行列 $X$（各行がデータ点、各列が特徴量）の中心化（平均を引く）を行った後、$X^T X$（共分散行列に比例）の固有値分解を行います。前述のように、$A^T A$ の固有ベクトルは $A$ の右特異ベクトルです。したがって、中心化されたデータ行列 $X_{centered}$ に対してSVDを行うことでもPCAと同じ主成分を得ることができます。
データ行列 $X$ ($N$ 個のデータ点 $\times$ $D$ 個の特徴量, $N \times D$) を中心化した $X_{centered}$ に対してSVDを実行：
$X_{centered} = U \Sigma V^T$
ここで、$V$ ($D \times D$) の列ベクトルが主成分の方向を示します。データを $k$ 次元に削減するには、$V$ の最初の $k$ 列 $V_k$ ($D \times k$) を使用し、元のデータを行列積により新しい空間に射影します。
$X_{reduced} = X_{centered} V_k$
$X_{reduced}$ は $N \times k$ の行列となり、各行が $k$ 次元に削減されたデータ点となります。
これは、PCAを手軽に行う方法としてSVDがよく用いられる理由です。NumPyのSVDは、このようにデータ行列の次元削減にも直接的に応用できます。
4.2. 画像圧縮・ノイズ除去
画像のデータも、NumPyを使えば数値の行列として扱うことができます。グレースケール画像は1つの行列として、カラー画像はR, G, B各チャンネルを別々の行列として扱うことができます。SVDによる低ランク近似は、画像の主要な情報（大きな特異値に対応する成分）を保持しつつ、細かい情報やノイズ（小さな特異値に対応する成分）を削減するため、画像圧縮やノイズ除去に応用できます。
例として、グレースケール画像のSVDによる圧縮を実装してみましょう。画像ファイル（例えばPNG形式）を読み込むには、Pillow (PIL) やOpenCVといったライブラリが便利です。ここではPillowを使います。
まず、必要なライブラリをインストールします。
pip install Pillow
“`python
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt # 画像表示のため
画像を読み込み、グレースケールに変換
try:
# 適当な画像ファイル名を指定してください
# 例: ‘test_image.png’ が存在すると仮定
img_path = ‘test_image.png’ # <– 実際の画像ファイル名に置き換えてください
img = Image.open(img_path).convert(‘L’) # グレースケールに変換
img_matrix = np.array(img)
except FileNotFoundError:
print(f”Error: 画像ファイル ‘{img_path}’ が見つかりません。”)
print(“サンプル画像を作成します。”)
# サンプル画像として簡単なグラデーション画像を作成
width, height = 200, 150
img_matrix = np.zeros((height, width), dtype=np.uint8)
for i in range(height):
img_matrix[i, :] = int(255 * (i / height)) # 縦方向グラデーション
img = Image.fromarray(img_matrix, ‘L’)
img.save(img_path)
print(f”サンプル画像 ‘{img_path}’ を作成しました。”)
print(f”\n元の画像の形状 (height, width): {img_matrix.shape}”)
元の画像を表示
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.imshow(img_matrix, cmap=’gray’)
plt.title(‘Original Image’)
plt.axis(‘off’)
SVDを実行
画像行列は一般に大きいので、full_matrices=False が効率的
U, s, Vh = np.linalg.svd(img_matrix, full_matrices=False)
print(f”SVD結果の形状: U={U.shape}, s={s.shape}, Vh={Vh.shape}”)
特異値の累積寄与率を計算し、プロット
explained_variance_ratio = np.cumsum(s2) / np.sum(s2)
plt.subplot(1, 2, 2)
plt.plot(explained_variance_ratio)
plt.xlabel(‘Number of Singular Values (k)’)
plt.ylabel(‘Cumulative Explained Variance Ratio’)
plt.title(‘Cumulative Explained Variance by Singular Values’)
plt.grid(True)
plt.show()
上位 k 個の特異値/特異ベクトルを使って近似画像を再構成
いくつかの異なる k の値で試す
k_values = [1, 5, 10, 20, 50] # 試すランクの値
plt.figure(figsize=(15, len(k_values)*2))
for i, k in enumerate(k_values):
if k > len(s):
k = len(s)
print(f”Warning: k={k_values[i]} は特異値の数 {len(s)} を超えています。k={k} を使用します。”)
# ランク k の近似行列を計算
U_k = U[:, :k]
s_k = s[:k]
Vh_k = Vh[:k, :]
Sigma_k = np.diag(s_k)
img_approx = U_k @ Sigma_k @ Vh_k
# 画素値が0-255の範囲に収まるようにクリッピング
img_approx = np.clip(img_approx, 0, 255).astype(np.uint8)
# 近似画像を表示
plt.subplot(len(k_values), 2, 2*i+1)
plt.imshow(img_approx, cmap='gray')
plt.title(f'Approximate Image (k={k})')
plt.axis('off')
# 圧縮率を計算 (元の行列要素数 / 近似に必要な要素数)
# 元: height * width
# 近似: U_k (height * k) + s_k (k) + Vh_k (k * width)
original_elements = img_matrix.shape[0] * img_matrix.shape[1]
compressed_elements = img_matrix.shape[0] * k + k + k * img_matrix.shape[1]
compression_ratio = original_elements / compressed_elements
# 元の行列との差のノルムを表示（オプション）
error_norm = np.linalg.norm(img_matrix - img_approx)
plt.subplot(len(k_values), 2, 2*i+2)
# エラー画像を表示
error_img = np.abs(img_matrix.astype(float) - img_approx.astype(float))
plt.imshow(error_img, cmap='gray', vmin=0, vmax=255)
plt.title(f'Error (k={k}), Ratio={compression_ratio:.2f}, Error Norm={error_norm:.2f}')
plt.axis('off')
plt.tight_layout()
plt.show()
“`
このコードを実行すると、元の画像、特異値の累積寄与率プロット、そして異なるランク $k$ で再構成された近似画像が表示されます。累積寄与率プロットを見ると、少数の大きな特異値が全体の情報の大部分を説明していることがわかります。ランク $k$ を大きくするほど、近似画像は元の画像に近づきますが、圧縮率は低下します。適切な $k$ を選択することで、画質と圧縮率のバランスを取ることができます。エラー画像を見ると、小さな $k$ の場合は大きな構造が捉えきれていないのに対し、大きな $k$ では細かいディテールが再現されていることがわかります。
カラー画像の場合は、通常R, G, Bの各チャンネル（それぞれが高さx幅の行列）に対して独立にSVDを適用し、それぞれを低ランク近似して再構成します。
4.3. 疑似逆行列（Pseudoinverse）の計算
正方行列でない行列や、正方行列であっても正則でない（逆行列を持たない）行列に対しては、通常の意味での逆行列は存在しません。しかし、多くの応用で「逆行列のようなもの」が必要になることがあります。そこで使われるのが疑似逆行列（Pseudoinverse）、あるいはMoore-Penrose逆行列です。
行列 $A$ の疑似逆行列 $A^+$ は、SVDを用いて以下のように定義されます。
$A = U \Sigma V^T$ のとき、$A^+ = V \Sigma^+ U^T$
ここで、$\Sigma^+$ は $\Sigma$ の疑似逆行列で、これは $\Sigma$ の非ゼロの特異値 $\sigma_i$ をその逆数 $1/\sigma_i$ に置き換え、転置した行列です。
NumPyには、疑似逆行列を計算するための関数 numpy.linalg.pinv がありますが、その内部実装はSVDを利用しています。numpy.linalg.pinv(A) は、svd(A) を計算し、得られた特異値と特異ベクトルを使って上記の式で $A^+$ を計算します。
疑似逆行列は、線形方程式 $Ax=b$ の最小二乗解を求める際などに非常に有用です。特に、解が複数存在する場合の最小ノルム解や、解が存在しない場合の最小二乗誤差解を求める際に使われます。$x = A^+ b$ がその解となります。
“`python
例として、線形従属な行を持つ行列を作成 (ランク落ちしている)
A_singular = np.array([[1, 2],
[2, 4],
[3, 6]]) # ランク1の行列
print(“\n特異な行列 A_singular:”)
print(A_singular)
NumPyの pinv を使って疑似逆行列を計算
A_singular_pinv = np.linalg.pinv(A_singular)
print(“\nA_singular の疑似逆行列:”)
print(A_singular_pinv)
SVDを使って疑似逆行列を計算し、結果を比較
U, s, Vh = np.linalg.svd(A_singular)
特異値 s を使って Sigma+ を作成
非ゼロの特異値の逆数を取る
小さすぎる特異値はゼロとみなす閾値を設けるのが一般的 (NumPy pinv の tolerance に相当)
ここでは簡単のため、ゼロでないものをそのまま使う (実際は許容誤差を設定すべき)
s_plus = np.zeros(A_singular.T.shape) # Sigma+ は n x m 形状
tolerance = np.max(s) * max(A_singular.shape) * np.finfo(s.dtype).eps # NumPy pinv のデフォルトに近い計算
s_reciprocal = np.zeros(s.shape)
s_reciprocal[s > tolerance] = 1.0 / s[s > tolerance] # 閾値より大きい特異値の逆数
Sigma_plus = np.zeros(A_singular.shape[::-1]) # A_singular.shape[::-1] は (n, m) 形状
Sigma_plus[:s_reciprocal.shape[0], :s_reciprocal.shape[0]] = np.diag(s_reciprocal)
疑似逆行列 A+ を計算: A+ = V @ Sigma+ @ U.T
NumPy svd は Vh (=V.T) を返すので、V = Vh.T
A_singular_pinv_svd = Vh.T @ Sigma_plus @ U.T
print(“\nSVDから計算した A_singular の疑似逆行列:”)
print(A_singular_pinv_svd)
両者の差を確認
print(“\nNumPy pinv と SVD から計算した疑似逆行列の差の最大絶対値:”)
print(np.max(np.abs(A_singular_pinv – A_singular_pinv_svd)))
``pinv
NumPyので計算した結果と、SVDを使って手動で計算した結果がほぼ一致することが確認できます。実際にはpinv` 関数を使う方が、数値的な安定性などを考慮して実装されているため推奨されます。しかし、SVDの理解という点では、疑似逆行列が特異値と特異ベクトルから構成されることを知っておくことは重要です。
4.4. 主成分分析 (PCA) との関連性
前述のように、PCAはデータ行列の中心化バージョンに対するSVD（あるいは共分散行列の固有値分解）と密接に関連しています。ここでは、NumPyのSVDを使ってPCAの主要なステップを実行する方法を示します。
“`python
サンプルデータの作成 (例えば、20個のデータ点、3つの特徴量)
np.random.seed(42)
data = np.random.rand(20, 3) * 10
特徴量間に線形相関を作るため、少し操作を加える
data[:, 1] = data[:, 0] * 0.5 + np.random.randn(20) # Feature 1 = 0.5 * Feature 0 + Noise
data[:, 2] = data[:, 0] * 1.2 + np.random.randn(20) # Feature 2 = 1.2 * Feature 0 + Noise
print(“\n元のデータ行列 (20×3):”)
print(data)
1. データの中心化 (各特徴量の平均を引く)
mean = np.mean(data, axis=0)
data_centered = data – mean
print(“\n中心化されたデータ行列:”)
print(data_centered)
2. 中心化されたデータ行列に対してSVDを適用
SVD (data_centered) = U @ Sigma @ Vh
ここで Vh の行 (= V の列) が主成分の方向ベクトルとなる
U, s, Vh = np.linalg.svd(data_centered)
print(f”\n中心化データに対するSVD結果: U={U.shape}, s={s.shape}, Vh={Vh.shape}”)
3. 特異値 s は、主成分によって説明される分散の平方根に比例
s の2乗が、対応する主成分によって説明される分散（特異値の2乗は A.T @ A の固有値）
explained_variance = s**2
total_variance = np.sum(explained_variance)
explained_variance_ratio = explained_variance / total_variance
print(“\n説明される分散 (特異値の2乗):”, explained_variance)
print(“合計分散:”, total_variance)
print(“説明される分散比率:”, explained_variance_ratio)
print(“累積説明分散比率:”, np.cumsum(explained_variance_ratio))
主成分は Vh の行ベクトル (V の列ベクトル)
Vh は転置されているので、V を得るには Vh.T を使う
principal_components = Vh.T # V は D x D 行列 (3×3)
print(“\n主成分の方向ベクトル (V の列):”)
print(principal_components)
4. データを上位 k 個の主成分空間に射影 (次元削減)
例えば、上位 2 つの主成分を使う (k=2)
k = 2
principal_components_k = principal_components[:, :k] # V の最初の k 列
print(f”\n上位 {k} 個の主成分ベクトル (V の最初の {k} 列):”)
print(principal_components_k)
中心化されたデータを上位 k 個の主成分に射影
射影されたデータ = data_centered @ principal_components_k
data_reduced = data_centered @ principal_components_k
print(f”\n{k} 次元に削減されたデータ (主成分スコア):”)
print(data_reduced)
5. 必要であれば、元の特徴量空間に戻す（近似）
data_reconstructed = data_reduced @ principal_components_k.T + mean
data_reconstructed = (data_reduced @ principal_components_k.T) + mean
print(f”\n{k} 次元情報を使って再構築された元のデータ:”)
print(data_reconstructed)
比較として、sklearn の PCA と結果を比べてみる
from sklearn.decomposition import PCA
pca = PCA(n_components=k)
data_reduced_sklearn = pca.fit_transform(data)
print(“\nscikit-learn PCA の結果:”)
print(data_reduced_sklearn)
Note: scikit-learn の PCA の主成分の符号は SVD のそれと異なる場合がありますが、方向は同じです。
また、scikit-learn は fit_transform で中心化も含めて処理します。
“`
このコードは、データの中心化、SVDの実行、特異値から説明される分散比率の計算、主成分方向の取得、そしてデータの上位主成分空間への射影（次元削減）という、PCAの主要なステップをNumPyのSVD機能だけで実現しています。特異値の大きさは、対応する主成分がどれだけ多くの分散を説明しているかを示しており、次元削減において保持すべき主成分の数を決定するのに役立ちます。累積説明分散比率を見ることで、例えば「全分散の95%を説明するためにはいくつの主成分が必要か」といった判断ができます。
4.5. レコメンデーションシステム（協調フィルタリング）への応用（概念）
SVDは、レコメンデーションシステムにおける協調フィルタリング手法、特に「行列分解（Matrix Factorization）」モデルの基盤としても利用されます。ユーザーとアイテムの評価行列 $R$（ユーザーを行、アイテムを列とする行列）は、通常非常に大きく、かつほとんどの要素が欠損（評価されていない）した疎行列です。行列分解モデルでは、この評価行列 $R$ を、潜在的な特徴（因子）空間におけるユーザー行列 $U$ とアイテム行列 $V$ の積として近似します。
$R \approx P Q^T$
ここで、$P$ はユーザーを行、潜在因子を列とする行列、$Q$ はアイテムを行、潜在因子を列とする行列です。$P Q^T$ の $(i, j)$ 成分は、ユーザー $i$ がアイテム $j$ を評価するであろう予測値となります。
SVDは理論的には任意の行列に適用できますが、評価行列のように非常に疎で大規模な行列にNumPyの linalg.svd を直接適用するのは、メモリや計算コストの観点から非現実的です。しかし、SVDの考え方（行列を特異値/特異ベクトルで分解する）は、ユーザーとアイテムを低次元の潜在因子空間に埋め込むという行列分解モデルの核心をなしています。実際にレコメンデーションシステムで使われる行列分解アルゴリズム（例: Funk SVD、NMFなど）は、SVDから派生した最適化手法（例: 確率的勾配降下法、交互最小二乗法）を用いて、欠損値がある疎行列に対して近似的な行列分解を行います。
したがって、NumPyの svd はレコメンデーションシステムの直接的な実装には向きませんが、その背後にある数学的概念（行列分解、低ランク近似、潜在因子の抽出）を理解する上で非常に重要です。
5. NumPy SVDの高度な利用と注意点
5.1. 大規模データと計算効率
NumPyの linalg.svd は、密な行列に対しては非常に効率的にSVDを計算できます。しかし、行列のサイズが非常に大きくなると、計算時間とメモリ使用量が増大します。特に、$m \times n$ 行列に対してFull SVDを計算する場合、計算量は概ね $O(\min(m^2 n, m n^2))$ のオーダーになります。また、大きな $U$ ($m \times m$) と $V^T$ ($n \times n$) をメモリに保持する必要があります。
大規模な疎行列の場合、numpy.linalg.svd は効率が悪くなります。このような場合は、SciPyライブラリの scipy.sparse.linalg サブモジュールに含まれるSVD関数を利用することを検討してください。例えば、scipy.sparse.linalg.svds は、疎行列に対して特定の数の最大の特異値とそれに対応する特異ベクトルを効率的に計算できます。
“`python
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import svds
サンプル疎行列を作成
例えば、1000×1000 の疎行列で、要素は1%だけ非ゼロ
sparse_matrix_shape = (1000, 1000)
nnz = 10000 # 非ゼロ要素数
row_indices = np.random.randint(0, sparse_matrix_shape[0], nnz)
col_indices = np.random.randint(0, sparse_matrix_shape[1], nnz)
values = np.random.rand(nnz)
A_sparse = csc_matrix((values, (row_indices, col_indices)), shape=sparse_matrix_shape)
print(f”\n作成した疎行列の形状: {A_sparse.shape}”)
print(f”非ゼロ要素数: {A_sparse.nnz}”)
密行列に変換してNumPyのsvdを試す (注意: メモリを大量消費する可能性あり)
A_dense = A_sparse.toarray()
try:
U_dense, s_dense, Vh_dense = np.linalg.svd(A_dense)
print(f”NumPy SVD (密行列) 結果の形状: U={U_dense.shape}, s={s_dense.shape}, Vh={Vh_dense.shape}”)
except MemoryError:
print(“NumPy SVD (密行列) でメモリ不足が発生しました。”)
SciPyの svds を使って上位 k 個の特異値/ベクトルを計算
k_sparse = 10 # 計算する特異値の数
try:
U_sparse, s_sparse, Vh_sparse = svds(A_sparse, k=k_sparse)
print(f”\nSciPy SVDS (疎行列) 結果の形状 (上位 {k_sparse} 個):”)
print(f”U={U_sparse.shape}, s={s_sparse.shape}, Vh={Vh_sparse.shape}”)
# svds の s は昇順で返されるため、降順に並べ替える
s_sparse = s_sparse[::-1]
# U と Vh の列/行も s に合わせて並べ替える必要がある
U_sparse = U_sparse[:, ::-1]
Vh_sparse = Vh_sparse[::-1, :]
print(“SciPy SVDS の特異値 (降順):”, s_sparse)
except Exception as e:
print(f”SciPy SVDS の実行中にエラーが発生しました: {e}”)
``scipy.sparse.linalg.svdsは、計算する特異値の数kを指定できるため、必要な情報だけを効率的に抽出できます。これは、大規模なデータセットにおけるSVDの適用範囲を広げます。注意点として、svdsは特異値を**昇順**で返すため、NumPyのsvd` と同様に降順で扱いたい場合は、戻り値を並べ替える必要があります。
5.2. 数値的安定性と許容誤差
浮動小数点演算には常に誤差が伴います。SVD計算も例外ではありません。特に非常に小さい特異値は、計算機上ではゼロと区別がつかなくなることがあります。疑似逆行列の計算や、非ゼロの特異値の数を数えて行列のランクを推定する場合など、小さい特異値をどのように扱うか（ゼロとみなすか）が重要になります。
numpy.linalg.pinv のように、多くの数値線形代数関数は、このために許容誤差（tolerance）パラメータを持っています。特異値がこの許容誤差以下の場合、その特異値はゼロとみなされます。手動でSVDの結果を利用して処理を行う場合も、同様の考慮が必要です。
5.3. 特異ベクトルの一意性
特異値 $\sigma_i$ が全て異なる場合、対応する特異ベクトル $u_i$ と $v_i$ の方向（符号を除いて）は一意的です。しかし、複数の特異値が等しい場合（縮退している場合）、それらに対応する特異ベクトルは、それらが張る部分空間内での任意の正規直交基底として選ぶことができます。NumPyを含むSVDアルゴリズムの実装は、このような場合にある特定の基底を選択しますが、その選択はアルゴリズムによって異なる可能性があり、必ずしも一意的ではありません。したがって、特異値が等しい場合の特異ベクトルの具体的な値や符号が、異なる環境やNumPyのバージョンでわずかに異なることがあっても驚く必要はありません。これは特異ベクトルが張る部分空間が同じであれば、多くの応用においては問題になりません。
6. まとめ
本記事では、PythonのNumPyライブラリを用いた特異値分解（SVD）の実装方法について、その数学的な基礎から具体的なコード例、そして様々な応用例まで、詳細に解説しました。
- SVDは、任意の行列 $A$ を $U \Sigma V^T$ の形に分解する強力な手法です。
- $U$ と $V^T$ は直交行列（回転）、$\Sigma$ は対角成分に特異値を並べた行列（スケーリング）を表します。
- 特異値は行列の情報を圧縮したものであり、大きい特異値ほどデータにおける主要なパターンを捉えています。
- NumPyの
numpy.linalg.svd関数を使うことで、手軽にSVDを計算できます。戻り値U,s,Vhを適切に扱うことで、元の行列の復元や、様々な応用が可能になります。 full_matrices引数でFull SVDとReduced SVDを切り替えられます。多くの応用では計算効率の良いReduced SVDが適しています。compute_uv引数で特異値のみの計算が可能です。- SVDは、低ランク近似、次元削減、画像圧縮・ノイズ除去、疑似逆行列の計算、そしてPCAなどの基盤となります。これらの応用をNumPyのSVD結果を使って実装する方法を示しました。
- 大規模な疎行列に対しては、SciPyの
scipy.sparse.linalg.svdsのような専門的な関数がより効率的です。 - 数値計算上の注意点として、小さい特異値の扱い（許容誤差）や、等しい特異値に対応する特異ベクトルの一意性について考慮が必要です。
SVDは線形代数の最も美しい結果の一つであり、その応用範囲は広大です。NumPyを使うことで、この強力な数学的ツールをPython上で簡単に利用することができます。本記事が、SVDの理論的な理解を深めるとともに、NumPyを使った実践的なスキルを習得するための一助となれば幸いです。
7. 付録：必要なライブラリのインポート
記事中の全てのコード例を実行するために必要なライブラリは以下の通りです。
“`python
import numpy as np
画像処理の例を実行する場合
try:
from PIL import Image
import matplotlib.pyplot as plt
except ImportError:
print(“Pillow および Matplotlib がインストールされていません。画像処理の例は実行できません。”)
print(“インストールするには: pip install Pillow matplotlib”)
SciPy sparse.linalg の例を実行する場合
try:
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import svds
except ImportError:
print(“SciPy がインストールされていません。SciPy sparse.linalg の例は実行できません。”)
print(“インストールするには: pip install scipy”)
“`
必要に応じてこれらのライブラリをインストールしてください。