Julia Language Programming

線形代数編

●ハウスホルダー変換の基本

ハウスホルダー変換は簡単に言うと「鏡映変換」のことになります。たとえば、ベクトル [x, y] を x 軸で鏡映変換すると [x, -y] になります。また、y 軸で鏡映変換すると [-x, y] になります。このとき、ベクトルの大きさは変わらないことに注意してください。

一般に、ベクトル \(X = [x_1, x_2, \ldots, x_n]\) をベクトル \(Y = [y_1, y_2, \ldots, y_n]\) に鏡映変換するとき、ベクトル \(V \ (= X - Y)\) に垂直で、\(X\) と \(Y\) の角度を二等分する超平面が鏡になります。

        Ｘ
      ／│
    ／  │
  ／    │
Ｏ───Ａ── 鏡映面
  ＼    │
    ＼  │
      ＼│
        Ｙ

\(V\) の中点を A, 原点を O とすると、O-X-A は直角三角形になるので、X から A までの長さは \(|X|\cos r\), (r は角 O-X-A) になります。内積 \((V, X)\) の値は \(|X||V| \cos r\) になるので、 ベクトル \(Y\) は次式のように表すことができます。

ここで、ベクトル \(X, Y, V\) を n 行 1 列の行列と考えると、 \((V, X)V = V(V, X) = V(V^{\mathrm{T}}X) = (VV^{\mathrm{T}})X\) になるので、\(Y\) は次のようになります。

\(H\) を「ハウスホルダー行列」といい、\(H\) で表される線形変換を「ハウスホルダー変換」といいます。\(H\) は対称行列でかつ直交行列になります。

ハウスホルター変換で特に重要なのが \(X = [x_1, x_2, \ldots, x_n]\) を \(Y = [y_1, 0, \ldots]\), (\(y_1\) は大きさ \(|X|\) で \(x_1\) の符号を反転したもの) に移す変換です。このとき、ベクトル \(V\) は次のようになります。

V は単位ベクトルにしたほうが都合が良いので、c を次のように定義します。

ハウスホルダー行列 \(H\) を作る場合、\(c = \dfrac{1}{\sqrt{|X|(|X| \pm x_1)}}\) とすれば、\(H = I - VV^{\mathrm{T}}\) になります。

●ハウスホルダー変換のプログラム

それではプログラムを作りましょう。ハウスホルダー行列を生成する関数 householer() は次のようになります。

リスト : ハウスホルダー変換

function householder(xs::Vector{Float64})
    n = length(xs)
    x = norm(xs)
    if xs[1] < 0 x = -x end
    vs = copy(xs)
    vs[1] += x
    vs /= sqrt(vs[1] * x)
    Array{Float64}(I, n, n) - reshape(vs, n, 1) * reshape(vs, 1, n)
end

ベクトル xs の大きさを norm() で求めて変数 x にセットします。xs[1] が負ならば x も負にします。あとはベクトル vs を作って、行列 I - vsvs^T を作成して返すだけです。

それでは実行してみましょう。

julia> h = householder([0., 0., 1.])
3×3 Matrix{Float64}:
  0.0  0.0  -1.0
  0.0  1.0   0.0
 -1.0  0.0   0.0

julia> h'
3×3 adjoint(::Matrix{Float64}) with eltype Float64:
  0.0  0.0  -1.0
  0.0  1.0   0.0
 -1.0  0.0   0.0

julia> h * h
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

householder() にベクトル [0, 0, 1] を渡すと、それを [-1, 0, 0] に変換するハウスホルダー行列を返します。ハウスホルダー行列は対称行列かつ直交行列なので、h.T と h は等しくなり、h.T と h^-1 は等しくなるので、h * h は単位行列になります。

ところで、行 (または列) を交換する置換行列 R は、行列 A の左側から掛け算 (RA) すると A の行を交換します。右側から掛け算 (AR) すると列を交換します。簡単な例を示しましょう。

julia> r = [0. 0. 1.; 0. 1. 0.; 1. 0. 0.]
3×3 Matrix{Float64}:
 0.0  0.0  1.0
 0.0  1.0  0.0
 1.0  0.0  0.0

julia> x = [0. 1. 2.; 3. 4. 5.; 6. 7. 8.]
3×3 Matrix{Float64}:
 0.0  1.0  2.0
 3.0  4.0  5.0
 6.0  7.0  8.0

julia> r * x
3×3 Matrix{Float64}:
 6.0  7.0  8.0
 3.0  4.0  5.0
 0.0  1.0  2.0

julia> x * r
3×3 Matrix{Float64}:
 2.0  1.0  0.0
 5.0  4.0  3.0
 8.0  7.0  6.0

r は 1 行 (列) と 3 行 (列) を交換する置換行列です。r * x は 1 行と 3 行を交換し、x * r は 1 列と 3 列を交換しています。同様に、ハウスホルダー行列 H は行列 A の左側から掛け算 (HA) すると行に対してハウスホルダー変換を行い、右側から掛け算 (AH) すると列に対してハウスホルダー変換を行います。

