【線形代数】特異値分解をゼロからわかりやすく解説！

どうも。こんにちは。
ケミカルエンジニアのこーしです。

本日は、「特異値分解（SVD：Singular Value Decomposition）」をゼロからわかりやすく解説していきます！

線形代数の教科書では「固有値分解」を中心に学びますが、scikit-learn の PCA や PLS の実装を覗いてみると、内部では特異値分解が使われていることがわかります。

これはなぜでしょうか？

本記事では、特異値分解についてゼロからわかりやすく解説するとともに、固有値分解と特異値分解の違いについて詳しく解説します。

なぜ特異値分解がプログラムの実装で多用されているのかを紐解いていきましょう！

この記事を書いた人

こーし(@mimikousi)

特異値分解とは

\( n \times m \) の任意の行列 \( A \) に対して、

$$
A = U \Sigma V^\top
$$

の形に分解することを特異値分解（SVD）と呼びます。

\( \Sigma \) の対角成分 \( \sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r \geq 0 \) を特異値と呼び、\( r \) は \( A \) の階数（ランク）を表します。

固有値分解 \( A = P \Lambda P^{-1} \) との大きな違いは、以下の3点です。

 

特異値分解の意味

固有値分解と同じく、特異値分解も行列 \( A \) による変換を3つのステップに分解しています。

$$
A\boldsymbol{x} = U\Sigma V^\top \boldsymbol{x}
$$

つまり、どんな行列の変換も「回転 → 伸縮 → 回転」という3ステップに分解できる、というのが特異値分解の本質的なメッセージです。

固有値分解では入力と出力で同じ基底（\( P \)）を使うため正方行列に限定されますが、特異値分解では入力側 \( V \) と出力側 \( U \) で別々の基底を使うため、長方形の行列でも自然に分解できるのです。

なお、ステップ1で \( V^\top \) を掛ける操作は「\( \boldsymbol{x} \) を右特異ベクトル基底で表したときの座標」を取り出す操作です。

例えば \( \boldsymbol{x} = c_1  \boldsymbol{v}_1 + \cdots + c_m \boldsymbol{v}_m\) と表したいとき、行列で書くと \( \boldsymbol{x} = V\boldsymbol{c} \) なので、両辺に左から \( V^\top \) を掛けて \( \boldsymbol{c} = V^\top \boldsymbol{x} \) となります（\( V^\top V = I \)）。

前回の固有値分解の記事で解説した \( P^{-1}\boldsymbol{x} \) と完全に同じ役割を果たしています。

 

右特異ベクトルと左特異ベクトルの関係

特異値分解の各ベクトルの関係を、定義式から直接導いてみましょう。

SVDの定義式 \(A = U \Sigma V^\top\) の両辺に右から \(\boldsymbol{v}_j\)（\(V\) の第 \(j\) 列）を掛けます。

$$
A \boldsymbol{v}_j = U \Sigma V^\top \boldsymbol{v}_j
$$

右辺を順に計算していきましょう。

変形①：\(V^\top \boldsymbol{v}_j\) を計算する

\(V\) は直交行列（\(V^\top V = I\)）で、第 \(i\) 列は \(\boldsymbol{v}_i\) なので、

$$
V^\top \boldsymbol{v}_j = \begin{pmatrix} \boldsymbol{v}_1^\top \\ \boldsymbol{v}_2^\top \\ \vdots \\ \boldsymbol{v}_m^\top \end{pmatrix} \boldsymbol{v}_j = \begin{pmatrix} \boldsymbol{v}_1^\top \boldsymbol{v}_j \\ \boldsymbol{v}_2^\top \boldsymbol{v}_j \\ \vdots \\ \boldsymbol{v}_m^\top \boldsymbol{v}_j \end{pmatrix}
$$

正規直交性より、

$$
\boldsymbol{v}_i^\top \boldsymbol{v}_j = \begin{cases} 1 & (i = j) \\ 0 & (i \neq j) \end{cases}
$$

なので、\(j\) 番目の成分だけ \(1\)、他はすべて \(0\) となり、

$$
V^\top \boldsymbol{v}_j = \begin{pmatrix} 0 \\ \vdots \\ 1 \\ \vdots \\ 0 \end{pmatrix} = \boldsymbol{e}_j
$$

