M.Hiroi's Home Page

Julia Language Programming

お気楽 Julia プログラミング超入門


Copyright (C) 2019-2021 Makoto Hiroi
All rights reserved.

線形代数編

●QR 分解と QR 法

前回は固有値と固有ベクトルを求める簡単な方法として「累乗法」と「ヤコビ法」を紹介しました。今回は行列の固有値を求める方法の一つである「QR 法 (QR algorithm)」と、その基礎になる「QR 分解 (QR decomposition)」について取り上げます。

QR 法の基本的な考え方はシンプルですが、これをそのままプログラムするのは効率的ではありません。このため、対称行列であれば「三重対角行列」に、それ以外の行列は「ヘッセンベルグ行列」に相似変換してから QR 法を適用するのが主流のようです。本稿では「実対称行列」の固有値を求める QR 法の基本について簡単に説明します。

●QR 分解の基本

まずは最初に QR 分解を説明しましょう。QR 分解は行列 \(A\) を直交行列 \(Q\) と上三角行列 \(R\) に分解することです。

\( A = QR \)

基本的な考え方は簡単です。\(A\) の左側から適当な直交行列 \(Q_i\) をいくつか掛け算したら上三角行列 \(R\) になったとしましょう。これを式で示すと次のようになります。

\( Q_n \cdots Q_2 Q_1 A = R \)

直交行列の逆行列は転置行列なので、上の式は次のように変換できます。

\( A = {Q_1}^{\mathrm{T}} {Q_2}^{\mathrm{T}} \cdots {Q_n}^{\mathrm{T}} R \)

直交行列の掛け算は直交行列になるので、行列 \(A\) は直交行列 \(Q\) と上三角行列 \(R\) に分解することができました。 

Q * A = [ cos(r)  sin(r)    *  [ a  b     = [ ...,                ... 
         -sin(r)  cos(r) ]       c  d ]       c*cos(r)-a*sin(r)   ... ]

c * cos(r) - a * sin(r) = 0

c / a = sin(r) / cos(r) = tan(r),  r = atan(c / a) になるが、

ここで e = √(a2 + c2) とすると

sin(r) = c / e, cos(r) = a / e

QR 分解は「ギブンス回転」、「ハウスホルダー変換 (Householder transformation)」、「グラム・シュミット分解」などの手法を使って行うことができます。ここでは実装が簡単なギブンス回転を使って説明します。ギブンス回転の場合、直交行列 \(Q_i\) を掛け算するとき、対角成分よりも下の成分の一つ \(a_{ij} \ (i \lt j)\) が 0 になるように角度を決めます。上図を見てください。

正方行列 \(A\) が二次元行列 [[a, b], [c, d]] の場合、左下の要素 c を 0 にすれば、\(A\) は上三角行列になります。この場合、条件は tan(r) = c / a になりますが、三辺の直角三角形 a, c, e (= √(a*a + c* c)) を考えると、sin(r) の値は c / e に、cos(r) の値は a / e になります。

それでは実際に試してみましょう。 

julia> function givens(a, c)
       d = sqrt(a * a + c * c)
       [a/d c/d; -c/d a/d]
       end
givens (generic function with 1 method)

julia> a = [2. 1.; 1. 3.]
2×2 Matrix{Float64}:
 2.0  1.0
 1.0  3.0

julia> q = givens(2, 1)
2×2 Matrix{Float64}:
  0.894427  0.447214
 -0.447214  0.894427

julia> q * a
2×2 Matrix{Float64}:
 2.23607  2.23607
 0.0      2.23607

確かにギブンス回転を使って正方行列を上三角行列に変換することができました。

●ギブンス回転による QR 分解

それではギブンス回転を使って正方行列を QR 分解するプログラムを作りましょう。次のリストを見てください。 

リスト : QR 分解 (ギブンス回転, 効率は悪い)

using LinearAlgebra

function qr_g(a::Matrix{Float64})
    n = size(a, 1)
    qs = Array{Float64}(I, n, n)   # 直交行列
    for y = 1 : n
        for x = y + 1 : n
            # println("$x, $y, $(a[x, y])")
            d = sqrt(a[x, y] ^ 2 + a[y, y] ^ 2)
            c = a[y, y] / d
            s = a[x, y] / d
            # 直交行列の生成
            xs = Array{Float64}(I, n, n)
            xs[x, x] = xs[y, y] = c
            xs[x, y] = -s
            xs[y, x] = s
            #
            a = xs * a
            qs = qs * xs'
        end
    end
    qs, a
end

直交行列は変数 qs に保持します。あとは二重の for ループの中で、直交行列 xs を作って対角線よりも下の成分を 0 にしていくだけです。

簡単な実行例を示します。ご参考までに、0 にする成分の座標と値を表示しています。 

julia> a = [2. 1.; 1. 3.]
2×2 Matrix{Float64}:
 2.0  1.0
 1.0  3.0

julia> x, y = qr_g(a)
2, 1, 1.0
([0.894427 -0.447214; 0.447214 0.894427], [2.23607 2.23607; 0.0 2.23607])

julia> x
2×2 Matrix{Float64}:
 0.894427  -0.447214
 0.447214   0.894427

