標準固有値問題の精度保証付き数値解法

固有値問題とは

固有値問題は与えられた行列の固有値および固有ベクトルを求めること。つまり、 $ \newcommand{\C}{\mathbb{C}} \newcommand{\IC}{\mathbb{IC}} \newcommand{\bX}{\bar{X}} $

$$ Ax=\lambda x,\quad A\in \C^{n\times n} $$

を満たす $\lambda\in \C$, $x\in \C^{n}$ ($x\neq 0$) を求める問題を指す。

Juliaで固有値と固有ベクトルを求める際には、LinearAlgebraパッケージに含まれている関数eigen・eigvals・eigvecsで求めることができる。

関数 役割
eigen() 固有値・固有ベクトル両方
eigvals() 固有値
eigvecs() 固有ベクトル

これらは行列の構造・サイズをみて適宜、Lapackの関数を呼び出している（はず）。

密行列の全固有値の精度保証

密行列とは、行列内のほとんどの要素が $0$ でない行列のことを指す。ここでは、行列$A\in\C^{n\times n}$ $(i=1,2,\cdots,n)$ の全ての固有値に対する精度保証アルゴリズムを考える。

実装したいアルゴリズム

(1) ある行列 $A$ の全ての近似固有ベクトル $\bX\in\C^{n\times n}$ $(i=1,2,\cdots,n)$ を求める。

(2) $A\cdot\bX$ の包含 $C\in\IC^{n\times n}$ を求める。

(3) 区間連立1次方程式 $\bX G=C$ の解集合の包含 $G\supset \bX^{-1}\cdot C$ ($G\in\IC^{n\times n}$) を求める。

(4) $G\in\IC^{n\times n}$ に対してゲルシュゴリンの定理を適用することを考える。

(1) 行列 $A$ の全ての近似固有ベクトル $\bX\in\C^{n\times n}$ $(i=1,2,\cdots,n)$ を求める。

行列 $A\in\C^{n\times n}$ $(i=1,2,\cdots,n)$ が対角化可能であるとするとき、$\lambda_i$ に対する固有ベクトル$x^{(i)}\in\C^{n}$ $(i=1,2,\cdots,n)$ を並べた行列

$$ X :=[x^{(1)},x^{(2)},\cdots,x^{(n)}]\in\C^{n\times n} $$

によって、$A$ は

$$ X^{-1}AX = D = \rm{diag}(\lambda_1,\lambda_2,\cdots,\lambda_n) $$

と対角化される。

しかし、実際のところ上記の数値計算で得られる固有ベクトル $X$ は近似の値 $\bX\in \C^{n\times n}$ (以下、近似固有ベクトルという) である。しかし、近似固有ベクトルを並べた行列を使って、対角化したものについては、対角優位性が期待できる。この行列を $G$ とすると、

$$ G := \bX^{-1}A\bX $$

となり、行列 $G$ の対角成分は行列 $A$ の固有値の近似となる。

また、対角優位性とは、ある行列$A=(a_{ij})$について、

$$ \left|a_{i i}\right|>\sum_{j \neq i}\left|a_{i j}\right| \quad(i=1,2, \ldots, n) $$

であること。つまり、行列のどの対角成分の絶対値も、それと同じ列にある非対角成分の和よりも大きい時。

(2) $A\cdot\bX$ の包含 $C\in\IC^{n\times n}$ を求める。

方法1各演算を区間演算に変更した方式(参考：interval_dot-mul.ipynb)、行列のサイズが小さい場合はこちらの方が早い（たぶん）

方法2BLASを使う。IntervalFunctions.jlの呼び出し。大きな行列に対してはこの方が早い。

(3 )区間連立1次方程式 $\hat XG=C$ の解集合の包含 $G\supset\hat X^{-1}\cdot C$ を求める。

4) $G\in\IC^{n\times n}$ に対してゲルシュゴリンの定理を適用することを考える。

ゲルシュゴリンの包含定理

行列 $A=(A_{ij})\in \C^{n\times n}$ について、$A$ の全ての固有値 $\lambda_i$ は、$\bigcup_{1\le i \le n} U_i$ の内部に存在する。

$$ A\subseteq\bigcup_{1\le i \le n} U_i $$

ただし、

$$ U_i = \left\{z\in\mathbb{C}:|z-a_{ii}|\le \sum_{j\neq i}|a_{ij}|\right\} $$

である。円盤領域 $U_i$ をゲルシュゴリン円板といい、$A$ が強い優対角性を持つとき、$A$ の対角成分が $A$ の固有値の良い近似となる。今回の場合、$G=(g_{ij})_{1\le i,j\le n}\in\IC^{n\times n}$ の対角成分 $(g_{ii})_{1\le i\le n}$ を $(c_{i})_{1\le i\le n}\in \IC^{n}$ とし、各成分の半径を $\mathrm{rad}(c_i)$ とすると、

$$ \begin{array}{ll} \mbox{ゲルシュゴリン円の中心:}&\mathrm{mid}(c_i)\\ \mbox{ゲルシュゴリン円の半径:}&r_i=\sum_{j\neq i}\rm{mag}(g_{ij})+\rm{rad}(c_{i})\\ \end{array} $$

となる。

以上により、全固有値（ $n=100$ ）の精度保証がゲルシュゴリンの包含定理と区間演算によって成功した。

TODO　標準・一般化固有値問題を対象に非線形方程式の精度保証法を利用することで、一部の固有値を精度保証する方法（verifyeig）を今後、紹介したい。

本資料は以下のような文献・Web ページ等を参考にこの文章は書いています.

参考文献

  1. 大石進一編著, 精度保証付き数値計算の基礎, コロナ社, 2018.
    精度保証付き数値計算の教科書、3.4.1章にある「密行列に対する精度保証法」を実装した。

  2. 齊藤宣一, 計算数理演習, (最終閲覧日 2020年12月9日).
    2019年度に行われた、計算数理演習（東京大学理学部数学科，教養学部統合自然科学科）の講義資料．ゲルシュゴリンの包含定理について、MATLABで実装した例が掲載されている.

井藤佳奈子, 高安亮紀, 2020年12月9日, 2022年2月21日（最終更新）