つまり、\(V^\top\) で写すと \(\boldsymbol{v}_j\) は標準基底ベクトル \(\boldsymbol{e}_j\) になります。

 

変形②：\(\Sigma \boldsymbol{e}_j\) を計算する

対角行列 \(\Sigma\) に \(\boldsymbol{e}_j\) を掛けることは、\(\Sigma\) の第 \(j\) 列を取り出す操作です。

$$
\Sigma \boldsymbol{e}_j = \begin{pmatrix} 0 \\ \vdots \\ \sigma_j \\ \vdots \\ 0 \end{pmatrix} = \sigma_j \boldsymbol{e}_j
$$

 

変形③：\(U(\sigma_j \boldsymbol{e}_j)\) を計算する

行列 \(U\) に標準基底ベクトル \(\boldsymbol{e}_j\) を掛けると、\(U\) の第 \(j\) 列 \(\boldsymbol{u}_j\) が取り出されます。

$$
U(\sigma_j \boldsymbol{e}_j) = \sigma_j \cdot U\boldsymbol{e}_j = \sigma_j \boldsymbol{u}_j
$$

 

まとめ：基本関係式

以上をつなげると、

$$
A \boldsymbol{v}_j = \sigma_j \boldsymbol{u}_j
$$

これは固有値分解の基本関係式 \(A\boldsymbol{v} = \lambda \boldsymbol{v}\) の長方形行列版です。

固有値分解では「入力と同じ方向のまま \(\lambda\) 倍」されますが、特異値分解では「入力空間の \(\boldsymbol{v}_j\) が、出力空間の \(\boldsymbol{u}_j\) の方向に \(\sigma_j\) 倍されて出てくる」のです。

入力と出力で別々の基底を使うため、長方形行列でも自然に成立する関係式となります。

 

左特異ベクトルの求め方

基本関係式 \(A\boldsymbol{v}_j = \sigma_j \boldsymbol{u}_j\) の両辺を \(\sigma_j\) で割ると（\(\sigma_j > 0\) のとき）、

$$
\boldsymbol{u}_j = \frac{A \boldsymbol{v}_j}{\sigma_j}
$$

つまり、右特異ベクトル \(\boldsymbol{v}_j\) と特異値 \(\sigma_j\) さえ求めれば、左特異ベクトル \(\boldsymbol{u}_j\) は \(A\) で写すだけで自動的に得られるのです。

ただし、\(\sigma_j = 0\) のときは0で割れないため、別の方法（\(AA^\top\) の固有ベクトルとして求める方法）が必要になります。

これは具体例で扱います。

 

特異値分解の求め方

特異値分解を実際に計算するには、対称行列 \(A^\top A\) と \(AA^\top\) の固有値分解を使います。

なぜこの2つの行列が出てくるのか、SVDの定義式から導いてみましょう。

 

右特異ベクトルは \(A^\top A\) の固有ベクトル

SVDの定義式 \(A = U \Sigma V^\top\) から、\(A^\top A\) を計算すると、

$$
\begin{aligned} A^\top A &= (U \Sigma V^\top)^\top (U \Sigma V^\top) \\[5pt] &= V \Sigma^\top U^\top \cdot U \Sigma V^\top \\[5pt] &= V \Sigma^\top \Sigma V^\top \quad (U^\top U = I) \end{aligned}
$$

\(\Sigma^\top \Sigma\) は \(m \times m\) の対角行列で、対角成分は \(\sigma_1^2, \sigma_2^2, \ldots, \sigma_r^2, 0, \ldots, 0\) となります。

これは前回学んだ対称行列のスペクトル分解 \(A = V \Lambda V^\top\) そのものの形なので、

 

左特異ベクトルは \(AA^\top\) の固有ベクトル

同様に、\(AA^\top\) を計算すると、

$$
\begin{aligned} AA^\top &= (U \Sigma V^\top)(U \Sigma V^\top)^\top \\[5pt] &= U \Sigma V^\top \cdot V \Sigma^\top U^\top \\[5pt] &= U \Sigma \Sigma^\top U^\top \quad (V^\top V = I) \end{aligned}
$$

\(\Sigma \Sigma^\top\) は \(n \times n\) の対角行列で、対角成分は \(\sigma_1^2, \sigma_2^2, \ldots, \sigma_r^2, 0, \ldots, 0\) となります。