julia> y
2×2 Matrix{Float64}:
 2.23607  2.23607
 0.0      2.23607

julia> x * y
2×2 Matrix{Float64}:
 2.0  1.0
 1.0  3.0

julia> b = [2. 1.; 1. 2.]
2×2 Matrix{Float64}:
 2.0  1.0
 1.0  2.0

julia> x, y = qr_g(b)
2, 1, 1.0
([0.894427 -0.447214; 0.447214 0.894427], [2.23607 1.78885; 0.0 1.34164])

julia> x
2×2 Matrix{Float64}:
 0.894427  -0.447214
 0.447214   0.894427

julia> y
2×2 Matrix{Float64}:
 2.23607  1.78885
 0.0      1.34164

julia> x * y
2×2 Matrix{Float64}:
 2.0  1.0
 1.0  2.0

julia> c = [1. 4. 5.; 4. 2. 6.; 5. 6. 3.]
3×3 Matrix{Float64}:
 1.0  4.0  5.0
 4.0  2.0  6.0
 5.0  6.0  3.0

julia> x, y = qr_g(c)
2, 1, 4.0
3, 1, 5.0
3, 2, 1.5718104959867523
([0.154303 0.801784 0.57735; 0.617213 -0.534522 0.57735; 0.771517 0.267261 -0.57735], 
 [6.48074 6.48074 6.78935; 1.86555e-16 3.74166 1.60357; -4.03004e-16 0.0 4.6188])

julia> x
3×3 Matrix{Float64}:
 0.154303   0.801784   0.57735
 0.617213  -0.534522   0.57735
 0.771517   0.267261  -0.57735

julia> y
3×3 Matrix{Float64}:
  6.48074      6.48074  6.78935
  1.86555e-16  3.74166  1.60357
 -4.03004e-16  0.0      4.6188

julia> x * y
3×3 Matrix{Float64}:
 1.0  4.0  5.0
 4.0  2.0  6.0
 5.0  6.0  3.0

julia> d = ones(4, 4)
4×4 Matrix{Float64}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> d += diagm(0 => [5,6,7,8])
4×4 Matrix{Float64}:
 6.0  1.0  1.0  1.0
 1.0  7.0  1.0  1.0
 1.0  1.0  8.0  1.0
 1.0  1.0  1.0  9.0

julia> x, y = qr_g(d)
2, 1, 1.0
3, 1, 1.0
4, 1, 1.0
3, 2, 0.640056896475198
4, 2, 0.6234292001593574
4, 3, 0.45309952467475356
([0.960769 -0.192327 -0.151767 -0.13; 0.160128 0.972948 -0.126473 -0.108334;
  0.160128 0.0905068 0.978542 -0.0928575; 0.160128 0.0905068 0.0585347 0.981194], 
 [6.245 2.40192 2.56205 2.72218; 1.04511e-17 6.79932 1.59518 1.68569;
  1.10272e-16 -1.08046e-16 7.60863 1.22711; -7.54246e-18 -3.42625e-17
 -3.0464e-18 8.49955])

julia> x
4×4 Matrix{Float64}:
 0.960769  -0.192327   -0.151767   -0.13
 0.160128   0.972948   -0.126473   -0.108334
 0.160128   0.0905068   0.978542   -0.0928575
 0.160128   0.0905068   0.0585347   0.981194

julia> y
4×4 Matrix{Float64}:
  6.245         2.40192       2.56205     2.72218
  1.04511e-17   6.79932       1.59518     1.68569
  1.10272e-16  -1.08046e-16   7.60863     1.22711
 -7.54246e-18  -3.42625e-17  -3.0464e-18  8.49955

julia> x * y
4×4 Matrix{Float64}:
 6.0  1.0  1.0  1.0
 1.0  7.0  1.0  1.0
 1.0  1.0  8.0  1.0
 1.0  1.0  1.0  9.0

julia> e = ones(5, 5)
5×5 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0

julia> e += diagm(0 => [6,7,8,9,10])
5×5 Matrix{Float64}:
 7.0  1.0  1.0   1.0   1.0
 1.0  8.0  1.0   1.0   1.0
 1.0  1.0  9.0   1.0   1.0
 1.0  1.0  1.0  10.0   1.0
 1.0  1.0  1.0   1.0  11.0

julia> x, y = qr_g(e)
2, 1, 1.0
3, 1, 1.0
4, 1, 1.0
5, 1, 1.0
3, 2, 0.6931032800836721
4, 2, 0.6796436827080401
5, 2, 0.666696903039527
4, 3, 0.5170969529394737
5, 3, 0.5035182907397784
5, 4, 0.39948447185327024
([0.961524 -0.175085 … -0.118963 -0.106152; 0.137361 0.973758 … -0.101968 -0.0909873;
 … ;
  0.137361 0.0839447 … 0.982748 -0.0707679; 0.137361 0.0839447 … 0.0411789 0.984432],
 [7.28011 2.47249 … 2.74721 2.88457; -1.09771e-16 7.86682 … 1.80601 1.88995; … ;
  9.33911e-18 3.8958e-17 … 9.55851 1.12556; 8.33341e-18 7.22394e-17 … 1.32137e-17
 10.4812])

