線形代数学 I / 基礎 / II

第39回講義：特異値分解の基礎

1. 講義情報と予習ガイド

講義回: 第39回
関連項目: 特異値分解（SVD）、線形代数の応用
予習すべき内容: 行列の対角化、固有値と固有ベクトル、対称行列の性質

2. 学習目標

特異値分解（SVD）の基本概念と数学的定義を理解する
特異値と特異ベクトルの計算方法を習得する
特異値分解の幾何学的解釈を理解する
特異値分解のデータサイエンスにおける重要性を説明できるようになる
Python（NumPy）を用いて特異値分解を実装できるようになる

3. 基本概念

3.1 特異値分解とは

特異値分解（Singular Value Decomposition, SVD）は、任意の行列を三つの行列の積に分解する手法です。これは行列分解の一種であり、データサイエンスにおいて非常に重要なツールとなっています。

定義: \(m \times n\)の行列 \(A\) に対する特異値分解とは、以下の形式で表される分解のことです：

\(A = U \Sigma V^T\)

ここで、 - \(U\) は \(m \times m\) の直交行列（左特異ベクトル） - \(\Sigma\) は \(m \times n\) の対角行列（対角成分に特異値を持つ） - \(V^T\) は \(n \times n\) の直交行列の転置（右特異ベクトル）

特異値分解は任意の行列に対して存在し、この普遍性が特異値分解をデータサイエンスにおいて強力なツールにしています。

3.2 特異値とは

特異値は、行列の「重要度」や「スケール」を表す非負の実数です。

定義: 行列 \(A\) の特異値は、\(A^TA\) の固有値の平方根として定義されます。

すなわち、\(\sigma_i = \sqrt{\lambda_i(A^TA)}\) です。

特異値は常に非負であり、通常は大きい順に並べられます：\(\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r > 0 = \sigma_{r+1} = \cdots = \sigma_{\min(m,n)}\)

ここで、\(r\) は行列 \(A\) の階数（rank）です。

3.3 特異ベクトルとは

特異値分解では、二種類の特異ベクトルが登場します：左特異ベクトルと右特異ベクトルです。

定義: - 右特異ベクトル \(v_i\) は、\(A^TA\) の固有ベクトルとして定義されます。 - 左特異ベクトル \(u_i\) は、\(AA^T\) の固有ベクトル、または \(u_i = \frac{Av_i}{\sigma_i}\) として定義されます（\(\sigma_i \neq 0\) の場合）。

行列 \(U\) の列は左特異ベクトル、行列 \(V\) の列は右特異ベクトルです。これらの特異ベクトルは、それぞれ直交基底を形成します。

4. 理論と手法

4.1 特異値分解の存在定理

特異値分解は任意の行列に対して存在します。これは特異値分解が持つ重要な性質の一つです。

定理: 任意の \(m \times n\) 行列 \(A\) に対して、特異値分解 \(A = U \Sigma V^T\) が存在します。

この普遍性によって、様々な形状や性質を持つ行列に対して特異値分解を適用することができます。

4.2 特異値分解の計算手順

特異値分解を計算するための一般的な手順は以下の通りです：

行列 \(A^TA\) を計算する
\(A^TA\) の固有値と固有ベクトルを求める
固有値 \(\lambda_i\) が \(A^TA\) の固有値
固有ベクトル \(v_i\) が右特異ベクトル
特異値 \(\sigma_i = \sqrt{\lambda_i}\) を計算する
左特異ベクトル \(u_i\) を \(u_i = \frac{Av_i}{\sigma_i}\) として計算する（\(\sigma_i \neq 0\) の場合）
\(U\), \(\Sigma\), \(V^T\) を構成する

4.3 特異値と固有値の関係

特異値と固有値には重要な関係があります：

\(A^TA\) の固有値は \(AA^T\) の非ゼロ固有値と一致する
行列 \(A\) の特異値は \(A^TA\) の固有値の平方根
行列 \(A\) のランク \(r\) は、非ゼロ特異値の数に等しい