こちらも対称行列のスペクトル分解の形なので、

 

\(A^\top A\) と \(AA^\top\) の固有値の関係

ここで興味深い性質に気づきます。

\(A^\top A\) は \(m \times m\) 行列で固有値が \(m\) 個、\(AA^\top\) は \(n \times n\) 行列で固有値が \(n\) 個ありますが、

これはなぜでしょうか？

\(A^\top A\) の固有値 \(\lambda\) と固有ベクトル \(\boldsymbol{v} \neq \boldsymbol{0}\) を考えます。

$$
A^\top A \boldsymbol{v} = \lambda \boldsymbol{v}
$$

両辺に左から \(A\) を掛けると、

$$
A A^\top (A \boldsymbol{v}) = \lambda (A \boldsymbol{v})
$$

ここで \(A\boldsymbol{v} \neq \boldsymbol{0}\) であれば（\(\lambda \neq 0\) のときに成立）、

となります。

逆方向（\(AA^\top\) の固有ベクトル \(\boldsymbol{u}\) から \(A^\top \boldsymbol{u}\) が \(A^\top A\) の固有ベクトルになる）も同様に示せます。

つまり、\(A^\top A\) と \(AA^\top\) の非ゼロ固有値は1対1で対応するということです。

\(m \neq n\) のときは行列のサイズが違うので、余った分の固有値は \(0\) になります。

例えば \(n > m\) なら、\(AA^\top\)（\(n \times n\) 行列）の固有値のうち \(m\) 個は \(A^\top A\) と共通で、残り \(n - m\) 個は \(0\) です。

念のため、\(A\boldsymbol{v}_i \neq \boldsymbol{0}\)（\(\lambda_i \neq 0\) のとき成立）を確認しておきます。

$$
\| A\boldsymbol{v}_i \|^2 = \boldsymbol{v}_i^\top A^\top A \boldsymbol{v}_i = \lambda_i \boldsymbol{v}_i^\top \boldsymbol{v}_i = \lambda_i
$$

つまり \(\lambda_i > 0\) なら \(\| A\boldsymbol{v}_i \| = \sqrt{\lambda_i} = \sigma_i > 0\) なので、\(A\boldsymbol{v}_i\) は確かに非ゼロベクトルです。

このノルムの計算結果から、\(\boldsymbol{u}_i = A\boldsymbol{v}_i / \sigma_i\) を作ると単位ベクトルになることもわかります。

 

求め方のまとめ

以上を踏まえると、特異値分解は次の手順で求められます。

なお、数値安定性の観点から、実際のプログラムではこのように \(A^\top A\) を経由せず、\(A\) を直接ハウスホルダー変換などで分解します。

ここではあくまで「特異値分解の中身を理解する」ための手順として紹介します。

 

具体例（3×2 行列）

それでは、長方形の行列を使って実際に特異値分解を計算してみましょう。

$$
A = \begin{pmatrix} 1 & 1 \\ 1 & 0 \\ 0 & 1 \end{pmatrix}
$$

\(A\) は \(3 \times 2\) の長方形行列なので、固有値分解は使えませんが、特異値分解は問題なく適用できます。

 

Step 1：\(A^\top A\) を計算する

$$
\begin{aligned} A^\top A &= \begin{pmatrix} 1 & 1 & 0 \\ 1 & 0 & 1 \end{pmatrix}\begin{pmatrix} 1 & 1 \\ 1 & 0 \\ 0 & 1 \end{pmatrix} \\[5pt] &= \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \end{aligned}
$$

 

Step 2：\(A^\top A\) の固有値・固有ベクトルを求める

特性方程式を立てます。

$$
\begin{aligned}
\det(A^\top A - \lambda I) &= \det\begin{pmatrix} 2-\lambda & 1 \\ 1 & 2-\lambda \end{pmatrix}\\[5pt]
&= (2-\lambda)^2 - 1\\[5pt]
&= 0
\end{aligned}
$$

$$
(\lambda - 1)(\lambda - 3) = 0 \quad \Rightarrow \quad \lambda_1 = 3, \quad \lambda_2 = 1
$$

（特異値の慣例に従い、大きい順に並べます。）

\(\lambda_1 = 3\) のとき：

