コンテンツにスキップ

QR 分解

QR 分解は、対称行列と非対称行列を含む一般の行列を対角化するための、最も効率的で広く使われている方法の一つです。QR アルゴリズムは、行列 AA を直交行列 QQ と上三角行列 RR の積に反復的に分解することで機能します。この分解を繰り返し適用し、行列を A′=RQA^{\prime} = RQ として再構成することで、行列は対角形または三角形に収束し、そこから固有値を取り出すことができます。固有ベクトルは QQ 行列の累積積から得られます。

数学的基礎

行列 AA に対して、QR 分解は次で与えられます。

A=QR A = QR

ここで、

  • QQ は直交行列(QTQ=IQ^T Q = I)、
  • RR は上三角行列です。

QR アルゴリズムは、この分解を反復的に適用して AA を対角形または三角形に収束させます。

  1. A0=AA_0 = A から始めます。
  2. 各反復 kk について、
    • QR 分解を計算します:Ak=QkRkA_k = Q_k R_k。
    • 行列を再構成します:Ak+1=RkQkA_{k+1} = R_k Q_k。
  3. AkA_k が対角行列または三角行列に収束するまで繰り返します。

AA の固有値は最終的な行列 AkA_k の対角成分として得られ、固有ベクトルはすべての QkQ_k 行列の積から得られます。

アルゴリズム

  1. 初期化:

    • 行列 A0=AA_0 = A から始めます。
  2. QR 分解:

    • AkA_k を QkQ_k と RkR_k に分解します。 Ak=QkRk A_k = Q_k R_k
  3. 再構成:

    • 行列 Ak+1A_{k+1} を次のように再構成します。 Ak+1=RkQk A_{k+1} = R_k Q_k
  4. 変換の累積:

    • 固有ベクトル行列 PP を次のように更新します。 Pk+1=PkQk P_{k+1} = P_k Q_k
    • P0=IP_0 = I(単位行列)と初期化します。
  5. 収束の確認:

    • AkA_k が十分に対角形または三角形になる(すなわち、非対角要素が指定した許容誤差を下回る)まで、この過程を繰り返します。
  6. 固有値と固有ベクトルの抽出:

    • 固有値は最終的な AkA_k の対角要素です。
    • 固有ベクトルは最終的な PkP_k の各列です。

例

例として、単精度演算を用いて実対称行列の QR 分解を行います。

ステップ 1:初期化

行列 AA から始めます。

A=(4.000000−1.0000003.000000 −1.0000003.000000−1.000000 3.000000−1.0000005.000000). A = \begin{pmatrix} 4.000000 & -1.000000 & 3.000000 \\\ -1.000000 & 3.000000 & -1.000000 \\\ 3.000000 & -1.000000 & 5.000000 \end{pmatrix}.

ステップ 2:1 回目の QR 反復

ステップ 2.1:QQ と RR の計算