4.4 計算例

具体的な計算例を通して特異値分解の手順を見てみましょう。

例題: 以下の \(2 \times 2\) 行列 \(A\) の特異値分解を求めなさい。

\[A = \begin{bmatrix} 3 & 2 \\ 2 & 2 \end{bmatrix}\]

解答:

ステップ 1: \(A^TA\) を計算します。

\[A^TA = \begin{bmatrix} 3 & 2 \\ 2 & 2 \end{bmatrix}^T \begin{bmatrix} 3 & 2 \\ 2 & 2 \end{bmatrix} = \begin{bmatrix} 3 & 2 \\ 2 & 2 \end{bmatrix} \begin{bmatrix} 3 & 2 \\ 2 & 2 \end{bmatrix} = \begin{bmatrix} 13 & 10 \\ 10 & 8 \end{bmatrix}\]

ステップ 2: \(A^TA\) の固有値と固有ベクトルを求めます。

特性多項式は：

\[\det(A^TA - \lambda I) = \det\begin{bmatrix} 13-\lambda & 10 \\ 10 & 8-\lambda \end{bmatrix} = (13-\lambda)(8-\lambda) - 10 \cdot 10 = \lambda^2 - 21\lambda + 4 = 0\]

この二次方程式を解くと：

\[\lambda_1 = \frac{21 + \sqrt{21^2 - 16}}{2} \approx 20.39, \quad \lambda_2 = \frac{21 - \sqrt{21^2 - 16}}{2} \approx 0.61\]

それぞれの固有値に対応する固有ベクトルを求めます：

\(\lambda_1 = 20.39\) に対して： \(\(\begin{bmatrix} 13-20.39 & 10 \\ 10 & 8-20.39 \end{bmatrix} \begin{bmatrix} v_{11} \\ v_{12} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}\)\)

これを解くと \(v_1 \approx \begin{bmatrix} 0.8 \\ 0.6 \end{bmatrix}\) （正規化後）

\(\lambda_2 = 0.61\) に対して： \(\(\begin{bmatrix} 13-0.61 & 10 \\ 10 & 8-0.61 \end{bmatrix} \begin{bmatrix} v_{21} \\ v_{22} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}\)\)

これを解くと \(v_2 \approx \begin{bmatrix} -0.6 \\ 0.8 \end{bmatrix}\) （正規化後）

ステップ 3: 特異値を計算します。

\[\sigma_1 = \sqrt{\lambda_1} = \sqrt{20.39} \approx 4.52, \quad \sigma_2 = \sqrt{\lambda_2} = \sqrt{0.61} \approx 0.78\]

ステップ 4: 左特異ベクトルを計算します。

\[u_1 = \frac{Av_1}{\sigma_1} = \frac{1}{4.52} \begin{bmatrix} 3 & 2 \\ 2 & 2 \end{bmatrix} \begin{bmatrix} 0.8 \\ 0.6 \end{bmatrix} \approx \begin{bmatrix} 0.71 \\ 0.71 \end{bmatrix}\]

\[u_2 = \frac{Av_2}{\sigma_2} = \frac{1}{0.78} \begin{bmatrix} 3 & 2 \\ 2 & 2 \end{bmatrix} \begin{bmatrix} -0.6 \\ 0.8 \end{bmatrix} \approx \begin{bmatrix} -0.71 \\ 0.71 \end{bmatrix}\]

ステップ 5: 特異値分解を構成します。

\[A = U\Sigma V^T = \begin{bmatrix} 0.71 & -0.71 \\ 0.71 & 0.71 \end{bmatrix} \begin{bmatrix} 4.52 & 0 \\ 0 & 0.78 \end{bmatrix} \begin{bmatrix} 0.8 & -0.6 \\ 0.6 & 0.8 \end{bmatrix}\]

以上が行列 \(A\) の特異値分解です。

4.5 特異値分解の幾何学的解釈

特異値分解は線形変換の幾何学的な意味を明らかにします。行列 \(A\) による線形変換は以下の3つのステップに分解できます：