$$
(A^\top A - 3I)\boldsymbol{v} = \begin{pmatrix} -1 & 1 \\ 1 & -1 \end{pmatrix}\boldsymbol{v} = \boldsymbol{0}
$$

第1行より \(-x_1 + x_2 = 0\)、すなわち \(x_2 = x_1\)。 \(x_1 = 1\) とおいて正規化すると、

$$
\boldsymbol{v}_1 = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ 1 \end{pmatrix}
$$

\(\lambda_2 = 1\) のとき：

$$
(A^\top A - I)\boldsymbol{v} = \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}\boldsymbol{v} = \boldsymbol{0}
$$

第1行より \(x_1 + x_2 = 0\)、すなわち \(x_2 = -x_1\)。 \(x_1 = 1\) とおいて正規化すると、

$$
\boldsymbol{v}_2 = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ -1 \end{pmatrix}
$$

 

Step 3：特異値を求める

$$
\sigma_1 = \sqrt{\lambda_1} = \sqrt{3}, \quad \sigma_2 = \sqrt{\lambda_2} = 1
$$

 

Step 4：左特異ベクトル \(\boldsymbol{u}_1, \boldsymbol{u}_2\) を求める

特異値 \(\sigma_1, \sigma_2\) はどちらも正なので、先ほど導いた基本関係式 \(A\boldsymbol{v}_j = \sigma_j \boldsymbol{u}_j\) を変形した

$$
\boldsymbol{u}_j = \frac{A \boldsymbol{v}_j}{\sigma_j}
$$

をそのまま使えます。

\(\boldsymbol{u}_1\)：

$$
\begin{aligned} A \boldsymbol{v}_1 &= \begin{pmatrix} 1 & 1 \\ 1 & 0 \\ 0 & 1 \end{pmatrix} \cdot \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ 1 \end{pmatrix} \\[5pt] &= \frac{1}{\sqrt{2}}\begin{pmatrix} 2 \\ 1 \\ 1 \end{pmatrix} \end{aligned}
$$

$$
\boldsymbol{u}_1 = \frac{A\boldsymbol{v}_1}{\sigma_1} = \frac{1}{\sqrt{2}\cdot\sqrt{3}}\begin{pmatrix} 2 \\ 1 \\ 1 \end{pmatrix} = \frac{1}{\sqrt{6}}\begin{pmatrix} 2 \\ 1 \\ 1 \end{pmatrix}
$$

\(\boldsymbol{u}_2\)：

$$
\begin{aligned} A \boldsymbol{v}_2 &= \begin{pmatrix} 1 & 1 \\ 1 & 0 \\ 0 & 1 \end{pmatrix} \cdot \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ -1 \end{pmatrix} \\[5pt] &= \frac{1}{\sqrt{2}}\begin{pmatrix} 0 \\ 1 \\ -1 \end{pmatrix} \end{aligned}
$$

$$
\boldsymbol{u}_2 = \frac{A\boldsymbol{v}_2}{\sigma_2} = \frac{1}{\sqrt{2}}\begin{pmatrix} 0 \\ 1 \\ -1 \end{pmatrix}
$$

 

Step 4-2：\(\boldsymbol{u}_3\) を求める

ここまでで特異値2つ分の左特異ベクトル \(\boldsymbol{u}_1, \boldsymbol{u}_2\) が得られました。

しかし、\(U\) は \(3 \times 3\) の直交行列にしたいため、もう1本 \(\boldsymbol{u}_3\) が必要です。

このとき対応する特異値は \(\sigma_3 = 0\) となるため、基本関係式 \(\boldsymbol{u}_j = A\boldsymbol{v}_j / \sigma_j\) ではゼロ除算となり計算できません。

そこで、先ほど示した「左特異ベクトルは \(AA^\top\) の固有ベクトル」という性質を使い、\(\boldsymbol{u}_3\) を\(AA^\top\) の固有値 \(0\) に対応する固有ベクトルとして求めます。

まず \(AA^\top\) を計算すると、

$$
\begin{aligned} AA^\top &= \begin{pmatrix} 1 & 1 \\ 1 & 0 \\ 0 & 1 \end{pmatrix}\begin{pmatrix} 1 & 1 & 0 \\ 1 & 0 & 1 \end{pmatrix} \\[5pt] &= \begin{pmatrix} 2 & 1 & 1 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \end{pmatrix} \end{aligned}
$$