julia> x
5×5 Matrix{Float64}:
 0.961524  -0.175085   -0.139272   -0.118963   -0.106152
 0.137361   0.973758   -0.119376   -0.101968   -0.0909873
 0.137361   0.0839447   0.979687   -0.0892219  -0.0796139
 0.137361   0.0839447   0.0572978   0.982748   -0.0707679
 0.137361   0.0839447   0.0572978   0.0411789   0.984432

julia> y
5×5 Matrix{Float64}:
  7.28011       2.47249       2.60985      2.74721       2.88457
 -1.09771e-16   7.86682       1.72207      1.80601       1.88995
  1.09335e-17  -1.70749e-17   8.67313      1.35131       1.40861
  9.33911e-18   3.8958e-17   -3.15735e-17  9.55851       1.12556
  8.33341e-18   7.22394e-17  -4.5869e-17   1.32137e-17  10.4812

julia> x * y
5×5 Matrix{Float64}:
 7.0  1.0  1.0   1.0   1.0
 1.0  8.0  1.0   1.0   1.0
 1.0  1.0  9.0   1.0   1.0
 1.0  1.0  1.0  10.0   1.0
 1.0  1.0  1.0   1.0  11.0

正常に動作していますね。なお、0 の精度が足りない (もっと 0 に近づけたい) 場合は、上三角行列に qr_g() を繰り返し適用すれば OK です。 

julia> d
4×4 Matrix{Float64}:
 6.0  1.0  1.0  1.0
 1.0  7.0  1.0  1.0
 1.0  1.0  8.0  1.0
 1.0  1.0  1.0  9.0

julia> q, r = qr_g(d)
([0.960769 -0.192327 -0.151767 -0.13; 0.160128 0.972948 -0.126473 -0.108334;
  0.160128 0.0905068 0.978542 -0.0928575; 0.160128 0.0905068 0.0585347 0.981194], 
 [6.245 2.40192 2.56205 2.72218; 1.04511e-17 6.79932 1.59518 1.68569;
  1.10272e-16 -1.08046e-16 7.60863 1.22711; -7.54246e-18 -3.42625e-17 -3.0464e-18
  8.49955])

julia> q1, r1 = qr_g(r)
([1.0 -1.67352e-18 -1.76576e-17 1.20776e-18; 1.67352e-18 1.0 2.21285e-17 4.61245e-18;
  1.76576e-17 -2.21285e-17 1.0 -9.7332e-19; -1.20776e-18 -4.61245e-18 9.7332e-19 1.0],
 [6.245 2.40192 2.56205 2.72218; 2.72755e-49 6.79932 1.59518 1.68569;
  -1.2326e-32 0.0 7.60863 1.22711; 1.19971e-50 0.0 0.0 8.49955])

julia> q1
4×4 Matrix{Float64}:
  1.0          -1.67352e-18  -1.76576e-17   1.20776e-18
  1.67352e-18   1.0           2.21285e-17   4.61245e-18
  1.76576e-17  -2.21285e-17   1.0          -9.7332e-19
 -1.20776e-18  -4.61245e-18   9.7332e-19    1.0

julia> r1
4×4 Matrix{Float64}:
  6.245        2.40192  2.56205  2.72218
  2.72755e-49  6.79932  1.59518  1.68569
 -1.2326e-32   0.0      7.60863  1.22711
  1.19971e-50  0.0      0.0      8.49955

julia> q * q1 * r1
4×4 Matrix{Float64}:
 6.0  1.0  1.0  1.0
 1.0  7.0  1.0  1.0
 1.0  1.0  8.0  1.0
 1.0  1.0  1.0  9.0

実用的には、ギブンス回転よりも効率が良い「ハウスホルダー変換」を用いて QR 分解を行うのが主流のようです。ハウスホルダー変換は次回以降に取り上げたいと思います。なお、Julia には行列を QR 分解する関数 LinearAlgebra.qr() が用意されているので、私たちが作る必要はありません。簡単な使用例を示します。 

julia> c
3×3 Matrix{Float64}:
 1.0  4.0  5.0
 4.0  2.0  6.0
 5.0  6.0  3.0

julia> q, r = qr(c)
LinearAlgebra.QRCompactWY{Float64,Matrix{Float64}}
Q factor:
3×3 LinearAlgebra.QRCompactWYQ{Float64,Matrix{Float64}}:
 -0.154303   0.801784  -0.57735
 -0.617213  -0.534522  -0.57735
 -0.771517   0.267261   0.57735
R factor:
3×3 Matrix{Float64}:
 -6.48074  -6.48074  -6.78935
  0.0       3.74166   1.60357
  0.0       0.0      -4.6188

julia> q * r
3×3 Matrix{Float64}:
 1.0  4.0  5.0
 4.0  2.0  6.0
 5.0  6.0  3.0

julia> d
4×4 Matrix{Float64}:
 6.0  1.0  1.0  1.0
 1.0  7.0  1.0  1.0
 1.0  1.0  8.0  1.0
 1.0  1.0  1.0  9.0