直交変換 \(V^T\): 元の空間での回転
対角行列 \(\Sigma\) による拡大・縮小: 各軸方向への伸縮
直交変換 \(U\): 変換後の空間での回転

これは、任意の線形変換が「回転→伸縮→回転」という3つの基本的な操作の組み合わせで表現できることを意味します。

特異値分解の幾何学的な解釈では、単位球（または超球）の変換を考えるとわかりやすくなります：

\(V^T\) によって単位球を回転させる
\(\Sigma\) によって各軸方向に伸縮させ、楕円（または超楕円）にする
\(U\) によって楕円を回転させる

最終的に得られる楕円の半軸の長さが特異値 \(\sigma_i\) に対応します。

5. Pythonによる実装と可視化

特異値分解はNumPyを使用して簡単に計算することができます。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

# 行列の定義
A = np.array([[3, 2], 
              [2, 2]])

# 特異値分解の計算
U, sigma, VT = np.linalg.svd(A)

print("左特異ベクトル行列 U:")
print(U)
print("\n特異値:")
print(sigma)
print("\n右特異ベクトル行列 V^T:")
print(VT)

# 特異値分解の検証
Sigma = np.zeros((A.shape[0], A.shape[1]))
np.fill_diagonal(Sigma, sigma)
A_reconstructed = U @ Sigma @ VT
print("\n再構成された行列 A:")
print(A_reconstructed)
print("\n元の行列 A:")
print(A)
print("\n再構成誤差:")
print(np.linalg.norm(A - A_reconstructed))

# 特異値分解の幾何学的解釈の可視化
def plot_transformation(A, U, sigma, VT):
    # 座標の範囲を設定
    x = np.linspace(-1, 1, 100)
    y = np.linspace(-1, 1, 100)
    X, Y = np.meshgrid(x, y)

    # 単位円上の点を生成
    theta = np.linspace(0, 2*np.pi, 100)
    circle_x = np.cos(theta)
    circle_y = np.sin(theta)
    circle_points = np.vstack([circle_x, circle_y])

    # 変換された楕円の計算
    transformed_points = A @ circle_points

    # プロットの設定
    fig, axs = plt.subplots(1, 4, figsize=(20, 5))

    # 元の単位円
    axs[0].plot(circle_x, circle_y, 'b-')
    axs[0].set_xlim(-3, 3)
    axs[0].set_ylim(-3, 3)
    axs[0].set_aspect('equal')
    axs[0].grid(True)
    axs[0].set_title('元の単位円')

    # V^Tによる回転
    rotated_points = VT @ circle_points
    axs[1].plot(rotated_points[0], rotated_points[1], 'g-')
    axs[1].set_xlim(-3, 3)
    axs[1].set_ylim(-3, 3)
    axs[1].set_aspect('equal')
    axs[1].grid(True)
    axs[1].set_title('V^Tによる回転')

    # Σによる拡大・縮小
    scaled_points = np.diag(sigma) @ rotated_points
    axs[2].plot(scaled_points[0], scaled_points[1], 'r-')
    axs[2].set_xlim(-5, 5)
    axs[2].set_ylim(-5, 5)
    axs[2].set_aspect('equal')
    axs[2].grid(True)
    axs[2].set_title('Σによる拡大・縮小')

    # Uによる回転（最終結果）
    axs[3].plot(transformed_points[0], transformed_points[1], 'm-')
    axs[3].set_xlim(-5, 5)
    axs[3].set_ylim(-5, 5)
    axs[3].set_aspect('equal')
    axs[3].grid(True)
    axs[3].set_title('Aによる変換（最終結果）')

    plt.tight_layout()
    plt.show()

# 行列変換の可視化
plot_transformation(A, U, sigma, VT)

この例では、2×2行列の特異値分解を計算し、その結果を表示しています。また、特異値分解の幾何学的解釈を可視化するために、単位円がどのように変換されるかを段階的に示しています。

6. データサイエンスにおける特異値分解の応用