固有値 \(\lambda = 0\) に対応する固有ベクトルは \(AA^\top \boldsymbol{u} = \boldsymbol{0}\) を満たすので、

$$
\begin{pmatrix} 2 & 1 & 1 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \end{pmatrix}\boldsymbol{u} = \boldsymbol{0}
$$

第2行：\(x_1 + x_2 = 0 \Rightarrow x_2 = -x_1\)

第3行：\(x_1 + x_3 = 0 \Rightarrow x_3 = -x_1\)
（第1行 \(2x_1 + x_2 + x_3 = 2x_1 - x_1 - x_1 = 0\) も満たすので矛盾無し）

\(x_1 = 1\) とおいて正規化すると、

$$
\boldsymbol{u}_3 = \frac{1}{\sqrt{3}}\begin{pmatrix} 1 \\ -1 \\ -1 \end{pmatrix}
$$

念のため、\(\boldsymbol{u}_1, \boldsymbol{u}_2\) と直交していることを検算してみましょう。

$$
\boldsymbol{u}_1^\top \boldsymbol{u}_3 = \frac{1}{\sqrt{6}\sqrt{3}}(2 - 1 - 1) = 0 \quad
$$

$$
\boldsymbol{u}_2^\top \boldsymbol{u}_3 = \frac{1}{\sqrt{2}\sqrt{3}}(0 - 1 + 1) = 0 \quad
$$

きちんと直交していますね。

これで \(\boldsymbol{u}_1, \boldsymbol{u}_2, \boldsymbol{u}_3\) が \(\mathbb{R}^3\) の正規直交基底となり、\(U\) を直交行列として組み立てられます。

上述の「特異値分解の求め方」で示した手順の通り、\(A^\top A\) から右特異ベクトルと特異値を求め、基本関係式から左特異ベクトル \(\boldsymbol{u}_1, \boldsymbol{u}_2\) を導き、\(AA^\top\) から残りの \(\boldsymbol{u}_3\) を補う、という流れで特異値分解が完成しました。

 

Step 5：\(U, \Sigma, V\) を組み立てる

$$
U = \begin{pmatrix} \frac{2}{\sqrt{6}} & 0 & \frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{2}} & \frac{-1}{\sqrt{3}} \\ \frac{1}{\sqrt{6}} & \frac{-1}{\sqrt{2}} & \frac{-1}{\sqrt{3}} \end{pmatrix}
$$

$$
\Sigma = \begin{pmatrix} \sqrt{3} & 0 \\ 0 & 1 \\ 0 & 0 \end{pmatrix}, \quad V = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}
$$

\(\Sigma\) の3行目はすべて0です。

これは「3次元空間に埋め込んだとき、入力（2次元）からは到達できない方向がある」ことを表しています。

このような \(\Sigma\) の0行を含めない、コンパクトな形のSVD（thin SVD）もあり、scikit-learn の

np.linalg.svd(A, full_matrices=False)

ではこちらが返されます。

 

検算（\(U \Sigma V^\top = A\) の確認）

まず \(\Sigma V^\top\) を計算します。

$$
V^\top = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}
$$

$$
\begin{aligned}
\Sigma V^\top &= \begin{pmatrix} \sqrt{3} & 0 \\ 0 & 1 \\ 0 & 0 \end{pmatrix} \cdot \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\\[5pt]
&= \frac{1}{\sqrt{2}}\begin{pmatrix} \sqrt{3} & \sqrt{3} \\ 1 & -1 \\ 0 & 0 \end{pmatrix}
\end{aligned}
$$

次に \(U(\Sigma V^\top)\) を計算します。

代表的な要素をチェックしてみましょう。

(1,1) 要素：

$$
\begin{aligned} &\frac{2}{\sqrt{6}} \cdot \frac{\sqrt{3}}{\sqrt{2}} + 0 \cdot \frac{1}{\sqrt{2}} + \frac{1}{\sqrt{3}} \cdot 0 \\[5pt] &= \frac{2\sqrt{3}}{\sqrt{12}} = \frac{2\sqrt{3}}{2\sqrt{3}} = 1 \quad \end{aligned}
$$

(1,2) 要素：