julia> q, r = qr(d)
LinearAlgebra.QRCompactWY{Float64,Matrix{Float64}}
Q factor:
4×4 LinearAlgebra.QRCompactWYQ{Float64,Matrix{Float64}}:
 -0.960769   0.192327    0.151767   -0.13
 -0.160128  -0.972948    0.126473   -0.108334
 -0.160128  -0.0905068  -0.978542   -0.0928575
 -0.160128  -0.0905068  -0.0585347   0.981194
R factor:
4×4 Matrix{Float64}:
 -6.245  -2.40192  -2.56205  -2.72218
  0.0    -6.79932  -1.59518  -1.68569
  0.0     0.0      -7.60863  -1.22711
  0.0     0.0       0.0       8.49955

julia> q * r
4×4 Matrix{Float64}:
 6.0  1.0  1.0  1.0
 1.0  7.0  1.0  1.0
 1.0  1.0  8.0  1.0
 1.0  1.0  1.0  9.0

julia> e
5×5 Matrix{Float64}:
 7.0  1.0  1.0   1.0   1.0
 1.0  8.0  1.0   1.0   1.0
 1.0  1.0  9.0   1.0   1.0
 1.0  1.0  1.0  10.0   1.0
 1.0  1.0  1.0   1.0  11.0

julia> q, r = qr(e)
LinearAlgebra.QRCompactWY{Float64,Matrix{Float64}}
Q factor:
5×5 LinearAlgebra.QRCompactWYQ{Float64,Matrix{Float64}}:
 -0.961524   0.175085    0.139272    0.118963   -0.106152
 -0.137361  -0.973758    0.119376    0.101968   -0.0909873
 -0.137361  -0.0839447  -0.979687    0.0892219  -0.0796139
 -0.137361  -0.0839447  -0.0572978  -0.982748   -0.0707679
 -0.137361  -0.0839447  -0.0572978  -0.0411789   0.984432
R factor:
5×5 Matrix{Float64}:
 -7.28011  -2.47249  -2.60985  -2.74721  -2.88457
  0.0      -7.86682  -1.72207  -1.80601  -1.88995
  0.0       0.0      -8.67313  -1.35131  -1.40861
  0.0       0.0       0.0      -9.55851  -1.12556
  0.0       0.0       0.0       0.0      10.4812

julia> q * r
5×5 Matrix{Float64}:
 7.0  1.0  1.0   1.0   1.0
 1.0  8.0  1.0   1.0   1.0
 1.0  1.0  9.0   1.0   1.0
 1.0  1.0  1.0  10.0   1.0
 1.0  1.0  1.0   1.0  11.0

●QR 分解と連立一次方程式

ところで、行列 \(A\) が正則行列であれば、QR 分解を使って連立一次方程式を解くことができます。次の式を見てください。

\(\begin{array}{l} Ax = b \\ (QR)x = b \\ Qy = b \quad (ただし y = Rx) \end{array}\)

\(Q\) は直交行列なので、y の値は \(Q^{\mathrm{T}} b\) で求めることができます。\(R\) は上三角行列なので、ガウスの消去法と同様に後退代入で x を求めることができます。ガウスの消去法はピボット選択が必要になる場合がありますが、QR 分解だとその必要がないところが利点です。

プログラムと実行例を示します。 

リスト : 連立一次方程式の解法

function qr_solver(a::Matrix{Float64}, b::Vector{Float64})
    n = size(a, 1)
    q, r = qr_g(a)
    x = q' * b
    # 後退代入
    for i = size(a, 1) : -1 : 1
        for j = i + 1 : n
            x[i] -= r[i, j] * x[j]
        end
        x[i] /= r[i, i]
    end
    x
end

julia> qr_solver([1. 1.; 2. 4.], [100., 272.])
2-element Vector{Float64}:
 64.00000000000001
 36.0

julia> qr_solver([1. 1. 1.; 2. 4. 6.; 2. 0. 4.], [10., 38., 14.])
3-element Vector{Float64}:
 3.0
 5.0
 2.0

julia> qr_solver([0. 2. 4.; 1. 1. 1.; 4. 2. 6.], [14., 10., 38.])
3-element Vector{Float64}:
 5.0
 3.0
 2.0

julia> qr_solver([1. 1. 1. 1.; -1. 1. -1. 1.; 8. 4. 2. 1.; -8. 4. -2. 1.],
 [-5., -7., -31., -35.])
4-element Vector{Float64}:
 -2.142206015326958e-16
 -9.0
  1.0000000000000002
  3.0000000000000004

julia> qr_solver([1. -1. 1. -1. 1.; 12. -6. 2. 0. 0.; 1. 1. 1. 1. 1.;
 12. 6. 2. 0. 0.; 4. 3. 2. 1. 0.], [1.,0.,8.,0.,1.])
5-element Vector{Float64}:
  0.3125000000000002
 -5.499438741189563e-17
 -1.8750000000000009
  3.5
  6.062500000000002

●QR 法の基本

次は QR 法の基本を簡単に説明します。上 (下) 三角行列の固有値は、対角行列と同様に対角成分と一致します。QR 分解は相似変換ではないので、行列 A を QR 分解したとしても、R の対角成分は A の固有値とは一致しません。そこで、次に示すような相似変換を考えます。