グラム・シュミットの直交化法を用いて AA の QR 分解を行います。

  • QQ の第 1 列: AA の第 1 列を規格化します。

    a1=(4.000000 −1.000000 3.000000),∥a1∥=42+(−1)2+32=26≈5.099020. \mathbf{a}_1 = \begin{pmatrix} 4.000000 \\\ -1.000000 \\\ 3.000000 \end{pmatrix}, \quad \|\mathbf{a}_1\| = \sqrt{4^2 + (-1)^2 + 3^2} = \sqrt{26} \approx 5.099020.

    したがって、

    q1=15.099020(4.000000 −1.000000 3.000000)≈(0.784465 −0.196116 0.588349). \mathbf{q}_1 = \frac{1}{5.099020} \begin{pmatrix} 4.000000 \\\ -1.000000 \\\ 3.000000 \end{pmatrix} \approx \begin{pmatrix} 0.784465 \\\ -0.196116 \\\ 0.588349 \end{pmatrix}.
  • QQ の第 2 列: AA の第 2 列を q1\mathbf{q}_1 に対して直交化します。

    a2=(−1.000000 3.000000 −1.000000),a2⋅q1≈−1.960784. \mathbf{a}_2 = \begin{pmatrix} -1.000000 \\\ 3.000000 \\\ -1.000000 \end{pmatrix}, \quad \mathbf{a}_2 \cdot \mathbf{q}_1 \approx -1.960784.

    v2\mathbf{v}_2 を計算します。

    v2=a2−(a2⋅q1)q1≈(−1.000000 3.000000 −1.000000)−(−1.960784)(0.784465 −0.196116 0.588349). \mathbf{v}_2 = \mathbf{a}_2 - (\mathbf{a}_2 \cdot \mathbf{q}_1) \mathbf{q}_1 \approx \begin{pmatrix} -1.000000 \\\ 3.000000 \\\ -1.000000 \end{pmatrix} - (-1.960784) \begin{pmatrix} 0.784465 \\\ -0.196116 \\\ 0.588349 \end{pmatrix}.

    v2≈(−1.000000+1.538462 3.000000−0.384615 −1.000000+1.153846)=(0.538462 2.615385 0.153846). \mathbf{v}_2 \approx \begin{pmatrix} -1.000000 + 1.538462 \\\ 3.000000 - 0.384615 \\\ -1.000000 + 1.153846 \end{pmatrix} = \begin{pmatrix} 0.538462 \\\ 2.615385 \\\ 0.153846 \end{pmatrix}.

    v2\mathbf{v}_2 を規格化します。

    ∥v2∥=0.5384622+2.6153852+0.1538462≈2.672612. \|\mathbf{v}_2\| = \sqrt{0.538462^2 + 2.615385^2 + 0.153846^2} \approx 2.672612.

    したがって、

    q2≈(0.2014560.9785930.057553). \mathbf{q}_2 \approx \begin{pmatrix} 0.201456 \\ 0.978593 \\ 0.057553 \end{pmatrix}.
  • QQ の第 3 列: AA の第 3 列を q1\mathbf{q}_1 と q2\mathbf{q}_2 に対して直交化します。

    a3=(3.000000 −1.000000 5.000000),a3⋅q1≈5.882353, \mathbf{a}_3 = \begin{pmatrix} 3.000000 \\\ -1.000000 \\\ 5.000000 \end{pmatrix}, \quad \mathbf{a}_3 \cdot \mathbf{q}_1 \approx 5.882353,

    a3⋅q2≈0.000000. \mathbf{a}_3 \cdot \mathbf{q}_2 \approx 0.000000.

    v3\mathbf{v}_3 を計算します。

    v3=a3−(a3⋅q1)q1−(a3⋅q2)q2. \mathbf{v}_3 = \mathbf{a}_3 - (\mathbf{a}_3 \cdot \mathbf{q}_1) \mathbf{q}_1 - (\mathbf{a}_3 \cdot \mathbf{q}_2) \mathbf{q}_2.

    v3≈(3.000000 −1.000000 5.000000)−5.882353(0.784465 −0.196116 0.588349)−0.000000(0.201456 0.978593 0.057553). \mathbf{v}_3 \approx \begin{pmatrix} 3.000000 \\\ -1.000000 \\\ 5.000000 \end{pmatrix} - 5.882353 \begin{pmatrix} 0.784465 \\\ -0.196116 \\\ 0.588349 \end{pmatrix} - 0.000000 \begin{pmatrix} 0.201456 \\\ 0.978593 \\\ 0.057553 \end{pmatrix}.

    v3≈(3.000000−4.615385 −1.000000+1.153846 5.000000−3.461538)=(−1.615385 0.153846 1.538462). \mathbf{v}_3 \approx \begin{pmatrix} 3.000000 - 4.615385 \\\ -1.000000 + 1.153846 \\\ 5.000000 - 3.461538 \end{pmatrix} = \begin{pmatrix} -1.615385 \\\ 0.153846 \\\ 1.538462 \end{pmatrix}.

    v3\mathbf{v}_3 を規格化します。

    ∥v3∥=(−1.615385)2+0.1538462+1.5384622≈2.236068. \|\mathbf{v}_3\| = \sqrt{(-1.615385)^2 + 0.153846^2 + 1.538462^2} \approx 2.236068.

    したがって、

    q3≈(−0.722222 0.068783 0.688889). \mathbf{q}_3 \approx \begin{pmatrix} -0.722222 \\\ 0.068783 \\\ 0.688889 \end{pmatrix}.
  • QQ と RR の構築:

    Q=(0.7844650.201456−0.722222 −0.1961160.9785930.068783 0.5883490.0575530.688889), Q = \begin{pmatrix} 0.784465 & 0.201456 & -0.722222 \\\ -0.196116 & 0.978593 & 0.068783 \\\ 0.588349 & 0.057553 & 0.688889 \end{pmatrix},

    R=QTA≈(5.0990200.0000005.882353 02.6726120.000000 002.236068). R = Q^T A \approx \begin{pmatrix} 5.099020 & 0.000000 & 5.882353 \\\ 0 & 2.672612 & 0.000000 \\\ 0 & 0 & 2.236068 \end{pmatrix}.