$$
\begin{aligned} &\frac{2}{\sqrt{6}} \cdot \frac{\sqrt{3}}{\sqrt{2}} + 0 \cdot \frac{-1}{\sqrt{2}} + \frac{1}{\sqrt{3}} \cdot 0 \\[5pt] &= \frac{2\sqrt{3}}{2\sqrt{3}} = 1 \quad \end{aligned}
$$

(2,1) 要素：

$$
\begin{aligned} &\frac{1}{\sqrt{6}} \cdot \frac{\sqrt{3}}{\sqrt{2}} + \frac{1}{\sqrt{2}} \cdot \frac{1}{\sqrt{2}} + \frac{-1}{\sqrt{3}} \cdot 0 \\[5pt] &= \frac{\sqrt{3}}{2\sqrt{3}} + \frac{1}{2} = \frac{1}{2} + \frac{1}{2} = 1 \quad \end{aligned}
$$

(2,2) 要素：

$$
\begin{aligned} &\frac{1}{\sqrt{6}} \cdot \frac{\sqrt{3}}{\sqrt{2}} + \frac{1}{\sqrt{2}} \cdot \frac{-1}{\sqrt{2}} + \frac{-1}{\sqrt{3}} \cdot 0 \\[5pt] &= \frac{1}{2} - \frac{1}{2} = 0 \quad \end{aligned}
$$

残りの要素も同様に計算すると、

$$
U \Sigma V^\top = \begin{pmatrix} 1 & 1 \\ 1 & 0 \\ 0 & 1 \end{pmatrix} = A \quad
$$

確かに復元できました。

 

スペクトル分解との対応

前回の記事で、\( n \times n \) の対称行列\(A\)のスペクトル分解は

$$
A = \sum_{i=1}^{n} \lambda_i \boldsymbol{v}_i \boldsymbol{v}_i^\top
$$

と展開できることを学びました。

特異値分解にも同じ表現があります。

$$
A = \sum_{i=1}^{r} \sigma_i \boldsymbol{u}_i \boldsymbol{v}_i^\top
$$

各項 \(\sigma_i \boldsymbol{u}_i \boldsymbol{v}_i^\top\) はランク1の行列で、特異値 \(\sigma_i\) の大きい順から「\(A\) の主要な構造」を表しています。

上位 \(k\) 項までで打ち切れば、

$$
A \approx \sum_{i=1}^{k}\sigma_i \boldsymbol{u}_i \boldsymbol{v}_i^\top
$$

となり、これがエッカート・ヤング（Eckart-Young）の定理で最良の低ランク近似になることが知られています。

 

PCAでスペクトル分解を上位 \(k\) 項で打ち切ったのと、まったく同じ考え方ですね。

 

特異値分解のメリット／デメリット

特異値分解のメリットとデメリットを下記にまとめてみました。

メリットの方が圧倒的に大きく、特に統計学・機械学習の文脈ではほぼすべての場面でSVDが優位に立ちます。

 

なぜ特異値分解が好まれるのか

ここからが本記事のメインテーマです。

なぜ scikit-learn の PCA や PLS は、教科書通りの「共分散行列の固有値分解」ではなく、わざわざ特異値分解を使っているのでしょうか。

理由は次の4つに整理できます。

順に見ていきましょう。

 

理由①：固有値分解は正方行列にしか適用できない

そもそも実用データの多くは長方形の行列です。

例えば化学プラントのプロセスデータでは、

というように、\(n \gg m\) の長方形行列がほとんどです。

固有値分解は定義上、正方行列にしか適用できません。

これには2つの数学的な理由があります。

理由 (a)：固有値・固有ベクトルの定義式が成り立たない

固有値・固有ベクトルの定義式は、

$$
A \boldsymbol{v} = \lambda \boldsymbol{v}
$$

でした。

ここで \(A\) が \(n \times m\) の長方形行列だとすると、

となり、両辺の次元が合わないため、そもそも等式が成立しません。

つまり、\(m \neq n\) では「向きが変わらず長さだけ変わる」という固有ベクトルの概念自体が定義できないのです。

理由 (b)：特性方程式に必要な行列式が存在しない

固有値を求めるには特性方程式 \(\det(A - \lambda I) = 0\) を解く必要がありますが、行列式 \(\det\) は正方行列にしか定義されていません。