\( A_{i+1} = {Q_i}^{\mathrm{T}} A_i Q_i \)

ここで \(A_i\) を QR 分解して、\(A_i = Q_i R_i\) とおけば、上式は次のようになります。

\(\begin{eqnarray} A_{i+1} &=& {Q_i}^{\mathrm{T}} A_i Q_i \\ &=& {Q_i}^{\mathrm{T}} Q_i R_i Q_i \\ &=& R_i Q_i \end{eqnarray}\)

つまり、\(A_i\) を QR 分解して \(Q_i\) と \(R_i\) を求め、それを逆に掛け算したものが \(A_{i+1}\) になります。これを繰り返すと \(A_{i+1}\) の対角成分より下側の成分は 0 に収束することが知られています。つまり、\(A_{i+1}\) は上三角行列になるので、その対角成分が固有値になります。これは相似変換なので、\(A_{i+1}\) の固有値は \(A\) の固有値と一致します。

ただし、ヤコビ法と違って直交行列が固有ベクトルと一致するとは限りません。固有ベクトルを求めるときは他の方法 (たとえば、逆反復法など) を使用することになります。

これをそのままプログラムすると次のようになります。 

リスト : QR 法 (ナイーブ版)

function qr_eig(a::Matrix{Float64}; max_iter = 512, eps = 1e-8)
    d = diag(a)
    for i = 1 : max_iter
        # println("$i, $d")
        q, r = qr_g(a)
        a1 = r * q
        d1 = diag(a1)
        if sum(abs.((d1 - d) ./ d1)) < eps
            return d1
        end
        a, d = a1, d1
    end
end

対角成分 d1 と d の差分が許容誤差 eps 内に入ったら、収束したとして対角成分 (固有値) d1 を返します。あとは、QR 法の原理をそのままプログラムしただけなので、特に難しいところはないと思います。

それでは実行してみましょう。ご参考までに、収束の様子を表示しています。 

julia> a
2×2 Matrix{Float64}:
 2.0  1.0
 1.0  3.0

julia> qr_eig(a)
1, [2.0, 3.0]
2, [3.0, 2.0]
3, [3.5, 1.5]
4, [3.6, 1.4]
5, [3.61538, 1.38462]
6, [3.61765, 1.38235]
7, [3.61798, 1.38202]
8, [3.61803, 1.38197]
9, [3.61803, 1.38197]
10, [3.61803, 1.38197]
11, [3.61803, 1.38197]
12, [3.61803, 1.38197]
2-element Vector{Float64}:
 3.618033988205326
 1.3819660117946744

julia> b
2×2 Matrix{Float64}:
 2.0  1.0
 1.0  2.0

julia> qr_eig(b)
1, [2.0, 2.0]
2, [2.8, 1.2]
3, [2.97561, 1.02439]
4, [2.99726, 1.00274]
5, [2.9997, 1.0003]
6, [2.99997, 1.00003]
7, [3.0, 1.0]
8, [3.0, 1.0]
9, [3.0, 1.0]
10, [3.0, 1.0]
2-element Vector{Float64}:
 2.999999999426406
 1.0000000005735947

julia> c
3×3 Matrix{Float64}:
 1.0  4.0  5.0
 4.0  2.0  6.0
 5.0  6.0  3.0

julia> qr_eig(c)
1, [1.0, 2.0, 3.0]
2, [10.2381, -1.57143, -2.66667]
3, [12.0549, -2.5701, -3.4848]

    ・・・省略・・・

28, [12.176, -3.66868, -2.50729]
29, [12.176, -3.66868, -2.50729]
30, [12.176, -3.66868, -2.50729]
3-element Vector{Float64}:
 12.175971065046914
 -3.6686830908604544
 -2.5072879741864527

julia> d
4×4 Matrix{Float64}:
 6.0  1.0  1.0  1.0
 1.0  7.0  1.0  1.0
 1.0  1.0  8.0  1.0
 1.0  1.0  1.0  9.0

julia> qr_eig(d)
1, [6.0, 7.0, 8.0, 9.0]
2, [7.23077, 6.91233, 7.5172, 8.33971]
3, [9.00393, 6.26756, 6.94463, 7.78388]

    ・・・省略・・・

63, [10.8039, 7.50775, 6.39228, 5.29609]
64, [10.8039, 7.50775, 6.39228, 5.29609]
65, [10.8039, 7.50775, 6.39228, 5.29609]
4-element Vector{Float64}:
 10.803886359051218
  7.507748637722444
  6.392275357812142
  5.296089645414208

julia> e
5×5 Matrix{Float64}:
 7.0  1.0  1.0   1.0   1.0
 1.0  8.0  1.0   1.0   1.0
 1.0  1.0  9.0   1.0   1.0
 1.0  1.0  1.0  10.0   1.0
 1.0  1.0  1.0   1.0  11.0

julia> qr_eig(e)
1, [7.0, 8.0, 9.0, 10.0, 11.0]
2, [8.4717, 8.11519, 8.65509, 9.43996, 10.3181]
3, [10.8837, 7.45513, 8.03553, 8.87518, 9.75042]

    ・・・省略・・・