ステップ 2.2:AA の更新

A=RQA = RQ を計算します。

A=RQ≈(6.561553−0.7592570.000000 −0.7592573.000000−0.650791 0.000000−0.6507912.438447). A = RQ \approx \begin{pmatrix} 6.561553 & -0.759257 & 0.000000 \\\ -0.759257 & 3.000000 & -0.650791 \\\ 0.000000 & -0.650791 & 2.438447 \end{pmatrix}.

ステップ 3:2 回目の QR 反復

ステップ 3.1:QQ と RR の計算

更新された AA の QR 分解を行います。

  • QQ の第 1 列: AA の第 1 列を規格化します。

    a1=(6.561553 −0.759257 0.000000),∥a1∥≈6.617647. \mathbf{a}_1 = \begin{pmatrix} 6.561553 \\\ -0.759257 \\\ 0.000000 \end{pmatrix}, \quad \|\mathbf{a}_1\| \approx 6.617647.

    したがって、

    q1≈(0.990000 −0.114706 0.000000). \mathbf{q}_1 \approx \begin{pmatrix} 0.990000 \\\ -0.114706 \\\ 0.000000 \end{pmatrix}.
  • QQ の第 2 列: AA の第 2 列を q1\mathbf{q}_1 に対して直交化します。

    a2=(−0.759257 3.000000 −0.650791),a2⋅q1≈−1.139000. \mathbf{a}_2 = \begin{pmatrix} -0.759257 \\\ 3.000000 \\\ -0.650791 \end{pmatrix}, \quad \mathbf{a}_2 \cdot \mathbf{q}_1 \approx -1.139000.

    v2\mathbf{v}_2 を計算します。

    v2=a2−(a2⋅q1)q1≈(−0.759257 3.000000 −0.650791)−(−1.139000)(0.990000 −0.114706 0.000000). \mathbf{v}_2 = \mathbf{a}_2 - (\mathbf{a}_2 \cdot \mathbf{q}_1) \mathbf{q}_1 \approx \begin{pmatrix} -0.759257 \\\ 3.000000 \\\ -0.650791 \end{pmatrix} - (-1.139000) \begin{pmatrix} 0.990000 \\\ -0.114706 \\\ 0.000000 \end{pmatrix}.

    v2≈(−0.759257+1.127610 3.000000−0.130000 −0.650791+0.000000)=(0.368353 2.870000 −0.650791). \mathbf{v}_2 \approx \begin{pmatrix} -0.759257 + 1.127610 \\\ 3.000000 - 0.130000 \\\ -0.650791 + 0.000000 \end{pmatrix} = \begin{pmatrix} 0.368353 \\\ 2.870000 \\\ -0.650791 \end{pmatrix}.

    v2\mathbf{v}_2 を規格化します。

    ∥v2∥≈2.939000. \|\mathbf{v}_2\| \approx 2.939000.

    したがって、

    q2≈(0.125000 0.974000 −0.221000). \mathbf{q}_2 \approx \begin{pmatrix} 0.125000 \\\ 0.974000 \\\ -0.221000 \end{pmatrix}.
  • QQ の第 3 列: AA の第 3 列を q1\mathbf{q}_1 と q2\mathbf{q}_2 に対して直交化します。

    a3=(0.000000 −0.650791 2.438447),a3⋅q1≈0.000000,a3⋅q2≈−0.624695. \mathbf{a}_3 = \begin{pmatrix} 0.000000 \\\ -0.650791 \\\ 2.438447 \end{pmatrix}, \quad \mathbf{a}_3 \cdot \mathbf{q}_1 \approx 0.000000, \quad \mathbf{a}_3 \cdot \mathbf{q}_2 \approx -0.624695.

    v3\mathbf{v}_3 を計算します。

    v3=a3−(a3⋅q1)q1−(a3⋅q2)q2. \mathbf{v}_3 = \mathbf{a}_3 - (\mathbf{a}_3 \cdot \mathbf{q}_1) \mathbf{q}_1 - (\mathbf{a}_3 \cdot \mathbf{q}_2) \mathbf{q}_2.

    v3≈(0.000000 −0.650791 2.438447)−0.000000(0.990000 −0.114706 0.000000)−(−0.624695)(0.125000 0.974000 −0.221000). \mathbf{v}_3 \approx \begin{pmatrix} 0.000000 \\\ -0.650791 \\\ 2.438447 \end{pmatrix} - 0.000000 \begin{pmatrix} 0.990000 \\\ -0.114706 \\\ 0.000000 \end{pmatrix} - (-0.624695) \begin{pmatrix} 0.125000 \\\ 0.974000 \\\ -0.221000 \end{pmatrix}.

    v3≈(0.000000+0.078087 −0.650791−0.608000 2.438447+0.138000)=(0.078087 −1.258791 2.576447). \mathbf{v}_3 \approx \begin{pmatrix} 0.000000 + 0.078087 \\\ -0.650791 - 0.608000 \\\ 2.438447 + 0.138000 \end{pmatrix} = \begin{pmatrix} 0.078087 \\\ -1.258791 \\\ 2.576447 \end{pmatrix}.

    v3\mathbf{v}_3 を規格化します。

    ∥v3∥≈2.828427. \|\mathbf{v}_3\| \approx 2.828427.

    したがって、

    q3≈(0.027600 −0.445000 0.911000). \mathbf{q}_3 \approx \begin{pmatrix} 0.027600 \\\ -0.445000 \\\ 0.911000 \end{pmatrix}.
  • QQ と RR の構築:

    Q=(0.9900000.1250000.027600 −0.1147060.974000−0.445000 0.000000−0.2210000.911000), Q = \begin{pmatrix} 0.990000 & 0.125000 & 0.027600 \\\ -0.114706 & 0.974000 & -0.445000 \\\ 0.000000 & -0.221000 & 0.911000 \end{pmatrix},

    R=QTA≈(6.6176470.0000000.000000 02.9390000.000000 002.828427). R = Q^T A \approx \begin{pmatrix} 6.617647 & 0.000000 & 0.000000 \\\ 0 & 2.939000 & 0.000000 \\\ 0 & 0 & 2.828427 \end{pmatrix}.