長方形行列に対しては \(A - \lambda I\) という式自体が意味を持たないため（単位行列のサイズが決まらない）、特性方程式を立てることができません。

一方、特異値分解は \(A^\top A\)（\(m \times m\) の正方対称行列）の固有値分解を経由するため、長方形の \(A\) であっても問題なく適用できます。

これがSVDの最大の汎用性です。

 

理由②：特異値は常に実数・非負

固有値は一般の正方行列に対して複素数になることがあります。

例えば、

$$
A = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}
$$

は90度回転を表す行列ですが、特性方程式は

$$
\lambda^2 + 1 = 0 \quad \Rightarrow \quad \lambda = \pm i
$$

となり、固有値は虚数です。

これは「90度回転で向きが変わらないベクトルは（実数の世界では）存在しない」という事実と対応しています。

一方、特異値はどんな行列に対しても常に実数かつ非負です。

これを証明してみましょう。

特異値が実数・非負である証明

任意の実行列 \(X \in \mathbb{R}^{n \times m}\) に対して、\(X^\top X\) を考えます。

ステップ①：\(X^\top X\) は対称行列

$$
(X^\top X)^\top = X^\top (X^\top)^\top = X^\top X
$$

なので、\(X^\top X\) は実対称行列です。

前回紹介したスペクトル定理より、実対称行列の固有値はすべて実数です。

ステップ②：\(X^\top X\) は半正定値

任意のベクトル \(\boldsymbol{v}\) に対して、

$$
\boldsymbol{v}^\top (X^\top X) \boldsymbol{v} = (X\boldsymbol{v})^\top (X\boldsymbol{v}) = \| X\boldsymbol{v} \|^2 \geq 0
$$

ノルムの二乗は常に非負なので、\(X^\top X\) は半正定値行列です。

\(X^\top X \boldsymbol{v} = \lambda \boldsymbol{v}\) より、固有値 \(\lambda\) は

$$
\lambda = \frac{\boldsymbol{v}^\top X^\top X \boldsymbol{v}}{\boldsymbol{v}^\top \boldsymbol{v}} = \frac{\| X\boldsymbol{v} \|^2}{\| \boldsymbol{v} \|^2} \geq 0
$$

となり、すべての固有値が非負であることがわかります。

ステップ③：特異値は固有値の平方根

特異値の定義 \(\sigma_i = \sqrt{\lambda_i}\) より、\(\lambda_i\) が実数かつ非負であることから、

$$
\sigma_i = \sqrt{\lambda_i} \in \mathbb{R}, \quad \sigma_i \geq 0
$$

これで、特異値が常に実数かつ非負であることが示せました。

固有値のように複素数を扱う必要がなく、計算機上でもシンプルに扱えます。

これも特異値分解が好まれる大きな理由のひとつです。

 

理由③：計算コストは実用上ほぼ同等

「SVDの方が計算量が多そう」というイメージがありますが、実用上の \(n \gg m\) のケースでは固有値分解とほぼ同じ計算量です。

手法	計算量
固有値分解（\(X^\top X\) 経由）	\(O(nm^2 + m^3)\)
特異値分解（直接）	\(O(nm^2 + m^3)\)
打ち切りSVD（上位 \(k\) 個）	\(O(nmk)\)

\(n \gg m\) のとき \(nm^2\) の項が支配的になるため、固有値分解とSVDのオーダーは同じです。

定数係数を考慮するとSVDの方が若干重いですが、実用上は無視できる差です。

それよりも、scikit-learn の `PCA` や `TruncatedSVD` のように上位 \(k\) 個の特異値だけを計算する打ち切りSVDが使えるため、\(k \ll m\) の場面ではSVDの方が圧倒的に高速です。

例えばプロセスデータで200変数あったとして、上位5主成分だけ欲しい場合、固有値分解は必ず200成分すべて計算してしまいますが、打ち切りSVDなら5成分だけを効率的に計算できます。

 

理由④：数値安定性が圧倒的に高い

これが特異値分解が選ばれる、おそらく最大の理由です。

教科書の説明では、PCAは共分散行列 \(S = \dfrac{1}{n-1} X^\top X\) の固有値分解で求めます。

しかしこの \(X^\top X\) を計算する瞬間に、行列の条件数が二乗されてしまうという致命的な問題が発生します。

