部分固有値と付随する固有ベクトル（部分固有対）の精度保証付き数値計算¶

MATLABの区間演算パッケージであるINTLABには一般化固有値問題の固有値・固有ベクトルを精度保証付きで計算するverifyeigという関数が実装されている。本稿では、文献1および2を参考にverifyeig関数のアルゴリズムを説明し、一般化固有値問題

$$ Ax=\lambda Bx $$

を満たす固有値 $\lambda\in\mathbb{C}$ と固有ベクトル $x\in\mathbb{C}^n$ ($n$ は行列の次元) を求める実装を紹介する。

行列 $A,B\in \mathbb{C}^{n\times n}$ を与え、$B$ は正則、$B^{-1}A$ が対角化可能、かつ全ての固有値 $\lambda$ が重複固有値でないと仮定する。写像 $F:\mathbb{C}^{n+1}\to\mathbb{C}^{n+1}$ を

$$ F(\lambda, x) := \begin{bmatrix}x^Hx -1\\ Ax-\lambda Bx\end{bmatrix} $$

と定義すると $F(\lambda, x)=0$ をみたす $(\lambda, x)$ は一般化固有値問題の固有対となる（固有ベクトルは $\|x\|_2=1$ と正規化されている）。いま $\bar x\in \mathbb{C}^n$, $\bar \lambda \in \mathbb{C}$ を $F(\bar\lambda,\bar x)\approx 0$ とする近似解とし、 $R\in \mathbb{C}^{n+1\times n+1}$ を写像 $F$ の近似解におけるヤコビ行列の近似逆行列（$R\approx DF(\bar\lambda,\bar x)^{-1}$）とする。部分固有値・固有ベクトルの精度保証付き数値計算には次の定理を利用する。

定理 $\boldsymbol{y}\in \mathbb{IC}^n, \boldsymbol{d} \in \mathbb{IC}$ をある区間ベクトルと区間として、$\boldsymbol{w}:=(\boldsymbol{d},\boldsymbol{y}^T)^T\in \mathbb{IC}^{n+1}$ とする（${}^T$はベクトルの単なる転置を意味する）。 さらに

$$ g(\boldsymbol{w}):=z+\left(I-R\cdot DF\left(\bar\lambda + \boldsymbol{d}, \bar x + \boldsymbol{y}\right)\right)\boldsymbol{w} $$

と定義する。上の式のヤコビ行列 $DF(\lambda,x)$ と $z$ は

$$ \begin{align*} DF(\lambda,x):= \begin{bmatrix} 0 & 2x^T\\ -B\boldsymbol{x} & A-\boldsymbol{\lambda}B \\ \end{bmatrix},\quad z:=-R \left[\begin{array} a\bar x_1^2+\bar x_2^2+\dots +\bar x_n^2-1 \\ A\bar x-\bar\lambda B\bar x \\ \end{array}\right] \end{align*} $$

で与えられる。 このとき $\mathrm{int}(\boldsymbol{w})$ を区間ベクトル $\boldsymbol{w}$ の（要素ごとの）内部とすると、$g(\boldsymbol{w})\subset \mathrm{int}(\boldsymbol{w})$ ならば、一般化固有値問題 $Ax=\lambda Bx$ の真の固有値 $\lambda$ が $\bar \lambda +\boldsymbol{d}$ 内に唯一存在し、対応する固有ベクトルが $\bar x + \boldsymbol{y}$ に包含される。

この定理をもとに精度保証付き数値計算を次のように実装する。

  1. 近似固有値 $\bar{\lambda}$、近似固有ベクトル $\bar{x}$ を入力し、これらが実数であるときは実数の区間・区間ベクトルで、いずれかが複素数であるときは複素数の区間・区間ベクトルで $\boldsymbol{y}$, $\boldsymbol{d}$ を定義し、候補ベクトル $\boldsymbol{w}=(\boldsymbol{d},\boldsymbol{y}^T)^T$ を作成する。
  2. ヤコビ行列$DF$ と $z$、$g(\boldsymbol{w}):=z+\left(I-R\cdot DF\left(\bar\lambda + \boldsymbol{d}, \bar x + \boldsymbol{y}\right)\right)\boldsymbol{w}$ を区間演算により計算する。
  3. $g(\boldsymbol{w})\subsetneq \boldsymbol{w}$を満たすとき、$(\bar\lambda,\bar x)^T+g(\boldsymbol{w})$ の値を返す。