特異値分解はデータサイエンス、特に健康データ分析において多くの応用があります：

次元削減: 大きな特異値のみを使用してデータの近似表現を得ることができます
ノイズ除去: 小さな特異値を0にすることでノイズを取り除くことができます
画像圧縮: 画像データを少ない情報量で近似できます
推薦システム: ユーザー-アイテム評価行列の潜在構造を発見するのに役立ちます
健康データ解析: 医療画像や生体信号のパターン認識に使われます

例えば、医療画像の場合、SVDを用いることで画像の特徴を効率的に抽出し、ノイズを除去することが可能です。また、大量の健康データから潜在的なパターンを発見するためにも有効です。

# 健康データ分析における特異値分解の応用例
# 例：MRIデータの圧縮とノイズ除去

import numpy as np
import matplotlib.pyplot as plt
from skimage import color
from skimage import data

# サンプル画像の読み込み（実際のMRIデータの代わりに）
image = color.rgb2gray(data.astronaut())
plt.figure(figsize=(5, 5))
plt.imshow(image, cmap='gray')
plt.title('元の画像')
plt.axis('off')
plt.show()

# 特異値分解
U, sigma, VT = np.linalg.svd(image, full_matrices=False)

# 特異値のプロット
plt.figure(figsize=(10, 4))
plt.semilogy(sigma)
plt.title('特異値の大きさ')
plt.xlabel('インデックス')
plt.ylabel('特異値 (対数スケール)')
plt.grid(True)
plt.show()

# 異なる数の特異値を使用した画像の再構成
ranks = [5, 20, 50, 100]
fig, axes = plt.subplots(1, len(ranks), figsize=(20, 5))

for i, r in enumerate(ranks):
    # r個の特異値を使用して画像を再構成
    reconstructed = U[:, :r] @ np.diag(sigma[:r]) @ VT[:r, :]

    # 再構成された画像の表示
    axes[i].imshow(reconstructed, cmap='gray')
    axes[i].set_title(f'上位 {r} 特異値での再構成')
    axes[i].axis('off')

    # 情報圧縮率の計算
    original_size = image.shape[0] * image.shape[1]
    compressed_size = r * (image.shape[0] + image.shape[1] + 1)
    compression_ratio = compressed_size / original_size * 100

    print(f'上位 {r} 特異値を使用: 圧縮率 = {compression_ratio:.2f}%, RMSE = {np.sqrt(np.mean((image - reconstructed)**2)):.4f}')

plt.tight_layout()
plt.show()

このコード例では、画像の特異値分解を計算し、異なる数の特異値を使用して画像を再構成する方法を示しています。これは、データの次元削減や圧縮の基本的な例です。

7. 演習問題

基本問題

次の行列の特異値分解を求めなさい。 \(\(A = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}\)\)
特異値分解 \(A = U\Sigma V^T\) において、\(U\) と \(V\) が直交行列であることを証明しなさい。
行列 \(A\) が正方行列かつ対称行列である場合、その特異値分解と固有値分解の関係を説明しなさい。
次の行列の特異値を求め、ランクを決定しなさい。 \(\(B = \begin{bmatrix} 2 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 0 \end{bmatrix}\)\)
行列 \(A\) の特異値が \(\sigma_1, \sigma_2, \ldots, \sigma_n\) である場合、\(A^TA\) の固有値は何か？

応用問題