条件数とは 行列 \(X\) の条件数 \(\kappa(X)\) は、特異値の最大値と最小値の比

$$
\kappa(X) = \frac{\sigma_{\max}}{\sigma_{\min}}
$$

で定義され、入力の微小な誤差が出力にどれだけ増幅されるかを表します。

\(X^\top X\) を経由すると条件数が二乗される \(X\) の特異値を \(\sigma_1 \geq \cdots \geq \sigma_n \geq 0\) とすると、

$$
\kappa(X) = \frac{\sigma_1}{\sigma_n}
$$

ここで、\(X^\top X\) の固有値は \(\sigma_i^2\) であることを思い出してください（特異値の定義より）。

すると、\(X^\top X\) の条件数は

$$
\kappa(X^\top X) = \frac{\sigma_1^2}{\sigma_n^2} = \left( \frac{\sigma_1}{\sigma_n} \right)^2 = \kappa(X)^2
$$

となり、条件数が二乗されてしまうのです。

 

多重共線性との関係

ここで、化学プラントのプロセスデータの特徴を思い出してください。

多重共線性が強いとは、ある変数がほかの変数の線形結合でほぼ表現できる状態のこと。

線形代数的に言えば、データ行列 \(X\) の列ベクトルが「ほぼ線形従属」、つまり\(X\) が「ほぼランク落ち」している状態です。

このとき、\(X\) の特異値の中には非常に小さいもの（\(\sigma_n \approx 0\) に近い値）が含まれます。

すると、

$$
\kappa(X) = \frac{\sigma_1}{\sigma_n} \quad \text{が大きくなる}
$$

$$
\kappa(X^\top X) = \kappa(X)^2 \quad \text{はさらに大きくなる}
$$

となり、共分散行列の固有値分解は数値的にきわめて不安定になります。

 

有効桁数で見る具体的な影響

倍精度浮動小数点（float64）の精度は約 \(\epsilon \approx 10^{-16}\) です。

数値計算で得られる結果の有効桁数の目安は、

$$
\text{有効桁数} \approx 16 - \log_{10}(\kappa)
$$

となります。

状況	手法	有効桁数
\(\kappa(X) = 10^3\)（軽度の共線性）	SVD直接	13桁
	\(X^\top X\) 経由	10桁
\(\kappa(X) = 10^6\)（強い共線性、プロセス典型）	SVD直接	10桁
	\(X^\top X\) 経由	4桁
\(\kappa(X) = 10^8\)（極度の共線性）	SVD直接	8桁
	\(X^\top X\) 経由	0桁（信頼不能）

プロセスデータで典型的な \(\kappa(X) = 10^6\) のケースを見ると、SVD直接では10桁の精度を保てますが、\(X^\top X\) 経由ではわずか4桁にまで劣化します。

特に、データのばらつきが小さい方向（共線性方向）ほど、共分散行列の固有値分解で復元される第2、第3主成分の方向がわずかにズレるという形で影響が現れます。

ソフトセンサーやMSPC（多変量統計的プロセス管理）では、上位主成分よりむしろ下位の小さな分散方向が異常検知に重要な情報を持つため、この精度劣化はモデルの信頼性に直結します。

 

まとめ〜scikit-learn が SVD を使う理由〜

これらすべての理由を踏まえると、scikit-learn の PCA・PLS が内部実装で特異値分解を使う理由は、

という、まさに「実装する側にとって優位な点ばかり」の手法だからです。

特に化学プラントのように多重共線性が強いプロセスデータを扱う場合、この数値安定性の差は無視できません。

教科書では「共分散行列の固有値分解」と教わりますが、実装では \(X\) を直接SVDで分解する方が、結果として得られる主成分の方向もより正確になります。

 

参考文献

１．スモールデータ解析と機械学習

本記事でも引き続き、こちらの書籍を参考にしました。

機械学習を言葉多めでわかりやすく解説してくれる書籍で、線形代数のTipsも盛りだくさんです。

 

２．意味がわかる線形代数

初心者向けの線形代数の参考書です。

「意味」を理解することに重きを置いているため、統計学や機械学習へ応用する際にとても重宝します。

統計学や機械学習で詰まったらこの本に立ち戻ると良いと思います。