74, [13.3905, 9.54039, 8.43474, 7.35663, 6.2777]
75, [13.3905, 9.54039, 8.43474, 7.35663, 6.2777]
76, [13.3905, 9.54039, 8.43474, 7.35663, 6.2777]
5-element Vector{Float64}:
 13.390541233048985
  9.540394342247811
  8.43473663344427
  7.356631971237264
  6.277695820021704

ナイーブな QR 法の場合、固有値は対角線上に絶対値の大きい順に並びます。結果を見ればお分かりのように、ナイーブな QR 法は収束回数が多くなるという欠点があります。これを改良する方法に「原点移動法」と「減次 (デフレーション)」があります。

●原点移動法とデフレーション

QR 法の原点移動法は次式のように行います。

\(\begin{array}{l} A_{i} - \mu_{i}I = Q_{i}R_{i} \\ A_{i+1} = R_{i}Q_{i} + \mu_{i}I \end{array}\)
\(A_{i}' = A_{i} - \mu_{i}I = Q_{i}R_{i}\) とすると

\(\begin{eqnarray} A_{i+1}' &=& R_{i}Q_{i} \\ &=& Q_{i}^{\mathrm{T}} A_{i}' Q_{i} \\ &=& Q_{i}^{\mathrm{T}}(A_{i} - \mu_{i}I)Q_{i} \\ &=& Q_{i}^{\mathrm{T}} A_{i} Q_{i} - \mu_{i}Q_{i}^{\mathrm{T}}IQ_{i} \\ &=& A_{i+1} - μ_{i}I \end{eqnarray}\)

よって \(A_{i+1} = A_{i+1}' + \mu_{i}I = R_{i}Q_{i} + \mu_{i}I \)

\(A_{i}\) の対角成分から移動量 \(\mu_{i}\) を引き算して QR 分解を行い、\(A_{i+1}\) を求めるとき \(R_{i} @ Q_{i}\) の対角成分に \(\mu_{i}\) を足し算します。移動量 \(\mu\) は、右下隅 2 行 2 列の行列の固有値で、右下の要素に近い値を選びます。これを「ウィルキンソンの移動法」といいます。

プログラムは簡単に修正できるので、さっそく試してみたのですが、ウィルキンソンの移動法だけでは QR 法の収束回数を大幅に減らすことはできませんでした。右下隅の対角成分の収束は速くなるのですが、それ以外の対角成分の収束は速くなりません。

そこで「減次 (デフレーション)」という方法と組み合わせることにします。デフレーションは n * n の行列 \(A\) の固有値 \(\lambda_{n}\) をひとつ求めたら、行列の次数を一つ減らして n-1 * n-1 の行列 \(A'\) で固有値 \(\lambda_{n-1}\) を求めるという方法です。

相似変換で n * n の行列 \(A\) が上三角行列に近づいていくとき、対角成分は固有値に近づいていきます。一番下の行の成分 \(a_{n,0}, \ldots, a_{n,n-1}\) が 0 に収束すると、対角成分 \(a_{n,n}\) は固有値 \(\lambda_{n}\) に収束します。他の非対角成分が固有値 \(\lambda_{n}\) に影響を与えることはありません。

残りの固有値 \(\lambda_{1}, \ldots, \lambda_{n-1}\) は n-1 * n-1 の行列 \(A'\) を上三角行列に相似変換すれば求めることができるはずです。

そしてここがポイントですが、最初にウィルキンソンの移動法を適用し、そのあとデフレーションしたときも、その行列での右下隅 2 * 2 行列の固有値を求めてウィルキンソンの移動法を適用します。これで右下隅の対角成分 \(a_{n-1,n-1}\) の収束が速くなります。

一番下の行の成分 \(a_{n-1,0}, \ldots, a_{n-1,n-2}\) が 0 に収束したら、対角成分 \(a_{n-1,n-1}\) の値が固有値 \(\lambda_{n-1}\) になります。あとはこれを繰り返すことで、すべての固有値を高速に求めることができます。

●プログラムの作成

それではプログラムを作りましょう。最初に行列 [a b; c d] の固有値を求め、d に近い値を返す関数 eig22() を作ります。これは二次方程式を解くだけなので簡単です。

\(\begin{array}{l} \left|\begin{pmatrix} a & b \\ c & d \end{pmatrix} - \lambda I\right| = 0 \\ (a - \lambda) (d - \lambda) - bc = 0 \\ \lambda^2 - (a + d)\lambda + ad - bc \\ \lambda = \dfrac{a + d \pm \sqrt{(a + d)^2 - 4(ad - bc)}}{2} \end{array}\) 
リスト : シフト値を求める

function eig22(a, c, b, d)
    e = a + d
    f = sqrt(e * e - 4 * (a * d - b * c))
    k1 = (e + f) / 2
    k2 = (e - f) / 2
    abs(d - k1) < abs(d - k2) ? k1 : k2
end

次はウィルキンソンの移動法とデフレーションを行う関数 qr_eig_shiftd() を作ります。 

リスト : ウィルキンソンの移動法とデフレーション

function qr_eig_shiftd(a::Matrix{Float64}; max_iter = 1024, e = 1e-14)
    n = size(a, 1)
    ks = Float64[]
    i = 0
    for s in n : -1 : 2
        a = @view a[1:s, 1:s]
        while true
            # シフト値の計算
            # [a b; c d]... は (a, c, b, d) に展開される
            k = eig22(a[end-1:end, end-1:end]...)
            dk = diagm(0 => fill(k, s)) 
            q, r = qr_g(a - dk)
            a = r * q + dk
            # println("$i, $s, $k")
            # println(a)
            i += 1
            if all(x -> abs(x) < e, @view a[end, 1:end-1])
                push!(ks, a[end, end])
                break
            elseif i >= max_iter
                error("Repeated Over!!")
            end
        end
    end
    push!(ks, a[1, 1])
    ks
end

引数 a が固有値を求める行列、max_iter が繰り返しの最大値、e が 0 とみなす値です。変数 ks は固有値を格納する配列です。最初の for ループで、行列 a の次元数 (s) を一つずつ減らしていきます。Julia の場合、a の大きさは a[1:s, 1:s] で簡単に変更することができます。次の while ループで相似変換を行います。まず、シフト値を関数 eig22() で求めて変数 k にセットし、変数 dk に対角成分が k の対角行列をセットします。それから、gr_g() で行列 a - dk を QR 分解して、r * q + dk で相似変換します。

次の if 文で、最終行の対角成分以外の要素の絶対値が e 未満なったら収束と判定し、固有値 a[end, end] を ks に追加してから break で while ループを脱出します。なお、a[end, end] の左隣の要素 a[end - 1, end] だけで収束を判定する方法もあるようです。最初の for ループを抜けたとき、行列 a は 2 行 2 列になっています。a[2, 2] の固有値は ks に追加していますが、a[1, 1] の固有値は追加していません。push!() で a[1, 1] を追加してから ks を返します。

●簡単な実行例

プログラムはこれで完成です。簡単な実行例を示します。ご参考までに、収束の様子を表示しています。 

julia> a
2×2 Matrix{Float64}:
 2.0  1.0
 1.0  3.0

julia> qr_eig_shiftd(a)
0, 2, 3.618033988749895
[1.38197 0.0;
 5.83679e-17 3.61803]
2-element Vector{Float64}:
 3.618033988749895
 1.381966011250105

julia> b
2×2 Matrix{Float64}:
 2.0  1.0
 1.0  2.0

julia> qr_eig_shiftd(b)
0, 2, 1.0
[3.0 0.0;
 0.0 1.0]
2-element Vector{Float64}:
 1.0
 2.9999999999999996

julia> c
3×3 Matrix{Float64}:
 1.0  4.0  5.0
 4.0  2.0  6.0
 5.0  6.0  3.0

julia> qr_eig_shiftd(c)
0, 3, 8.520797289396148
[-2.18468 0.36261 2.69197;
 0.36261 -2.23785 4.17335;
 2.69197 4.17335 10.4225]
1, 3, 11.674435824128324
[-2.71117 -0.441417 0.101376;
 -0.441417 -3.46269 0.152193;
 0.101376 0.152193 12.1739]
2, 3, 12.17533752692046
[-2.73873 -0.463931 4.31264e-6;
 -0.463931 -3.43724 5.92369e-6;
 4.31264e-6 5.92369e-6 12.176]
3, 3, 12.175971065045754
[-2.76824 -0.48474 -6.0748e-16;
 -0.48474 -3.40773 1.87499e-17;
 3.32874e-19 4.17058e-19 12.176]
4, 2, -3.668683097953256
[-2.50729 4.996e-16;
 3.026e-16 -3.66868]
3-element Vector{Float64}:
 12.175971065046904
 -3.6686830979532563
 -2.5072879670936308

julia> d
4×4 Matrix{Float64}:
 6.0  1.0  1.0  1.0
 1.0  7.0  1.0  1.0
 1.0  1.0  8.0  1.0
 1.0  1.0  1.0  9.0

julia> qr_eig_shiftd(d)
0, 4, 9.618033988749895
[5.39664 0.322661 0.322317 -0.372566;
 0.322661 6.55309 0.552499 -0.638635; 
 0.322317 0.552499 8.08194 -1.25061;
 -0.372566 -0.638635 -1.25061 9.96834]
1, 4, 10.59155240899714
[5.32345 0.170189 0.0689349 0.0174575;
 0.170189 6.39147 0.158565 0.040156;
 0.0689349 0.158565 7.48633 0.123162;
 0.0174575 0.040156 0.123162 10.7987]
2, 4, 10.803318782718549
[5.31262 0.133156 0.0412196 -1.80854e-6;
 0.133156 6.39022 0.120797 -5.30007e-6;
 0.0412196 0.120797 7.49327 -2.16427e-5;
 -1.80854e-6 -5.30007e-6 -2.16427e-5 10.8039]
3, 4, 10.80388635904191
[5.30634 0.105319 0.0248049 1.55219e-16;
 0.105319 6.39008 0.0918731 -2.46633e-16;
 0.0248049 0.0918731 7.49969 -8.65056e-17;
 3.07492e-18 1.13884e-17 6.19408e-17 10.8039]
4, 3, 7.507241194881107
[5.29858 0.0522317 5.73353e-6;
 0.0522317 6.38978 4.27867e-5;
 5.73353e-6 4.27867e-5 7.50775]
5, 3, 7.507748705336563
[5.29673 0.0263884 -7.41171e-17;
 0.0263884 6.39164 1.4495e-15;
 7.02825e-17 1.04309e-15 7.50775]
6, 2, 6.392275290272989
[5.29609 1.52656e-16;
 1.38154e-16 6.39228]
4-element Vector{Float64}:
 10.803886359051248
  7.5077487053636505
  6.392275290272984
  5.2960896453121205

julia> e
5×5 Matrix{Float64}:
 7.0  1.0  1.0   1.0   1.0
 1.0  8.0  1.0   1.0   1.0
 1.0  1.0  9.0   1.0   1.0
 1.0  1.0  1.0  10.0   1.0
 1.0  1.0  1.0   1.0  11.0

julia> qr_eig_shiftd(e)
0, 5, 11.618033988749895
[6.40993 0.362173 0.363786 0.432553 0.391458;
 0.362173 7.54974 0.552183 0.656562 0.594186;
 0.363786 0.552183 8.89704 1.06661 0.96528;
 0.432553 0.656562 1.06661 10.7144 1.5515;
 0.391458 0.594186 0.96528 1.5515 11.4289]
1, 5, 12.663751162618677
[6.31314 0.189469 0.102528 0.0532872 0.0623671;
 0.189469 7.3679 0.199083 0.10347 0.121101;
 0.102528 0.199083 8.44433 0.230932 0.270282;
 0.0532872 0.10347 0.230932 9.6815 0.797625;
 0.0623671 0.121101 0.270282 0.797625 13.1931]
2, 5, 13.365811481974685
[6.30051 0.153776 0.0665952 0.0210645 0.000224358;
 0.153776 7.36025 0.156011 0.0493473 0.000525598;
 0.0665952 0.156011 8.42724 0.135139 0.00143937;
 0.0210645 0.0493473 0.135139 9.52147 0.00555418;
 0.000224358 0.000525598 0.00143937 0.00555418 13.3905]
3, 5, 13.390540613068227
[6.29326 0.127775 0.0458363 0.0114741 1.96132e-11;
 0.127775 7.35818 0.128488 0.0321642 5.498e-11;
 0.0458363 0.128488 8.42903 0.107399 1.83582e-10;
 0.0114741 0.0321642 0.107399 9.52899 9.04228e-10;
 1.96132e-11 5.49798e-11 1.83582e-10 9.04228e-10 13.3905]
4, 5, 13.390541233048944
[6.28851 0.106947 0.031705 0.00623607 1.1716e-17;
 0.106947 7.35706 0.105853 0.0208203 -2.77678e-16;
 0.031705 0.105853 8.43037 0.0846492 2.48374e-16;
 0.00623607 0.0208203 0.0846492 9.53352 -8.33369e-17;
 2.03421e-26 6.60595e-26 2.69392e-25 1.57631e-24 13.3905]
5, 4, 9.539981182616382
[6.28228 0.0700965 0.010614 -7.94558e-7;
 0.0700965 7.35481 0.0537255 -4.02187e-6;
 0.010614 0.0537255 8.43197 -3.23372e-5;
 -7.94558e-7 -4.02187e-6 -3.23372e-5 9.54039]
6, 4, 9.540394425673295
[6.27972 0.0467165 0.00360448 1.516e-16;
 0.0467165 7.35531 0.0274146 3.00144e-16;
 0.00360448 0.0274146 8.43403 8.33901e-16;
 3.61641e-18 2.75049e-17 4.35461e-16 9.54039]
7, 3, 8.434726012260077
[6.2782 0.0233284 1.7822e-8;
 0.0233284 7.35613 2.72068e-7;
 1.7822e-8 2.72068e-7 8.43474]
8, 3, 8.434736666495779
[6.27782 0.0116638 1.57133e-16;
 0.0116638 7.35651 1.47204e-16;
 3.66213e-23 1.2472e-21 8.43474]
9, 2, 7.356631854844217
[6.2777 -5.72459e-17;
 -4.45688e-17 7.35663]
5-element Vector{Float64}:
 13.390541233048951
  9.540394425688127
  8.434736666495784
  7.35663185484422
  6.2776958199229265

結果を見ればお分かりのように、ウィルキンソンの移動法とデフレーションを適用することで、QR 法の収束回数を大幅に減少させることができました。ここまで効果が高いとは M.Hiroi も大変驚きました。なお、実用的には行列 A をそのまま QR 分解するのではなく、「三重対角行列」に相似変換してから QR 法を適用したほうが良いようです。これは次回以降に取り上げたいと思います。

●参考文献・URL

  1. 奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
  2. 行列の固有値問題 (PDF), (桂田祐史さん)
  3. 固有値問題ノートの補足, (桂田祐史さん)
  4. 固有値問題 (1) 対称行列の固有値, (fussy さん)
  5. QR法 - Wikipedia
  6. QR分解 - Wikipedia
初版 2019 年 1 月 12 日
改訂 2021 年 12 月 5 日