以下の \(3 \times 4\) 行列の特異値分解を計算しなさい。 \(\(C = \begin{bmatrix} 4 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & 2 & 0 \end{bmatrix}\)\) 特異値、左特異ベクトル、右特異ベクトルをそれぞれ求め、\(C = U\Sigma V^T\) の形で表しなさい。また、行列の階数を特異値から決定しなさい。
以下の \(2 \times 3\) 行列の特異値分解を計算しなさい。 \(\(D = \begin{bmatrix} 2 & 1 & 0 \\ 1 & 2 & 1 \end{bmatrix}\)\) 特異値、左特異ベクトル、右特異ベクトルをそれぞれ求め、\(D = U\Sigma V^T\) の形で表しなさい。
特異値分解を用いて、以下の長方形行列の擬似逆行列（ムーア・ペンローズ擬似逆行列）を計算しなさい。 \(\(E = \begin{bmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{bmatrix}\)\) ヒント：特異値分解 \(E = U\Sigma V^T\) が与えられたとき、擬似逆行列 \(E^+\) は \(E^+ = V\Sigma^+U^T\) として与えられる。ここで \(\Sigma^+\) は \(\Sigma\) の非ゼロ要素の逆数を対応する位置に持つ行列である。
NumPyを使用して、ランダムな100×50の行列を生成し、その特異値分解を計算しなさい。上位10個の特異値のみを使用して元の行列を近似したときの再構成誤差を計算し、フロベニウスノルムの意味での相対誤差を求めなさい。
健康データ分析における実用的な問題：100人の患者から収集された10種類のバイオマーカー（血圧、血糖値、コレステロール値など）の時系列データがあるとします。このデータは100×10の行列として表現できます。このようなデータに対して特異値分解を適用することで、どのような洞察が得られるか説明しなさい。特に、主要な特異値とそれに対応する特異ベクトルの解釈に焦点を当て、患者グループの分類や潜在的な健康パターンの発見にどのように役立つか考察しなさい。

8. よくある質問と解答

Q1: 特異値分解と固有値分解の違いは何ですか？

A1: 固有値分解は正方行列にのみ適用可能であり、必ずしもすべての行列に対して存在するわけではありません。一方、特異値分解は任意の行列（正方、長方形を問わず）に適用可能で、常に存在します。また、固有値分解では \(A = PDP^{-1}\) という形式で分解しますが、特異値分解では \(A = U\Sigma V^T\) という形式で分解します。ここで、\(U\) と \(V\) は直交行列であり、\(\Sigma\) は対角行列です。

Q2: 特異値分解の計算コストはどのくらいですか？

A2: \(m \times n\) 行列の特異値分解の計算量は一般的に \(O(\min(mn^2, m^2n))\) です。大規模な行列に対しては、計算コストが高くなる可能性があります。しかし、部分的な特異値分解（上位k個の特異値と対応する特異ベクトルのみを計算）を使用することで、計算コストを削減できる場合があります。

Q3: 特異値がゼロに近い場合、どのように処理すべきですか？

A3: 特異値がゼロに非常に近い場合、数値計算の誤差によって問題が生じる可能性があります。このような場合、閾値を設定し、それ以下の特異値をゼロとして扱うことが一般的です。これは「打ち切り特異値分解（truncated SVD）」と呼ばれ、ノイズの除去や数値安定性の向上に役立ちます。

Q4: 特異値分解はどのようにして健康データの解析に役立ちますか？

A4: 健康データの解析において、特異値分解は以下のように役立ちます： - 大量の生体信号データ（脳波、心電図など）から重要なパターンを抽出 - 医療画像のノイズ除去と特徴抽出 - 患者-症状データの潜在構造の発見 - 遺伝子発現データの次元削減 - 時系列健康データの傾向分析

例えば、複数の患者から収集された時系列の健康指標データに特異値分解を適用することで、共通のパターンや異常を検出することができます。

Q5: 特異値分解と主成分分析（PCA）の関係は何ですか？

A5: 主成分分析（PCA）は、データの共分散行列に対する固有値分解を基にしています。中心化されたデータ行列 \(X\) に対して、\(X^TX/(n-1)\) が共分散行列になります。このデータ行列 \(X\) に直接特異値分解を適用すると、得られる右特異ベクトル \(V\) は、共分散行列の固有ベクトル（主成分）と一致します。また、特異値の二乗を \((n-1)\) で割ったものが、対応する主成分の分散になります。つまり、PCAは特異値分解の特殊なケースと見なすことができます。

Q6: 特異値分解におけるランク低減はどのような意味を持ちますか？

A6: ランク低減（rank reduction）とは、小さな特異値を0とし、対応する特異ベクトルを無視することで、元の行列の低ランク近似を得る手法です。これにより、データの「ノイズ」を除去し、「信号」部分を保持することができます。エッカート・ヤングの定理によれば、上位k個の特異値に基づく近似は、ランクkの行列の中で元の行列に最も近いものとなります（フロベニウスノルムの意味で）。これはデータ圧縮やノイズ除去に非常に有用です。

応用問題5の解答

特異値分解（SVD）は、100人の患者から収集された10種類のバイオマーカーの時系列データ（100×10行列）の分析において、以下のような重要な洞察を提供することができます。（実例

1. データの潜在構造の抽出

特異値分解を適用すると、データの潜在的な構造を特異値と特異ベクトルによって表現することができます。100×10の行列 \(X\) を特異値分解すると：

\[X = U\Sigma V^T\]

ここで： - \(U\) は100×100の直交行列で、列（左特異ベクトル）は患者の「潜在パターン」を表す - \(\Sigma\) は100×10の対角行列で、特異値を対角成分に持つ - \(V^T\) は10×10の直交行列で、行（右特異ベクトル）はバイオマーカーの「特徴パターン」を表す

2. 主要な変動パターンの特定

特異値の大きさは、対応する特異ベクトルの重要度を示します。最大の特異値 \(\sigma_1\) に対応する右特異ベクトル \(v_1\) は、バイオマーカー間の最も支配的な相関パターンを表します。例えば：

\(v_1\) の成分が血圧と心拍数に大きな値を持つ場合、これらの指標が強く相関していることを示唆します
\(v_1\) の成分が正負の値を持つ場合、それらのバイオマーカー間に逆相関があることを示します

同様に、左特異ベクトル \(u_1\) は、このパターンに最も強く反応する患者グループを特定します。

3. 患者の分類とサブグループの発見

左特異ベクトル \(U\) の列は、患者を異なる「潜在的健康状態」に分類するための基底を提供します。例えば：

特定の左特異ベクトル \(u_i\) における値に基づいて患者をクラスタリングすることで、類似した健康プロファイルを持つサブグループを特定できます
これらのサブグループは、異なる病態生理学的プロセスや治療反応パターンを表す可能性があります

実際には、各患者を上位2〜3個の左特異ベクトルの空間にプロットすることで、視覚的に患者クラスターを識別することができます。

4. 次元削減とノイズ除去

特異値の大きさが急速に減少する場合、上位k個の特異値と対応する特異ベクトルのみを使用してデータを近似することができます：

\[X_k = \sum_{i=1}^{k} \sigma_i u_i v_i^T\]

この低ランク近似により： - データの本質的な構造を保持しながら次元を削減できます - 測定ノイズや偶然の変動を除去できます - 小さい特異値に対応する成分は、多くの場合ノイズを表しています

例えば、上位3つの特異値で全変動の80%を説明できる場合、10次元のデータを3次元に効果的に削減できます。

5. バイオマーカーの重要度評価

右特異ベクトル \(V\) の行は、各バイオマーカーがデータの主要な変動パターンにどのように寄与しているかを示します：

複数の右特異ベクトルで大きな成分を持つバイオマーカーは、患者群のさまざまな健康側面に重要な役割を果たしています
特定の右特異ベクトルで大きな値を持つバイオマーカーのグループは、共変動する可能性が高く、共通の生理学的経路に関連している可能性があります

これにより、診断や治療モニタリングに最も有用なバイオマーカーの組み合わせを特定できます。

6. 時間的変化の検出と予測

時系列データにSVDを適用することで、時間とともに変化するパターンを特定できます：

主要な特異ベクトルの時間的変化をモニタリングすることで、疾患の進行や治療への反応を追跡できます
患者の時間的軌跡を特異ベクトル空間で分析することで、予後予測モデルを開発できる可能性があります

例えば、治療に良く反応する患者は特異ベクトル空間で類似した軌跡をたどる傾向があるかもしれません。