コンテンツにスキップ

Jacobi 回転

Jacobi 回転法は、対称行列を対角化するための古典的な反復アルゴリズムです。小規模から中規模の行列に特に適しており、その単純さと頑健さで知られています。この方法は、一連の直交変換(回転)によって非対角要素を系統的に消去していき、最終的に元の行列の固有値を要素とする対角行列に収束させるものです。

Jacobi 法は、行列の最大の非対角要素を対象とし、それをゼロにする回転を適用します。この過程を、すべての非対角要素が十分小さくなるまで反復的に繰り返すと、対角行列が得られます。元の行列の固有値はその対角成分として得られ、固有ベクトルは過程中に適用したすべての回転行列の積から得られます。

原理

対称行列 AA に対して、次を満たす直交行列 PP を求めることが目標です。

D=PTAP D = P^T A P

ここで DD は AA の固有値を含む対角行列であり、PP の各列は対応する固有ベクトルです。

Jacobi 法は、AA に一連の直交変換(回転)を適用することでこれを実現します。各回転は特定の非対角要素 AijA_{ij} を対象とし、それをゼロにします。

回転行列

Jacobi 回転行列 RR は、単位行列と 4 つの要素だけが異なる直交行列です。

R=(1⋯0⋯0⋯0 ⋮⋱⋮⋮⋮ 0⋯cos⁡θ⋯−sin⁡θ⋯0 ⋮⋮⋱⋮⋮ 0⋯sin⁡θ⋯cos⁡θ⋯0 ⋮⋮⋮⋱⋮ 0⋯0⋯0⋯1) R = \begin{pmatrix} 1 & \cdots & 0 & \cdots & 0 & \cdots & 0 \\\ \vdots & \ddots & \vdots & & \vdots & & \vdots \\\ 0 & \cdots & \cos \theta & \cdots & -\sin \theta & \cdots & 0 \\\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\\ 0 & \cdots & \sin \theta & \cdots & \cos \theta & \cdots & 0 \\\ \vdots & & \vdots & & \vdots & \ddots & \vdots \\\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 1 \end{pmatrix}

ここで、cos⁡θ\cos \theta と sin⁡θ\sin \theta は第 ii 行・第 jj 行と第 ii 列・第 jj 列の交点に配置されます。角度 θ\theta は、非対角要素 AijA_{ij} がゼロになるように選ばれます。

アルゴリズム

  1. 最大の非対角要素を特定する:

    • 行列 AA の中で(絶対値が)最大の非対角要素 AijA_{ij} を見つけます。
  2. 回転角 θ\theta を計算する:

    • 角度 θ\theta は次を満たすように選びます。 tan⁡(2θ)=2AijAii−Ajj \tan(2\theta) = \frac{2A_{ij}}{A_{ii} - A_{jj}}
    • これから cos⁡θ\cos \theta と sin⁡θ\sin \theta を計算します。
  3. 回転行列 RR を構築する:

    • cos⁡θ\cos \theta と sin⁡θ\sin \theta を用いて回転行列 RR を構築します。
  4. 回転を適用する:

    • 行列 AA を次のように更新します。 A′=RTAR A^{\prime} = R^T A R
    • この変換により AijA_{ij} と AjiA_{ji} がゼロになります。
  5. 変換を累積する:

    • 固有ベクトル行列 PP を次のように更新します。 P′=PR P^{\prime} = P R
    • これにより回転が累積され、最終的な固有ベクトル行列が得られます。
  6. 収束するまで繰り返す:

    • すべての非対角要素が指定した許容誤差 ϵ \epsilon より小さくなるまで、この過程を繰り返します。

例

対称行列 AA を考えます。

A=(412 131 215) A = \begin{pmatrix} 4 & 1 & 2 \\\ 1 & 3 & 1 \\\ 2 & 1 & 5 \end{pmatrix}

ステップ 1:初期化

行列 AA と、回転を累積するための単位行列 PP から始めます。

A=(4.0000001.0000002.000000 1.0000003.0000001.000000 2.0000001.0000005.000000),P=(1.0000000.0000000.000000 0.0000001.0000000.000000 0.0000000.0000000.000000). A = \begin{pmatrix} 4.000000 & 1.000000 & 2.000000 \\\ 1.000000 & 3.000000 & 1.000000 \\\ 2.000000 & 1.000000 & 5.000000 \end{pmatrix}, \quad P = \begin{pmatrix} 1.000000 & 0.000000 & 0.000000 \\\ 0.000000 & 1.000000 & 0.000000 \\\ 0.000000 & 0.000000 & 0.000000 \end{pmatrix}.