ステップ 3.2:AA の更新

A=RQA = RQ を計算します。

A=RQ≈(6.6176470.0000000.000000 0.0000002.9390000.000000 0.0000000.0000002.828427). A = RQ \approx \begin{pmatrix} 6.617647 & 0.000000 & 0.000000 \\\ 0.000000 & 2.939000 & 0.000000 \\\ 0.000000 & 0.000000 & 2.828427 \end{pmatrix}.

ステップ 4:収束の確認

すべての非対角要素が許容誤差 10−610^{-6} を下回っているかを確認します。この場合、非対角要素はすでにゼロなので、行列 AA は対角化されています。

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

A≈(6.6176470.0000000.000000 0.0000002.9390000.000000 0.0000000.0000002.828427). A \approx \begin{pmatrix} 6.617647 & 0.000000 & 0.000000 \\\ 0.000000 & 2.939000 & 0.000000 \\\ 0.000000 & 0.000000 & 2.828427 \end{pmatrix}.

最終結果

単精度の QR 分解を用いて計算した AA の固有値は次のとおりです。

λ1≈6.617647,λ2≈2.939000,λ3≈2.828427. \lambda_1 \approx 6.617647, \quad \lambda_2 \approx 2.939000, \quad \lambda_3 \approx 2.828427.

要点

単精度の QR 分解法で得られる固有値は真の値に近いものの、丸め誤差のためにわずかに異なることがあります。より高い精度を得るには、倍精度演算を用いるか、より厳しい収束条件のもとで反復回数を増やすことをお勧めします。

利点

  • 効率性:QR アルゴリズムは大きな行列に対して非常に効率的です。
  • 汎用性:対称行列と非対称行列の両方に使えます。
  • 安定性:このアルゴリズムは数値的に安定で頑健です。

制限

  • 計算コスト:非常に大きな行列では、QR 分解のステップの計算コストが高くなることがあります。
  • 非対称行列に対する収束の遅さ:非対称行列では、多くの反復が必要になることがあります。