Jacobi法による固有値分解の解説とプログラミングの実装

はじめに

固有値問題は、数値計算や機械学習の分野で頻繁に登場する。本記事では、行列の固有値と固有ベクトルを求めるための数値計算手法である Jacobi法 を解説する。Pythonのコードと数式を組み合わせながら、アルゴリズムの動作を詳しく説明する。

---

Jacobi法とは？

Jacobi法は、回転行列 を使って対象の行列を対角化し、その結果として固有値と固有ベクトルを得る方法である。対象となる行列 $A$ に対して、以下のように行列の対角化を繰り返し行う。：

$$
A' = Q^T A Q
$$

ここで、$Q$ は回転行列である。回転を繰り返すことで、非対角要素が十分小さくなるまで処理を行う。

---

Jacobi法のアルゴリズム

Jacobi法は以下の手順で進行する：

1. 行列 $A$ の最大非対角要素 $(i, j)$ を特定する。
2. 回転行列 $Q$ を次の式に基づいて構築する：
   • 回転角度 $\phi$ は次式で決定：
     $$
     \phi = \frac{1}{2} \arctan\left(\frac{2A_{ij}}{A_{ii} - A_{jj}}\right)
     $$
   • 回転行列 $Q$ の成分：

     \displaylines{
         Q_{ii} = \cos{\phi} \\
         Q_{ij} = - \sin{\phi} \\
         Q_{ji} = \sin{\phi} \\
         Q_{jj} = \cos{\phi}
     }
   

3. 行列の更新：
   $$
   A' = Q^T A Q
   $$
4. 非対角要素が収束条件を満たすまで繰り返し。

---

Pythonコードの解説

以下は、PythonによるJacobi法の実装である。

import numpy as np
import numpy.linalg as LA

def jacobi(A, eps=1e-10):
    n = len(A)
    A = A.astype(np.float64)
    eigenvectors = np.eye(n)
    while True:
        max_idx = None
        max_val = 0
        for i in range(n):
            for j in range(i+1, n):
                if abs(A[i, j]) > max_val:
                    max_val = abs(A[i, j])
                    max_idx = (i, j)
        if max_val < eps:
            break
        i, j = max_idx
        if A[i, i] != A[j, j]:
            phi = 0.5 * np.arctan2(2 * A[i, j], A[i, i] - A[j, j])
        else:
            phi = np.pi / 4

        Q = np.eye(n)
        Q[i, i] = np.cos(phi)
        Q[j, j] = np.cos(phi)
        Q[i, j] = -np.sin(phi)
        Q[j, i] = np.sin(phi)
        
        A = Q.T @ A @ Q
        eigenvectors = eigenvectors @ Q
        
    return np.diag(A), eigenvectors

---

テストケース

以下の行列を使って、Jacobi法の動作を確認する。

A = np.array([[2, 1], [1, 3]])

jacobi関数の実行

eigvals, eigvecs = jacobi(A)
print('Eigenvalues:\n', eigvals)
print('Eigenvectors:\n', eigvecs)

numpy.linalgを使った検証

import numpy.linalg as LA
print('Check with Numpy:')
print('Eigenvalues:\n', LA.eig(A)[0])
print('Eigenvectors:\n', LA.eig(A)[1])

実行結果

Eigenvalues:
 [3.61803399 1.38196601]
Eigenvectors:
 [[ 0.52573111 -0.85065081]
 [ 0.85065081  0.52573111]]

Check with NumPy:
Eigenvalues:
 [1.38196601 3.61803399]
Eigenvectors:
 [[-0.85065081 -0.52573111]
 [ 0.52573111 -0.85065081]]

---

検証結果

• Jacobi法で求めた固有値と固有ベクトルは、NumPyの関数を使って求めた結果と一致していることが確認できる
• 一見すると、固有ベクトルの符号が異なるように見えるが、固有ベクトルは定数倍されても同じベクトルとして扱われるため、正しい結果である。また、固有値の順番が逆になっているが、これは固有値の順番を定める規則を持たないためである

---

まとめ

Jacobi法は、行列の対角化を行うための数値計算手法で、回転行列を用いて非対角要素を小さくしていく方法である。本記事では、Pythonコードを通じてその仕組みを説明した。さらに深掘りしたい場合は、QR法やパワー法など他の固有値問題のアルゴリズムも調べてみると良いだろう。

---