簡単な例を示しましょう。

julia> [0. 0. 1.]
1×3 Matrix{Float64}:
 0.0  0.0  1.0

julia> [0. 0. 1.] * h
1×3 Matrix{Float64}:
 -1.0  0.0  0.0

julia> a = [0. 0. 1.; 0. 0. 1.; 0. 0. 1.]
3×3 Matrix{Float64}:
 0.0  0.0  1.0
 0.0  0.0  1.0
 0.0  0.0  1.0

julia> a * h
3×3 Matrix{Float64}:
 -1.0  0.0  0.0
 -1.0  0.0  0.0
 -1.0  0.0  0.0

julia> b = a'
3×3 adjoint(::Matrix{Float64}) with eltype Float64:
 0.0  0.0  0.0
 0.0  0.0  0.0
 1.0  1.0  1.0

julia> h * b
3×3 Matrix{Float64}:
 -1.0  -1.0  -1.0
  0.0   0.0   0.0
  0.0   0.0   0.0

ハウスホルダー変換は鏡映なので、[-1, 0, 0] に H を適用すれば [0, 0, 1] に変換されます。つまり、H を二回適用すれば元に戻ります。H * H = I になるので、これは当然ですね。

julia> a
3×3 Matrix{Float64}:
 0.0  0.0  1.0
 0.0  0.0  1.0
 0.0  0.0  1.0

julia> a * h
3×3 Matrix{Float64}:
 -1.0  0.0  0.0
 -1.0  0.0  0.0
 -1.0  0.0  0.0

julia> (a * h) * h
3×3 Matrix{Float64}:
 0.0  0.0  1.0
 0.0  0.0  1.0
 0.0  0.0  1.0

julia> b
3×3 adjoint(::Matrix{Float64}) with eltype Float64:
 0.0  0.0  0.0
 0.0  0.0  0.0
 1.0  1.0  1.0

julia> h * b
3×3 Matrix{Float64}:
 -1.0  -1.0  -1.0
  0.0   0.0   0.0
  0.0   0.0   0.0

julia> h * (h * b)
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 1.0  1.0  1.0

●ハウスホルダー変換による QR 分解

次はハウスホルダー変換を使って行列を QR 分解するプログラムを作りましょう。基本的な考え方は簡単です。n * n の行列 \(A\) を QR 分解する場合、左端の列が \([x, 0, \ldots, 0]\) になるようにハウスホルダー行列 \(H\) を生成して、それを \(A\) に適用します。