ステップ 2:最大の非対角要素を見つける

AA の最大の非対角要素は A13=2.000000A_{13} = 2.000000 です。

ステップ 3:回転パラメータを計算する

(1,3)(1, 3) 平面内の回転について、

  • 角度 θ\theta を計算します。

    θ=12arctan⁡(2A13A11−A33). \theta = \frac{1}{2} \arctan\left(\frac{2A_{13}}{A_{11} - A_{33}}\right).

    値を代入すると、

    θ=12arctan⁡(2⋅2.0000004.000000−5.000000)=12arctan⁡(−4.000000). \theta = \frac{1}{2} \arctan\left(\frac{2 \cdot 2.000000}{4.000000 - 5.000000}\right) = \frac{1}{2} \arctan(-4.000000).

    単精度演算を用いると、

    θ≈−0.674741 radians. \theta \approx -0.674741 \, \text{radians}.
  • c=cos⁡(θ)c = \cos(\theta) と s=sin⁡(θ)s = \sin(\theta) を計算します。

    c≈0.780869,s≈−0.624695. c \approx 0.780869, \quad s \approx -0.624695.

ステップ 4:回転を適用する

回転行列 JJ を構築します。

J=(c0s 010 −s0c)≈(0.7808690.000000−0.624695 0.0000001.0000000.000000 0.6246950.0000000.780869). J = \begin{pmatrix} c & 0 & s \\\ 0 & 1 & 0 \\\ -s & 0 & c \end{pmatrix} \approx \begin{pmatrix} 0.780869 & 0.000000 & -0.624695 \\\ 0.000000 & 1.000000 & 0.000000 \\\ 0.624695 & 0.000000 & 0.780869 \end{pmatrix}.

AA と PP を更新します。

A=JTAJ,P=PJ. A = J^T A J, \quad P = P J.

回転後は次のようになります。

A≈(5.5615530.7808690.000000 0.7808693.0000000.624695 0.0000000.6246952.438447), A \approx \begin{pmatrix} 5.561553 & 0.780869 & 0.000000 \\\ 0.780869 & 3.000000 & 0.624695 \\\ 0.000000 & 0.624695 & 2.438447 \end{pmatrix},

P≈(0.7808690.000000−0.624695 0.0000001.0000000.000000 0.6246950.0000000.780869). P \approx \begin{pmatrix} 0.780869 & 0.000000 & -0.624695 \\\ 0.000000 & 1.000000 & 0.000000 \\\ 0.624695 & 0.000000 & 0.780869 \end{pmatrix}.

ステップ 5:他の非対角要素について繰り返す

すべての非対角要素が十分小さく(例えば許容誤差 10−610^{-6} 未満に)なるまで、次に大きい非対角要素についてこの過程を繰り返します。

ステップ 6:最終的な対角化行列

収束後、対角化された行列 AA は次のようになります。

A≈(6.0000000.0000000.000000 0.0000003.0000000.000000 0.0000000.0000003.000000). A \approx \begin{pmatrix} 6.000000 & 0.000000 & 0.000000 \\\ 0.000000 & 3.000000 & 0.000000 \\\ 0.000000 & 0.000000 & 3.000000 \end{pmatrix}.

対応する固有ベクトル行列 PP は次のようになります。

P≈(0.7071070.000000−0.707107 0.0000001.0000000.000000 0.7071070.0000000.707107). P \approx \begin{pmatrix} 0.707107 & 0.000000 & -0.707107 \\\ 0.000000 & 1.000000 & 0.000000 \\\ 0.707107 & 0.000000 & 0.707107 \end{pmatrix}.

利点

  • 単純さ:アルゴリズムの実装が簡単です。
  • 頑健さ:対称行列に対しては収束が保証されています。
  • 精度:非常に高精度な固有値と固有ベクトルが得られます。

制限

  • 収束の遅さ:大きな行列では多くの反復が必要です。
  • 大きな行列に対する非効率性:非常に大きな行列やスパース行列には適していません。
  • 計算コスト:各回転で行列全体を更新する必要があり、大きな系ではコストが高くなることがあります。