上記の手順に従って次元 $n=30$ のある行列 $A,B\in \mathbb{C}^{n\times n}$ に対し、$Ax=\lambda Bx$ を考える。 まず近似固有値$\bar{\lambda}$、近似固有ベクトル$\bar{x}$を計算する。

In [1]:
versioninfo()

Julia Version 1.10.0
Commit 3120989f39b (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 14 × Apple M3 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
  Threads: 2 on 10 virtual cores
In [2]:
using Pkg
Pkg.status("IntervalArithmetic")

Status `~/.julia/environments/v1.10/Project.toml`
  [d1acc4aa] IntervalArithmetic v0.22.6
In [3]:
using LinearAlgebra
n = 30
A, B = randn(n, n), randn(n, n)
λ, x = eigen(A, B)
Out[3]:
GeneralizedEigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
30-element Vector{ComplexF64}:
   -8.807235191171742 + 0.0im
  -1.8720739353881435 - 0.4211957322307501im
  -1.8720739353881435 + 0.42119573223075013im
  -1.2450164394864036 - 2.338661769764113im
  -1.2450164394864036 + 2.338661769764113im
  -1.1558402670064905 + 0.0im
   -0.789551856115187 - 1.002224359974505im
   -0.789551856115187 + 1.002224359974505im
  -0.5393502183046783 - 0.7116749651587355im
  -0.5393502183046783 + 0.7116749651587354im
  -0.3675440335623449 + 0.0im
 -0.22852573744671337 - 0.1592805432423365im
 -0.22852573744671337 + 0.15928054324233654im
                      ⋮
    0.281692718115899 - 0.545904155982315im
    0.281692718115899 + 0.5459041559823149im
   0.2965924537980811 - 0.8603906140438535im
   0.2965924537980811 + 0.8603906140438535im
   0.5212264660922051 + 0.0im
   0.8420124984325313 - 0.4056142428796466im
   0.8420124984325315 + 0.4056142428796465im
   0.8668146997815379 + 0.0im
   1.2621456553888593 + 0.0im
    1.560528109473799 - 1.3546291706469042im
    1.560528109473799 + 1.354629170646904im
   186.07116084856375 + 0.0im
vectors:
30×30 Matrix{ComplexF64}:
    0.20217+0.0im   0.0318037+0.0280408im   …     0.234905+0.0im
  -0.217436+0.0im    0.364885+0.312187im          -0.20174+0.0im
       -1.0+0.0im    0.211263+0.346095im              -1.0+0.0im
   -0.26971+0.0im   0.0893835+0.155071im         -0.300353+0.0im
 -0.0561925+0.0im   0.0777207+0.195086im         0.0945225+0.0im
   0.253824+0.0im   -0.176184-0.131949im    …     0.250935+0.0im
   -0.16729+0.0im    0.252789+0.191767im         0.0614289+0.0im
  -0.321578+0.0im   0.0476998+0.0385842im         -0.31054+0.0im
  -0.348441+0.0im    0.126015+0.237507im         -0.238223+0.0im
   -0.18753+0.0im    0.491601+0.191111im         0.0205939+0.0im
 -0.0727752+0.0im  -0.0939492-0.101164im    …  -0.00795106+0.0im
  -0.438972+0.0im   0.0789246+0.187319im         -0.349904+0.0im
  0.0625889+0.0im   0.0358692+0.00392634im       0.0791018+0.0im
           ⋮                                ⋱  
 -0.0121096+0.0im    0.059365+0.0364191im       -0.0724114+0.0im
  -0.146576+0.0im   -0.130795-0.116416im          -0.12067+0.0im
  -0.311531+0.0im   0.0906312+0.390606im    …    -0.176047+0.0im
 0.00415899+0.0im    0.333781+0.246889im         0.0806722+0.0im
   0.224384+0.0im   -0.270805-0.0337405im         0.210771+0.0im
  -0.207086+0.0im    0.120852+0.105772im         -0.331818+0.0im
   0.310787+0.0im    0.236172+0.0814331im         0.527194+0.0im
  -0.296509+0.0im    0.311386+0.265348im    …     -0.14192+0.0im
  -0.540491+0.0im    0.507163+0.492837im         -0.233728+0.0im
  -0.341711+0.0im  -0.0338126-0.0721759im        -0.508083+0.0im
   0.243982+0.0im  -0.0719414-0.250216im         0.0968205+0.0im
  -0.039231+0.0im    0.161358+0.112677im         0.0785697+0.0im
  1. 近似固有値 $\bar{\lambda}$、近似固有ベクトル $\bar{x}$ を入力し、これらが実数であるときは実数の区間・区間ベクトルで、いずれかが複素数であるときは複素数の区間・区間ベクトルで $\boldsymbol{y}$, $\boldsymbol{d}$ を定義し、候補ベクトル $\boldsymbol{w}=(\boldsymbol{d},\boldsymbol{y}^T)^T$ を作成する。

候補ベクトルの区間幅は大きすぎたり、小さすぎたりすると定理の検証が失敗する。そのためある程度小さい大きさを選ぶ必要がある。今回の実装では $\epsilon=10^{-9}$ とした。

In [4]:
using IntervalArithmetic
lam = λ[1]
x1 = x[:, 1] # 対応する固有ベクトル
x = x1 ./ sqrt(x1' * x1)
ysize = length(x)

ϵ = 1e-9 # size of candidate vector

if isreal(lam) && isreal(x)
    lam = real(lam)
    x = real(x)
    id = interval(0, ϵ; format=:midpoint)
    iy = interval.(zeros(ysize), ϵ; format=:midpoint)
    iI = interval(Matrix{Float64}(I, n + 1, n + 1))
    izero = interval(0)
else
    id = Complex(interval.(0, ϵ; format=:midpoint), interval.(0, ϵ; format=:midpoint))
    iy = Complex.(interval.(zeros(ysize), ϵ; format=:midpoint), interval.(zeros(ysize), ϵ; format=:midpoint))
    iI = interval(Matrix{Complex{Float64}}(I, n + 1, n + 1))
    izero = interval(0 + 0im)
end
iw = [id; iy]
Out[4]:
31-element Vector{Interval{Float64}}:
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
   ⋮
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
 [-1.00001e-09, 1.00001e-09]_com
  1. ヤコビ行列$DF$ と $z$、$g(\boldsymbol{w}):=z+\left(I-R\cdot DF\left(\bar\lambda + \boldsymbol{d}, \bar x + \boldsymbol{y}\right)\right)\boldsymbol{w}$ を区間演算により計算する。
In [5]:
ix = interval(x);
ilam = interval(lam);
iA = interval(A);
iB = interval(B);
iDF = [izero transpose(interval(2) * (ix + iy)); -iB*(ix+iy) iA-(ilam.+id)*iB]
R = inv(mid.(iDF));
iR = interval(R);
z = -iR * [dot(ix, ix) - interval(1); iA * ix - ilam * iB * ix]
gw = z + (iI - iR * iDF) * iw
Out[5]:
31-element Vector{Interval{Float64}}:
 [-1.53421e-12, 1.60364e-12]_com
 [-6.28847e-15, 5.87119e-15]_com
 [-6.2841e-15, 6.24573e-15]_com
 [-1.33737e-14, 1.35622e-14]_com
 [-5.34187e-15, 5.95928e-15]_com
 [-1.65682e-14, 1.55644e-14]_com
 [-4.44868e-15, 4.00426e-15]_com
 [-2.1826e-14, 2.10103e-14]_com
 [-5.66461e-15, 6.17318e-15]_com
 [-9.08141e-15, 8.89336e-15]_com
 [-1.68061e-14, 1.59398e-14]_com
 [-3.91177e-15, 3.74838e-15]_com
 [-3.26079e-15, 3.35143e-15]_com
   ⋮
 [-3.99718e-15, 4.32441e-15]_com
 [-9.42529e-15, 9.6057e-15]_com
 [-1.38155e-14, 1.34268e-14]_com
 [-1.58001e-14, 1.49484e-14]_com
 [-3.53401e-15, 3.64503e-15]_com
 [-1.14521e-14, 1.23865e-14]_com
 [-2.50687e-14, 2.35591e-14]_com
 [-1.37185e-14, 1.32005e-14]_com
 [-2.83772e-14, 2.69667e-14]_com
 [-2.06638e-14, 2.20965e-14]_com
 [-1.31275e-14, 1.3718e-14]_com
 [-1.19503e-14, 1.15808e-14]_com
  1. $g(\boldsymbol{w})\subsetneq \boldsymbol{w}$を満たすとき、$(\bar\lambda,\bar x)^T+g(\boldsymbol{w})$ の値を返す。
In [6]:
all(issubset_interval.(gw, iw))
Out[6]:
true
In [7]:
[ilam; ix] + gw
Out[7]:
31-element Vector{Interval{Float64}}:
 [-8.80724, -8.80723]_com
 [0.121617, 0.121618]_com
 [-0.130802, -0.130801]_com
 [-0.601561, -0.60156]_com
 [-0.162248, -0.162246]_com
 [-0.0338032, -0.0338031]_com
 [0.15269, 0.152691]_com
 [-0.100636, -0.100635]_com
 [-0.193449, -0.193448]_com
 [-0.209609, -0.209608]_com
 [-0.112811, -0.11281]_com
 [-0.0437788, -0.0437787]_com
 [-0.264069, -0.264068]_com
   ⋮
 [-0.00728467, -0.00728466]_com
 [-0.0881746, -0.0881745]_com
 [-0.187405, -0.187404]_com
 [0.00250188, 0.00250189]_com
 [0.13498, 0.134981]_com
 [-0.124575, -0.124574]_com
 [0.186956, 0.186958]_com
 [-0.178369, -0.178368]_com
 [-0.325139, -0.325138]_com
 [-0.205561, -0.205559]_com
 [0.146769, 0.14677]_com
 [-0.0235999, -0.0235998]_com

この第一成分は対象の固有値の厳密な包含である。

In [8]:
ilam + gw[1]
Out[8]:
[-8.80724, -8.80723]_com

この手順をverifyeig関数として定義する。

In [9]:
using IntervalArithmetic
function verifyeig(A, lam, x, B=B = Matrix(I,size(A)))
    x = x ./ sqrt(x' * x)
    ysize = length(x)

    ϵ = 1e-9 # size of candidate vector

    if isreal(lam) && isreal(x)
        lam = real(lam)
        x = real(x)
        id = interval(0, ϵ; format=:midpoint)
        iy = interval.(zeros(ysize), ϵ; format=:midpoint)
        iI = interval(Matrix{Float64}(I, ysize + 1, ysize + 1))
        izero = interval(0)
    else
        id = Complex(interval.(0, ϵ; format=:midpoint), interval.(0, ϵ; format=:midpoint))
        iy = Complex.(interval.(zeros(ysize), ϵ; format=:midpoint), interval.(zeros(ysize), ϵ; format=:midpoint))
        iI = interval(Matrix{Complex{Float64}}(I, ysize + 1, ysize + 1))
        izero = interval(0 + 0im)
    end
    iw = [id; iy]

    # DF(w) = [0 transpose(2*(x+w[2:end])) ; -B*(x+w[2:end]) A-(lam+w[1]).*B]
    ix = interval(x)
    ilam = interval(lam)
    iA = interval(A)
    iB = interval(B)

    iDF(w) = [izero transpose(interval(2) * (ix + w[2:end])); -iB*(ix+w[2:end]) iA-(ilam.+w[1])*iB]
    R = inv(mid.(iDF(zeros(ysize + 1))))
    iR = interval(R)
    z = -iR * [dot(ix, ix) - interval(1); iA * ix - ilam * iB * ix]
    g(w) = z + (iI - iR * iDF(w)) * w
    gw = g(iw)
    if all(issubset_interval.(gw, iw))
        while maximum(radius, gw) / norm([lam; x], 1) > 5e-13
            iw = gw
            gw = g(iw)
        end
        return ilam + gw[1]
    else
        return NaN
    end
end
Out[9]:
verifyeig (generic function with 2 methods)

部分固有対の精度保証方法であるので、全固有値を精度保証するには向かない。以下のように全固有値を精度保証するには $O(n^4)$ のオーダーが必要になる。全固有値を精度保証するには、例えば、標準固有値問題の精度保証付き数値解法のような方法が有効である。

In [10]:
using LinearAlgebra
n = 30
A, B = randn(n, n), randn(n, n)
λ, x = eigen(A, B)
Out[10]:
GeneralizedEigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
30-element Vector{ComplexF64}:
   -3.927248809102562 + 0.0im
   -2.342175403774986 + 0.0im
  -1.8204436211401438 + 0.0im
  -0.9825359650919164 + 0.9062641005675697im
  -0.9825359650919163 - 0.9062641005675698im
  -0.9369434243815118 + 0.0im
  -0.7678548617673059 + 0.0im
  -0.6336042505433659 - 4.4581246542775945im
  -0.6336042505433657 + 4.4581246542775945im
 -0.41284829356200076 + 0.0im
  -0.3207959087791363 + 0.0im
 -0.26666982104990433 - 0.36793157410745986im
 -0.26666982104990433 + 0.36793157410745986im
                      ⋮
 -0.04156652673607428 + 0.8154422472066382im
   0.2577455745317516 - 0.18455474115985465im
   0.2577455745317516 + 0.18455474115985465im
   0.6689537858084744 - 0.24860338555344522im
   0.6689537858084744 + 0.24860338555344522im
   0.7723122085046161 + 0.6492318019823247im
   0.7723122085046162 - 0.6492318019823247im
   1.1473309779693959 - 0.10819027220057749im
   1.1473309779693959 + 0.10819027220057749im
   1.7003772980863652 + 0.0im
   1.7723005862476457 - 2.661353859339048im
   1.7723005862476457 + 2.661353859339048im
vectors:
30×30 Matrix{ComplexF64}:
  -0.696811+0.0im   -0.920033+0.0im  …    0.867718+0.0608822im
  -0.637245+0.0im   -0.236322+0.0im        0.24356+0.228255im
  -0.630751+0.0im    -0.86956+0.0im       0.486714-0.0649732im
 0.00190353+0.0im   -0.102257+0.0im      -0.334326-0.375381im
 -0.0135138+0.0im    0.452453+0.0im      -0.128828+0.209189im
  0.0991277+0.0im    0.461829+0.0im  …   -0.179949-0.13841im
  0.0410025+0.0im   -0.674928+0.0im       0.408206+0.591794im
  0.0469751+0.0im    0.187833+0.0im      -0.276895+0.386155im
 0.00228834+0.0im   -0.833872+0.0im      0.0625323-0.136356im
   0.123516+0.0im  -0.0498425+0.0im      -0.090858+0.091819im
   -0.54912+0.0im    0.768362+0.0im  …    0.138197-0.169806im
   0.241703+0.0im   -0.568661+0.0im       0.172141+0.26653im
   0.317403+0.0im   -0.551095+0.0im       0.236936-0.0286702im
           ⋮                         ⋱  
        1.0+0.0im    0.742088+0.0im       -0.20981+0.116056im
  -0.189429+0.0im    0.174489+0.0im      -0.130892-0.10924im
  -0.010899+0.0im   0.0251151+0.0im  …    0.289061-0.191932im
   0.150996+0.0im   -0.278616+0.0im     -0.0761204-0.171493im
  -0.110771+0.0im         1.0+0.0im      0.0302184+0.598693im
  -0.229355+0.0im     0.44826+0.0im        0.44877+0.289908im
   0.852991+0.0im    0.715839+0.0im      -0.261972+0.551214im
   0.540259+0.0im    0.099366+0.0im  …   -0.718619-0.13813im
   -0.70162+0.0im   -0.413567+0.0im       0.349212+0.152074im
   0.133896+0.0im  -0.0737264+0.0im       0.358246-0.0728801im
   0.246352+0.0im     0.18205+0.0im     0.00756773+0.0964924im
   0.317812+0.0im    -0.35158+0.0im      -0.473799-0.141664im
In [11]:
# ilam = zeros(Interval,length(λ))
for i = 1:length(λ)
    @show ilam = verifyeig(A, λ[i], x[:, i], B)
end

ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-3.9272488091026454, -3.9272488091024935, com)
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-2.3421754037750544, -2.342175403774924, com)
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-1.8204436211401773, -1.820443621140108, com)
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.9825359650919452, -0.9825359650918856, com) + Interval{Float64}(0.9062641005675408, 0.9062641005675977, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.9825359650919446, -0.9825359650918858, com) + Interval{Float64}(-0.9062641005675972, -0.9062641005675407, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.9369434243815297, -0.9369434243814961, com)
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.7678548617673198, -0.7678548617672895, com)
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.6336042505435233, -0.6336042505432147, com) + Interval{Float64}(-4.458124654277739, -4.458124654277461, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.6336042505435259, -0.6336042505432139, com) + Interval{Float64}(4.458124654277459, 4.45812465427774, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.4128482935620064, -0.41284829356199426, com)
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.32079590877914593, -0.320795908779129, com)
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.26666982104991177, -0.2666698210498967, com) + Interval{Float64}(-0.3679315741074671, -0.3679315741074521, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.26666982104991177, -0.2666698210498967, com) + Interval{Float64}(0.3679315741074521, 0.3679315741074671, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.20595380010554623, -0.20595380010549347, com) + Interval{Float64}(1.0529535033506603, 1.0529535033507122, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.20595380010554612, -0.20595380010549386, com) + Interval{Float64}(-1.0529535033507118, -1.0529535033506605, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.09080177424957564, -0.0908017742495591, com) + Interval{Float64}(-0.18038183818252132, -0.1803818381825058, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.09080177424957567, -0.09080177424955913, com) + Interval{Float64}(0.1803818381825058, 0.1803818381825213, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.04156652673609379, -0.04156652673605269, com) + Interval{Float64}(-0.8154422472066599, -0.8154422472066175, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(-0.04156652673609379, -0.04156652673605269, com) + Interval{Float64}(0.8154422472066175, 0.8154422472066599, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(0.25774557453174884, 0.25774557453175456, com) + Interval{Float64}(-0.18455474115985773, -0.18455474115985138, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(0.25774557453174884, 0.25774557453175456, com) + Interval{Float64}(0.18455474115985138, 0.18455474115985773, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(0.668953785808457, 0.6689537858084947, com) + Interval{Float64}(-0.24860338555346076, -0.24860338555342834, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(0.668953785808457, 0.6689537858084947, com) + Interval{Float64}(0.24860338555342834, 0.24860338555346076, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(0.7723122085045907, 0.7723122085046473, com) + Interval{Float64}(0.6492318019822947, 0.6492318019823529, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(0.7723122085045908, 0.7723122085046472, com) + Interval{Float64}(-0.6492318019823528, -0.6492318019822946, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(1.1473309779693606, 1.1473309779694252, com) + Interval{Float64}(-0.10819027220060964, -0.10819027220054633, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(1.1473309779693606, 1.1473309779694252, com) + Interval{Float64}(0.10819027220054633, 0.10819027220060964, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(1.7003772980863459, 1.7003772980864003, com)
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(1.7723005862475718, 1.77230058624774, com) + Interval{Float64}(-2.6613538593391373, -2.661353859338972, com)im
ilam = verifyeig(A, λ[i], x[:, i], B) = Interval{Float64}(1.7723005862475718, 1.77230058624774, com) + Interval{Float64}(2.661353859338972, 2.6613538593391373, com)im

TODO　固有値が縮退し多重固有値である場合や複数の固有値がクラスター形成している場合、この方法は失敗する。この場合の精度保証付き数値計算方法は文献3にある。今後この方法を実装していきたい。

本資料は以下のような文献・Web ページ等を参考に書いている。

参考文献¶

  1. 大石進一編著, 精度保証付き数値計算の基礎, コロナ社, 2018.
    （精度保証付き数値計算の教科書、3.4.2章にある「非線形方程式を利用した精度保証法」を実装した）

  2. S.M. Rump. Guaranteed Inclusions for the Complex Generalized Eigenproblem. Computing, 42:225-238, 1989.
    （今回の方法を初めて発表した原著論文）

  3. S.M. Rump. Computational Error Bounds for Multiple or Nearly Multiple Eigenvalues. Linear Algebra and its Applications (LAA), 324:209–226, 2001.
    （多重固有値や多数の固有値のクラスターに対する精度保証付き数値計算方法。INTLABのverifyeig.mにはこの方法も実装されている）

近藤慎佑, 高安亮紀，2022年9月19日（最終更新：2024年2月5日）