次は、行列 \(A\) の次数を減じて n-1 * n-1 の行列 \(A'\) = \(A\) ?? (Drop 0, Drop 0) を考えて、\(A'\) の左端の列が \([x', 0, \ldots, 0]\) になるようにハウスホルダー行列 \(H'\) を生成します。このとき、単位行列 \(I\) の右下の部分を \(H'\) に置き換えたものが実際のハウスホルダー行列 \(H\) になります。

\(A \times H\) の結果は、1 行目は \(A\) と同じになり、1 列目は \(A\) の 1 列 2 行目以降が 0 なので、\(A \times H\) の 1 列 2 行目以降も 0 になります。この行列も対称行列でかつ直交行列になることに注意してください。あとはこれを繰り返していけば QR 分解を行うことができます。

簡単な例を示しましょう。

julia> a = [1. 2. 3.; 4. 5. 6.; 7. 8. 10.]
3×3 Matrix{Float64}:
 1.0  2.0   3.0
 4.0  5.0   6.0
 7.0  8.0  10.0

julia> h = householder(a[:, 1])
3×3 Matrix{Float64}:
 -0.123091  -0.492366  -0.86164
 -0.492366   0.784146  -0.377745
 -0.86164   -0.377745   0.338946

julia> a1 = h * a
3×3 Matrix{Float64}:
 -8.12404      -9.60114    -11.9399
  4.44089e-16  -0.0859656   -0.549676
 -4.44089e-16  -0.90044     -1.46193

julia> h1 = Array{Float64}(I, 3, 3)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> h1[2:end, 2:end] = householder(a1[2:end, 2])
2×2 Matrix{Float64}:
 -0.0950385  -0.995474
 -0.995474    0.0950385

julia> h1
3×3 Matrix{Float64}:
 1.0   0.0         0.0
 0.0  -0.0950385  -0.995474
 0.0  -0.995474    0.0950385

julia> h1'
3×3 adjoint(::Matrix{Float64}) with eltype Float64:
 1.0   0.0         0.0
 0.0  -0.0950385  -0.995474
 0.0  -0.995474    0.0950385

julia> h1 * h1
3×3 Matrix{Float64}:
 1.0  0.0          0.0
 0.0  1.0          1.11022e-16
 0.0  1.11022e-16  1.0

julia> a2 = h1 * a1
3×3 Matrix{Float64}:
 -8.12404      -9.60114    -11.9399
  3.99874e-16   0.904534     1.50756
 -4.84285e-16  -1.249e-16    0.408248

julia> h * h1
3×3 Matrix{Float64}:
 -0.123091   0.904534   0.408248
 -0.492366   0.301511  -0.816497
 -0.86164   -0.301511   0.408248

julia> h * h1 * a2
3×3 Matrix{Float64}:
 1.0  2.0   3.0
 4.0  5.0   6.0
 7.0  8.0  10.0

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

リスト : QR 分解 (ハウスホルダー変換)

function qr_h(a::Matrix{Float64})
    n = size(a, 1)
    xs = Array{Float64}(I, n, n)
    for i = 1 : n - 1
        h = Array{Float64}(I, n, n)
        h[i:end, i:end] = householder(a[i:end, i])
        a = h * a
        xs = xs * h
    end
    xs, a
end

for ループで行列 a の次数を一つずつ減らしながら変換行列 h を生成して、それを a の左側から掛け算していくだけです。h は対称行列でかつ直交行列なので、直交行列 xs を更新するときは h を転置する必要はありません。

簡単な実行例を示します。

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

julia> q, r = qr_h(a)
([-0.894427 -0.447214; -0.447214 0.894427], [-2.23607 -2.23607; 2.22045e-16
  2.23607])

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

julia> r
2×2 Matrix{Float64}:
 -2.23607      -2.23607
  2.22045e-16   2.23607

julia> q * r
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> q, r = qr_h(b)
([-0.894427 -0.447214; -0.447214 0.894427], [-2.23607 -1.78885; 2.22045e-16
  1.34164])

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

julia> r
2×2 Matrix{Float64}:
 -2.23607      -1.78885
  2.22045e-16   1.34164

julia> q * r
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> q, r = qr_h(c)
([-0.154303 0.801784 -0.57735; -0.617213 -0.534522 -0.57735; -0.771517
  0.267261 0.57735], 
 [-6.48074 -6.48074 -6.78935; 5.47064e-16 3.74166 1.60357; -3.08466e-16
  0.0 -4.6188])

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

julia> r
3×3 Matrix{Float64}:
 -6.48074      -6.48074  -6.78935
  5.47064e-16   3.74166   1.60357
 -3.08466e-16   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

正常に動作していますね。ところで、前回作成した固有値を求める関数 qr_eig_shiftd() は QR 分解にギブンス回転 qr_g() を使っていますが、これを qr_h() に変えると実行速度がぐっと速くなります。

julia> a = ones(50, 50) + diagm(0 => 51:100)

・・・省略・・・

ギブンス回転 (qr_g() を使用)

julia> @time qr_eig_shiftd(a)
  1.005769 seconds (168.77 k allocations: 1.563 GiB, 11.60% gc time,
 0.64% compilation time)
50-element Vector{Float64}:
  99.74906780730886
  98.69950593267251
  97.6661476787786
  96.63999546636549
  95.61808475164233
  94.59902628563236
  93.58204186546503
  92.56664498361492
  91.55250785178008
  90.53939725270898
  89.52714021324178
  88.51560415555075
  87.50468470674302
   ⋮
  61.29957929912226
  60.29156988164414
  59.28324320168637
  58.274512288814
  57.26525905248004
  56.25531558038271
  55.244428618584614
  54.23218454696027
  53.21782509008516
  52.199674276162
  51.17236607080133
 129.59687693462823

ハウスホルダー変換 (qr_h() を使用)

julia> @time qr_eig_shiftd(a)
  0.091721 seconds (33.38 k allocations: 109.927 MiB, 8.75% gc time)
50-element Vector{Float64}:
  99.74906780730886
  98.6995059326725
  97.66614767877859
  96.63999546636548
  95.61808475164234
  94.59902628563238
  93.58204186546497
  92.56664498361489
  91.55250785178006
  90.53939725270898
  89.52714021324162
  88.51560415555063
  87.504684706743
   ⋮
  61.29957929912308
  60.291569881644186
  59.28324320168674
  58.274512288813916
  57.2652590524798
  56.255315580382494
  55.24442861858475
  54.23218454696045
  53.21782509008545
  52.19967427616198
  51.17236607080121
 129.59687693462786

実行環境 : Julia 1.10.5, Ubuntu 22.04 (WSL2), Intel Core i5-6200U 2.30GHz

gr_h() のほうが約 11 倍速くなりました。もっとも、これは qr_g() の実装がナイーブだからで、工夫次第では qr_g() でも速くすることが可能です。これは「三重対角行列の QR 分解」で取り上げます。

●三重対角化

「三重対角行列 (tridiagonal matrix)」は、対角成分とそれに隣接する成分以外が 0 の行列のことです。簡単な例を下図に示します。

[ a  b  0  0  0 ;
  c  d  e  0  0 ;
  0  f  g  h  0 ;
  0  0  i  j  k ;
  0  0  0  l  m  ]

図 : 三重対角行列

実対称行列の固有値を求める場合、三重対角行列に相似変換してから QR 法を適用することが行われます。一般的には、このほうが速くなるといわれています。実対称行列の場合、三重対角化はハウスホルダー変換を使うと簡単です。プログラムは次のようになります。

リスト : 対称行列の三重対角化

function tridiag(a::Matrix{Float64})
    n = size(a, 1)
    for i = 1 : n - 2
        h = Array{Float64}(I, n, n)
        h[i+1:end, i+1:end] = householder(a[i+1:end, i])
        a = h * a * h
    end
    a
end

QR 分解と同じように、行列 a の次数を一つずつ減らしながらハウスホルダー変換を行います。このとき、左上隅の成分はハウスホルダー変換を適用しないことに注意してください。左端の列であれば 2 行目以降の成分に、一番上の行であれば 2 列目以降の成分にハウスホルダー変換を適用します。たとえば、行 (または列) の成分が [x1, x2, ..., xn] とすると、[x2, ..., xn] を [x2', 0, ..., 0] に変換するわけです。

対称行列の性質から、1 行目と 1 列目の成分は同じになるので、ハウスホルダー行列 h は同じものになります。つまり、a の左右に h を掛け算 (h * a * h) すれば、行と列を変換することができます。また、h は対称行列でかつ直交行列なので、h * a * h は相似変換になります。したがって、三重対角化した行列の固有値は元の行列の固有値と一致します。

簡単な実行例を示します。

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> tridiag(c)
3×3 Matrix{Float64}:
  1.0          -6.40312    8.88178e-16
 -6.40312       8.46341    0.829268
  8.88178e-16   0.829268  -3.46341

julia> d = [6. 1. 1. 1.; 1. 7. 1. 1.; 1. 1. 8. 1.; 1. 1. 1. 9.]
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> tridiag(d)
4×4 Matrix{Float64}:
  6.0          -1.73205       1.43158e-16   1.05314e-16
 -1.73205      10.0           0.816497      3.1554e-16
  1.43158e-16   0.816497      7.0          -0.57735
  1.05314e-16  -8.87036e-17  -0.57735       7.0

julia> e = [7. 1. 1. 1. 1.; 1. 8. 1. 1. 1.; 1. 1. 9. 1. 1.; 1. 1. 1. 10. 1.;
 1. 1. 1. 1. 11.]
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> tridiag(e)
5×5 Matrix{Float64}:
  7.0          -2.0          -4.34443e-16  2.77556e-16  -1.24127e-16
 -2.0          12.5          -1.11803      4.37198e-16  -4.19456e-17
 -4.34443e-16  -1.11803       8.5          0.894427     -5.11549e-17
  2.77556e-16  -2.18951e-17   0.894427     8.5           0.67082
 -1.24127e-16   8.03365e-17  -1.75389e-17  0.67082       8.5

●三重対角化の高速化

さきほどのプログラムは、変換行列 h の大きさを行列 a と同じ大きさのままで計算しました。変換行列 h にも減次 (デフレーション) を適用すると、処理を高速化することができます。

行列 a にハウスホルダー変換による相似変換を適用すると、\(a_{0,0}, a_{0,1}, a_{1,0}\) 以外の成分 \(a_{i,0}, a_{0,i}\) は 0 になります。次に相似変換を適用すると \(a_{0,0}, a_{0,1}, a_{1,0}\) の値は変化せず、それ以外の成分 \(a_{i,0}\) と \(a_{0,i}\) は 0 のままです。つまり、1 行目と 1 列目の値は変化しないので、次数を一つ減らした行列 a' = a[1:, 1:]に相似変換を適用すればいいことになります。これをプログラムすると次のようになります。

リスト : 三重対角化の高速化

function tridiag1(a::Matrix{Float64})
    n = size(a, 1)
    for i = 1 : n - 2
        h = Array{Float64}(I, n - i + 1, n - i + 1)
        h[2:end, 2:end] = householder(a[i+1:end, i])
        a[i:end, i:end] = h * a[i:end, i:end] * h
    end
    a
end

a[i:end, i:end] の相似変換の結果を a[i:end, i:end] に代入していくだけです。この場合、引数の行列 a の値を破壊することに注意してください。行列 h の大きさが一つずつ減っていくので、大きさが n のまま演算するよりも効率的です。

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

julia> a = ones(200, 200) + diagm(0 => 101:300)

    ・・・省略・・・

julia> @time tridiag(a)
  0.399839 seconds (3.63 k allocations: 242.478 MiB, 4.59% gc time)

    ・・・省略・・・

julia> @time tridiag1(a)
  0.138611 seconds (3.85 k allocations: 143.240 MiB, 3.96% gc time)

    ・・・省略・・・

実行環境 : Julia 1.10.5, Ubuntu 22.04 (WSL2), Intel Core i5-6200U 2.30GHz

200 * 200 の行列で約 2.8 倍速くなりました。なお、参考 URL 6 『三重対角化 - [物理のかぎしっぽ]』によると、もっと効率の良い方法もあるようです。興味のある方はいろいろ試してみてください。

●三重対角行列の QR 分解

三重対角行列に変換すると、QR 分解も簡単になります。対角成分 \(a_{i,i}\) の 1 行下の成分 \(a_{i+1,i}\) を 0 にするだけです。ギブンス回転を使うと、プログラムは次のようになります。

リスト : 三重対角行列の QR 分解

function qr_tri(a::Matrix{Float64})
    n = size(a, 1)
    qs = Array{Float64}(I, n, n)   # 直交行列
    for y = 1 : n - 1
        x = y + 1
        d = sqrt(a[x, y] ^ 2 + a[y, y] ^ 2)
        c = a[y, y] / d
        s = a[x, y] / d
        a1 = a[y, :]
        a2 = a[x, :]
        a[y, :] = a1 * c    + a2 * s
        a[x, :] = a1 * (-s) + a2 * c
        q1 = qs[:, y]
        q2 = qs[:, x]
        qs[:, y] = q1 * c    + q2 * s
        qs[:, x] = q1 * (-s) + q2 * c
    end
    qs, a
end

ギブンス回転で行列を QR 分解する関数 qr_g() は二重の for ループになりますが、三重対角行列は一重の for ループで済みます。それから、a[y + 1, y] を 0 にする場合、ギブンス回転で値が変わるのは y, y + 1 行だけ、直交行列 qs では y, y + 1 列だけです。ここだけ計算すれば、処理をさらに高速化することができます。たとえば、4 行 4 列の行列 a で a[0, 1] を 0 にする場合は下図のようになります。

s = sin(r), c = cos(r) とする

[[ c, s, 0, 0],    [[a11, a12, a13, a14],     [[ ca11 + sa21,  ca12 + sa22,  ca13 + sa23,  ca14 + sa24],
 [-s, c, 0, 0], *   [a21, a22, a23, a24],  =   [-sa11 + ca21, -sa12 + ca22, -sa13 + ca23, -sa14 + ca24],
 [ 0, 0, 1, 0],     [a31, a32, a33, a34],      [ ...                                                  ],
 [ 0, 0, 0, 1]]     [a41, a42, a43, a44]]      [ ...                                                  ]]

[[q11, q12, q13, q14],    [[c, -s, 0, 0],     [[q11c+q12s, -q11s+q12c, ... ],
 [q21, q22, q23, q24], *   [s,  c, 0, 0],  =   [q21c+q22s. -q21s+q22c, ... ],
 [q31, q32, q33, q34],     [0,  0, 1, 0],      [q31c+q32s. -q31s+q32c, ... ],
 [q41, q42, q43, q44]]     [0,  0, 0, 1]]      [q41c+q42s. -q41s+q42c, ... ]],

一般に、a[i+1, i] を 0 にする場合、a の要素と直交行列 q の要素は次のようになります。

ai,j = cos(r) * ai,j + sin(r) * ai+1,j
ai+1,j = -sin(r) * ai,j + cos(r) * ai+1,j
qj,i = cos(r) * qj,i + sin(r) * qj,i+1
qj,i+1 = -sin(r) * qj,i + cos(r) * qj,i+1
j = 0, 1, ..., n-1 (n は行列 a の次元数)

関数 qr_tri() は、これをそのままプログラムしています。ただし、a[y, :], q[:, y] を書き換えると、a[x, :], q[:, x] (x = y + 1) を計算できなくなるので、コピーしていることに注意してください。なお、三重対角行列の特徴を使えば、さらに計算量を減らすことができると思います。興味のある方はプログラムを改良してみてください。

簡単な実行例を示します。

julia> x = diagm(0 => [1.,1.,1.,1.,1.], 1 => [1.,1.,1.,1.], -1 => [1.,1.,1.,1.])
5×5 Matrix{Float64}:
 1.0  1.0  0.0  0.0  0.0
 1.0  1.0  1.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  0.0  1.0  1.0

julia> q, r = qr_tri(x)
([0.707107 0.0 … 0.288675 0.5; 0.707107 0.0 … -0.288675 -0.5; … ;
  0.0 0.0 … 0.288675 0.5; 0.0 0.0 … 0.866025 -0.5], 
 [1.41421 1.41421 … 0.0 0.0; 0.0 1.0 … 1.0 0.0; … ; 
  0.0 0.0 … 1.1547 1.1547; 0.0 0.0 … 0.0 0.0])

julia> q
5×5 Matrix{Float64}:
 0.707107  0.0  -0.408248   0.288675   0.5
 0.707107  0.0   0.408248  -0.288675  -0.5
 0.0       1.0   0.0        0.0        0.0
 0.0       0.0   0.816497   0.288675   0.5
 0.0       0.0   0.0        0.866025  -0.5

julia> r
5×5 Matrix{Float64}:
 1.41421  1.41421   0.707107     0.0       0.0
 0.0      1.0       1.0          1.0       0.0
 0.0      0.0       1.22474      0.816497  0.816497
 0.0      0.0      -5.55112e-17  1.1547    1.1547
 0.0      0.0      -9.61481e-17  0.0       0.0

julia> q * r
5×5 Matrix{Float64}:
 1.0  1.0  -1.54575e-16   4.39774e-17   4.39774e-17
 1.0  1.0   1.0          -4.39774e-17  -4.39774e-17
 0.0  1.0   1.0           1.0           0.0
 0.0  0.0   1.0           1.0           1.0
 0.0  0.0   0.0           1.0           1.0

●三重対角行列による QR 法

最後に、三重対角化による QR 法のプログラムを作ります。まず前提として、三重対角行列を QR 分解して相似変換 (R * Q) する場合、三重対角行列の性質は残ることに注意してください。つまり、主対角線と副対角線以外は 0 のままで、相似変換を繰り返すと副対角線が 0 に近づいていくことになります。したがって、QR 分解は qrTri で大丈夫です。収束判定も対角成分 \(a_{i,i}\) の左隣 \(a_{i.i-1}\) が 0 になったかチェックすれば OK です。

プログラムは次のようになります。

リスト : 三重対角化による QR 法

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

前回作成した関数 qr_eig_shiftd() とほぼ同じですが、最初に引数の行列 a を tridiag1() で三重対角化するところと、QR 分解に関数 qr_tri() を使うところが異なります。QR 分解は while ループの中で行っているので、三重対角行列専用の関数 qr_tri() を使うことで高速化が期待できます。これはあとで試してみましょう。

簡単な実行例を示します。ご参考までに、収束の様子を表示しています。

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_tri(c)
0, 3, -3.5207972893961483
[12.0145 -1.53144 3.33067e-16; -1.53144 -2.35012 -0.0708263; 4.59423e-16
 -0.0708263 -3.66439]
1, 3, -3.668197833676132
[12.1751 -0.11324 4.33681e-16; -0.11324 -2.50641 2.95525e-5; 4.55054e-16
 2.95525e-5 -3.66868]
2, 3, -3.6686830979527425
[12.176 -0.00830084 4.35744e-16; -0.00830084 -2.50728 -3.04512e-16;
 4.55042e-16 -1.00509e-17 -3.66868]
3, 2, -2.5072879670936414
[12.176 1.15533e-15; 8.40335e-20 -2.50729]
3-element Vector{Float64}:
 -3.668683097953265
 -2.5072879670936414
 12.175971065046898

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_tri(d)
0, 4, 6.422649730810372
[8.97741 -2.58338 1.2254e-15 -1.10427e-15; 
 -2.58338 7.13358 0.232519 1.95883e-16;
 -5.88644e-17 0.232519 7.49568 -0.0345458;
 -9.6328e-17 -1.5938e-16 -0.0345458 6.39333]
1, 4, 6.392247267757526
[10.6395 -0.934698 1.14503e-15 8.84871e-16;
 -0.934698 5.47713 0.196186 -6.75546e-16;
 -6.9953e-17 0.196186 7.49108 8.68092e-7;
 -4.23763e-17 -1.80015e-16 8.68092e-7 6.39228]
2, 4, 6.39227529027296
[10.7934 -0.239022 1.26967e-15 -1.00938e-15;
 -0.239022 5.32413 0.196764 -4.13258e-16;
 -5.40844e-17 0.196764 7.49016 3.55952e-16;
 2.6955e-18 -1.82016e-16 3.26023e-17 6.39228]
3, 3, 7.507887645813202
[10.7992 -0.159927 -1.25421e-15;
 -0.159927 5.30074 1.24529e-5;
 5.37264e-17 1.24529e-5 7.50775]
4, 3, 7.507748705363892
[10.8018 -0.107359 1.28037e-15;
 -0.107359 5.29818 5.09497e-16;
 -5.36629e-17 3.97996e-18 7.50775]
5, 2, 5.296089645312119
[10.8039 1.94289e-16;
 -3.58317e-17 5.29609]
4-element Vector{Float64}:
  6.392275290272984
  7.507748705363649
  5.296089645312121
 10.803886359051255

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_tri(e)
0, 5, 7.829179606750032
[10.2782 -3.51428 9.16652e-16 4.72201e-16 -2.74867e-16;
 -3.51428 9.40004 -0.30337 4.33617e-17 4.46966e-16;
 1.58278e-16 -0.30337 9.19523 0.67934 -5.93262e-16;
 6.42852e-17 3.32996e-16 0.67934 7.9693 0.542268;
 -8.42689e-18 9.23703e-17 6.67436e-18 0.542268 8.15726]
1, 5, 8.613630062785017
[12.28 -2.58071 4.84232e-16 1.17124e-17 -6.0711e-16;
 -2.58071 7.39104 -0.0957376 -1.7743e-16 2.31053e-16;
 -1.94822e-16 -0.0957376 9.17079 0.813238 -1.97839e-16;
 1.22741e-16 -1.56577e-17 0.813238 7.73074 -0.108287;
 -1.72109e-16 -2.84358e-17 -1.66219e-17 -0.108287 8.4274]
2, 5, 8.443844687781366
[13.1472 -1.29275 5.67687e-16 5.38907e-16 -5.81562e-16;
 -1.29275 6.52181 -0.0451635 6.6733e-17 1.62917e-16;
 -3.57179e-17 -0.0451635 9.18163 0.808674 -3.019e-16;
 -1.91727e-16 1.1286e-16 0.808674 7.71461 0.000912412;
 -1.46771e-16 1.32153e-16 -7.40306e-18 0.000912412 8.43474]
3, 5, 8.434737300472122
[13.3432 -0.578635 6.87025e-16 2.72997e-19 -6.03941e-16;
 -0.578635 6.32529 -0.0227408 -1.46325e-16 3.09889e-18;
 -1.82312e-16 -0.0227408 9.19696 0.794871 -1.33005e-16;
 1.19002e-16 4.36856e-17 0.794871 7.69986 -5.34351e-10;
 -1.76603e-16 -8.84816e-17 -6.20366e-18 -5.34351e-10 8.43474]
4, 5, 8.434736666495763
[13.3815 -0.253211 5.19705e-16 5.1723e-16 -6.0015e-16;
 -0.253211 6.28678 -0.0115807 1.23319e-16 6.90239e-17;
 -4.37249e-17 -0.0115807 9.21143 0.78109 -2.94269e-16;
 -2.08929e-16 5.50762e-17 0.78109 7.68554 9.89506e-17;
 -1.6503e-16 1.08608e-16 -3.504e-18 -3.65317e-18 8.43474]
5, 4, 7.356613159845852
[13.3903 -0.0453323 6.58627e-16 -2.77847e-16;
 -0.0453323 6.27813 -0.021589 -5.34568e-18;
 -1.2272e-16 -0.021589 9.54025 -7.87266e-6;
 1.77329e-16 4.33488e-17 -7.87266e-6 7.35663]
6, 4, 7.356631854844208
[13.3905 -0.00810684 6.69703e-16 2.77796e-16;
 -0.00810684 6.27829 -0.0436893 -1.99348e-17;
 -1.22831e-16 -0.0436893 9.53981 6.24595e-16;
 -1.76998e-16 4.46709e-17 9.11142e-19 7.35663]
7, 3, 9.540394428749392
[13.3905 -0.0068693 -6.63505e-16;
 -0.0068693 6.2777 -4.09991e-11;
 1.22797e-16 -4.09995e-11 9.54039]
8, 3, 9.540394425688119
[13.3905 -0.0058212 6.64168e-16;
 -0.0058212 6.2777 3.71109e-16;
 -1.22796e-16 2.1909e-19 9.54039]
9, 2, 6.2776958199229265
[13.3905 4.10262e-16;
 -1.02735e-18 6.2777]
5-element Vector{Float64}:
  8.434736666495777
  7.3566318548442124
  9.540394425688119
  6.2776958199229265
 13.390541233048957

正常に動作していますね。最後に、qr_eig_shiftd() (ハウスホルダー変換版) と qr_eig_tri() の実行時間を比較してみましょう。150 * 150 の行列の固有値を求めます。結果は次のようになりました。

julia> a = ones(150, 150) + diagm(0 => 101:250)

    ・・・省略・・・

julia> @time qr_eig_shiftd(a)
 10.942398 seconds (307.93 k allocations: 6.750 GiB, 7.11% gc time,
 0.06% compilation time)
150-element Vector{Float64}:
 249.79893116894553
 248.76416278114044
 247.74116744704168
 246.72324663493745
 245.70825981169224
 244.69522441879658
 243.68359905363496
 242.67304998141313
 241.66335515023613
 240.65435838318652
 239.64594509765175
 238.63802838744004
 237.63054054489413
 236.62342763915865
 235.61664590780654
   ⋮
 113.2385096185777
 112.23436854353221
 111.23001894672427
 110.2254242446019
 109.22053703961036
 108.21529420483611
 107.20960874273386
 106.20335544015916
 105.19634349465547
 104.18825837426269
 103.17851797702366
 102.16582074574536
 101.14591994723759
 337.79585315050826

julia> @time qr_eig_tri(a)
  0.224789 seconds (383.95 k allocations: 482.152 MiB, 7.37% gc time)
150-element Vector{Float64}:
 171.39192237836053
 172.39434240742892
 173.39677250811434
 170.38951177152444
 174.3992133358845
 175.40166555340755
 169.38710994250334
 176.40412983169202
 177.4066068512409
 178.40909730323452
 168.38471625074655
 179.41160189074273
 167.38233005846217
 180.41412132998246
 166.37995072950142
   ⋮
 107.20960874271861
 244.6952244187965
 106.20335544014831
 245.708259811692
 105.19634349464835
 246.7232466349375
 104.18825837425366
 247.74116744704153
 103.17851797701763
 248.7641627811405
 102.16582074574266
 249.7989311689453
 101.14591994723386
 337.79585315047484

実行環境 : Julia 1.10.5, Ubuntu 22.04 (WSL2), Intel Core i5-6200U 2.30GHz

qr_eig_tri() のほうが約 50 倍高速になりました。行列が大きくなると、その差はもっと大きくなるでしょう。興味のある方はいろいろ試してみてください。

●参考文献・URL

1. 奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
2. 行列の固有値問題 (PDF), (桂田祐史さん)
3. 固有値問題ノートの補足, (桂田祐史さん)
4. 固有値問題 (1) 対称行列の固有値, (fussy さん)
5. HouesHolder変換 - [物理のかぎしっぽ], (東條遼平さん)
6. 三重対角化 - [物理のかぎしっぽ], (東條遼平さん)
7. QR法 - [物理のかぎしっぽ], (東條遼平さん)
8. ハウスホルダー変換 - Wikipedia

●プログラムリスト

#
# qr.jl : QR 分解と QR 法
#
#         Copyright (C) 2019-2021 Makoto Hiroi
#
using LinearAlgebra

# ギブンス回転による QR 分解 (非効率的)
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

# 連立一次方程式の解法
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] -= dot(r[i, i+1:end], x[i+1:end])
        x[i] /= r[i, i]
    end
    x
end

# 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

# 2 行 2 列の固有値で左下の要素に近い値を選ぶ
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 法 (ウィルキンソンの移動法とデフレーション)
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
            # シフト値の計算
            k = eig22(a[end-1:end, end-1:end]...) # [a b; c d] => (a, c, b, d)
            dk = diagm(0 => fill(k, s)) 
            q, r = qr_h(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

# ハウスホルダー変換
function householder(xs::Vector{Float64})
    n = length(xs)
    x = norm(xs)
    if xs[1] < 0 x = -x end
    vs = copy(xs)
    vs[1] += x
    vs /= sqrt(vs[1] * x)
    Array{Float64}(I, n, n) - reshape(vs, n, 1) * reshape(vs, 1, n)
end

# QR 分解 (ハウスホルダー変換)
function qr_h(a::Matrix{Float64})
    n = size(a, 1)
    xs = Array{Float64}(I, n, n)
    for i = 1 : n - 1
        h = Array{Float64}(I, n, n)
        h[i:end, i:end] = householder(a[i:end, i])
        a = h * a
        xs = xs * h
    end
    xs, a
end

# 対称行列の三重対角化
function tridiag(a::Matrix{Float64})
    n = size(a, 1)
    for i = 1 : n - 2
        h = Array{Float64}(I, n, n)
        h[i+1:end, i+1:end] = householder(a[i+1:end, i])
        a = h * a * h
    end
    a
end

# 高速版
function tridiag1(a::Matrix{Float64})
    n = size(a, 1)
    for i = 1 : n - 2
        h = Array{Float64}(I, n - i + 1, n - i + 1)
        h[2:end, 2:end] = householder(a[i+1:end, i])
        a[i:end, i:end] = h * a[i:end, i:end] * h
    end
    a
end

# 三重対角行列の QR 分解
function qr_tri(a::Matrix{Float64})
    n = size(a, 1)
    qs = Array{Float64}(I, n, n)   # 直交行列
    for y = 1 : n - 1
        x = y + 1
        d = sqrt(a[x, y] ^ 2 + a[y, y] ^ 2)
        c = a[y, y] / d
        s = a[x, y] / d
        a1 = a[y, :]
        a2 = a[x, :]
        a[y, :] = a1 * c    + a2 * s
        a[x, :] = a1 * (-s) + a2 * c
        q1 = qs[:, y]
        q2 = qs[:, x]
        qs[:, y] = q1 * c    + q2 * s
        qs[:, x] = q1 * (-s) + q2 * c
    end
    qs, a
end

# 三重対角化による QR 法
function qr_eig_tri(a::Matrix{Float64}; max_iter = 1024, e = 1e-14)
    n = size(a, 1)
    a = tridiag1(copy(a))
    ks = Float64[]
    i = 0
    for s = n : -1 : 2
        a = @view a[1:s, 1:s]
        while true
            # シフト値の計算
            k = eig22(a[end-1:end, end-1:end]...) # [a b; c d] => (a, c, b, d)
            dk = diagm(0 => fill(k, s)) 
            q, r = qr_tri(a - dk)
            a = r * q + dk
            # println("$i, $s, $k")
            # println(a)
            i += 1
            if abs(a[end, end-1]) < e
                push!(ks, a[end, end])
                break
            elseif i >= max_iter
                error("Repeated Over!!")
            end
        end
    end
    push!(ks, a[1, 1])
    ks
end