この内容はNumerical Linear Algebra内の以下のLectureに相当します．

  • Lecture 33, The Arnoldi Iteration
  • Lecture 34, How Arnoldi Locates Eigenvalues

しかし，この本にはどのように多項式$p(z)$を求めるかが(明示的に)書いてありません．
また，自分のググり力では，多項式の求め方を見つけることができませんでした．

考えてみると，こうすれば多項式が求まるのでは，という方法がでました．
これが正しいのかは分かりませんが，結果をプロットしてみると，確かに性質を満たしていました．
もっとスマートな方法があるに違いないと思いますが，情報をまとめようと思いblogに書きました．

Arnoldi法

Arnoldi法とは，Krylov部分空間$\mathcal{K}(A, b)$に対して正規直交基底を求める方法です．
これは単純に$\mathcal{K}(A, b)$に対して，Modified Gram-Schmidtを適用し，直交基底を求める際の情報をHessenberg行列に格納するアルゴリズムです．

次の記号を導入します．$Q$は$n\times n$の正規直交基底です．$H$は$(n+1) \times n$のHessenberg行列で，上三角行列に下1行を追加した形です．

$$ Q_n = \left[ \begin{array}{c|c|c|c} & & & \\
& & & \\
q_1 & q_2 & \cdots & q_n \\
& & & \\
& & & \\
\end{array} \right], \tilde{H}_n = \left[\begin{array}{cccc} h_{11} & h_{12} & \cdots & h_{1n}\\
h_{21} & h_{22} & \cdots & h_{2n}\\
& \ddots & \ddots & \vdots\\
& & h_{n,n-1} & h_{n,n}\\
& & & h_{n+1,n} \end{array}\right] $$

Arnoldi法は以下の関係式を満たす，$Q_{n}, \tilde{H}_n$を求めます．

$$ A Q_n = Q_{n+1} \tilde{H}_n $$

具体的にどのように求めていくのか，それは次のコードのとおりです．

Q[:, 0] = b / la.norm(b)
for k in range(it):
    u = A.dot(Q[:, k])
    for j in range(k+1):
        H[j,k] = np.dot(Q[:, j], u)
        u = u - H[j,k] * Q[:, j]
    
    H[k+1, k] = la.norm(u)
    if H[k+1, k] <= EPS
        break
    Q[:, k+1] = u/H[k+1, k]

現在選ばれている全ての直交基底との射影成分を引くことで，新しい直交基底を選んでいます．
また基底$Q_j$方向のpowerをHessenberg行列の$H[j, k]$に格納しています．
$H[k+1, k]$に格納されるのは，現在の基底で表すことができなかったベクトルのノルムです．
つまり$H[k+1, k]$が$0$になれば，現在の基底で表す事が可能なので，新しい基底を作らずにループを抜けます．

射影での解釈

$Q^{*}_{n} Q_{n+1}$という行列を考えると，これは対角に$1$が並ぶ$n \times (n+1)$の行列です．
なぜなら$Q$は直交基底なので，$q_i$と$q_j$の内積は，$i == j$の時に$1$(正規)，それ以外の時に$0$(直交)だからです．
そのため$Q^{*}_{n} Q_{n+1}\hat{H}_{n}$は，$n+1$行目を削除した$n \times n$の行列になります．

$$ Q^*_n Q_{n+1} \tilde{H}_n = H_n = \left[\begin{array}{cccc} h_{11} & h_{12} & \cdots & h_{1n}\\
h_{21} & h_{22} & \cdots & h_{2n}\\
& \ddots & \ddots & \vdots\\
& & h_{n,n-1} & h_{n,n}\\
\end{array}\right] $$

先に述べた式に戻り，左から$Q^{*}_{n}$をかける．

$$ A Q_n = Q_{n+1} \tilde{H}_n \\
Q^{*}_{n} A Q_n = Q^{*}_n Q_{n+1} \tilde{H}_n\\
Q^{*}_{n} A Q_n = H_n \hspace{3.0em} $$

$Q_n$は$\mathcal{K}_n$の正規直交基底なので，$H_n$は$\mathcal{K}_n$上に$A$を射影した行列となります．
そのため$H_n$の固有値は，$A$の固有値と関係があります．
$H_n$の固有値は，$\textit{Arnoldi eigenvalue estimates}$，または$\textit{Ritz values}$と呼ばれ，$A$のスペクトル半径に近い固有値を近似します．

これを確かめるために，$n=25, 40$の時に$n \times n$の実数ランダム行列$A$を作成し，Arnoldi法を適用してみました．
右の図の青のxが$A$の固有値で，オレンジのoが$H_n$の固有値です．
左の図は$A$をimshowで出した結果です．
こうみると，徐々に固有値を近似していく($A$の固有値の上に収まる)様子が分かります．

Arnoldi法と多項式近似

Kyrolov部分空間$\mathcal{K}$上の任意のベクトルは次式で表せます．

$$ x = c_0 b + c_1 Ab + c_2 A^2 b + \cdots + c_{n-1} A^{n-1}b $$

ここで多項式$q(z)$を導入します．すると$x$は次のように書けます．

$$ q(z) = c_0 + c_1 z + \cdots + c_{n-1} z^{n-1}\\
x = q(A) b $$

$q(z)$の中で$P^n$をmonicな$n$次多項式($n$次の項の係数が$1$)とした時，Arnoldi法は次の問題を解いていると見れます．

$$ \text{Find } p^n \in P^n, \text{such that } \|p^n (A) b\| = minimum $$

また$p^n \in P^n$の場合，$p^n(A)b = A^n b - Q_n y$となるので，次のような言い換えができます．

$$ \text{Find } y, \|A^n b - Q_n y\| = minimum $$

これを幾何学的に理解するには，p.260に載っている図が分かりやすいです．(引用 p.260 Figure 34.1 The least squares polynomial approximation problem underlying the Arnoldi iteration.)

$A^n b$を見る時は，既に$n$本の直交基底$Q_n$で表現されるKyrolov部分空間$\mathcal{K}_n$があり，Krylov列/基底の$n+1$本目を見ている場合です．
つまり，この多項式$p(A)b$は，現在の基底$Q_n$で表せないベクトルを表しています．
そしてArnoldi法は，ベクトルの最小化を行ないます．

Arnoldi Lemniscates

Arnoldi法の収束過程は，複素平面に$\textit{lemniscates}$を描くことで幾何学的に見ることができます．
これは$\textit{Arnoldi lemniscates}$と呼ばれ，次の式で描くことができます．

$$ \{ z \in \mathbb{C} : |p(z)| = C\} \\
C = \frac{\|p^n (A) b \|}{\|b\|} $$

さてここで本題となりますが，どのようにすれば$p(z)$を導出できるのでしょうか．
$p(z)$が分かるというのは，係数$c_0, c_1, …$が分かっているということです．
この係数を求めるために，以下の関係を列形式で書いて，$n=1, 2$の場合を見ていきます．

$$ A Q_n = Q_{n+1} \tilde{H}_n $$ $$ A \left[ \begin{array}{c|c|c|c} & & & \\
& & & \\
q_1 & q_2 & \cdots & q_n \\
& & & \\
& & & \\
\end{array} \right] = \left[ \begin{array}{c|c|c|c|c} & & & &\\
& & & &\\
q_1 & q_2 & \cdots & q_n & q_{n+1} \\
& & & &\\
& & & &\\
\end{array} \right] \left[\begin{array}{cccc} h_{11} & h_{12} & \cdots & h_{1n}\\
h_{21} & h_{22} & \cdots & h_{2n}\\
& \ddots & \ddots & \vdots\\
& & h_{n,n-1} & h_{n,n}\\
& & & h_{n+1,n} \end{array}\right] $$

  • $n=1$の時，$q_1 = \frac{b}{\|b\|}$を左辺にだけ代入します． $$ A q_1 = h_{11} q_1 + h_{21} q_2 \\
    A \frac{b}{\|b\|} = h_{11} q_1 + h_{21} q_{2} \\
    A b - \|b\|h_{11} q_{1} = \|b\|h_{21} q_{2} $$

    • これは$y_1 = \left[\|b\| h_{11} \right]$と見れば左辺は $$ A b - \|b\|h_{11} q_{1} = Ab - Q_1 y_1 $$
    • となり，「$p^n \in P^n$の場合，$p^n(A)b = A^n b - Q_n y$」を用いれば，$Ab$の係数は$1$なので，右辺が$p^1(A)b$となります．
      lemniscatesを書く際に求める$C$は，$\|p^n(A)b\|$を$\|b\|$で割るので，自然な感じがします． $$ \|b\|h_{21} q_{2} = p^1(A)b $$
    • この右辺が$p^1(A)b$を表すことが分かったので，$c_0 b + c_1 Ab + \cdots + c_{n-1} A^{n-1}$の形にします． $$ A b - \|b\|h_{11} \frac{b}{\|b\|} = \|b\|h_{21} q_{2} \\
      p^1(A)b = - h_{11} b + A b $$
    • $c_0 = -h_{11}$，$c_1 = 1$，$C = h_{21}$となり，これで$p(z) = C$を満たすlemniscatesが書けます．

  • $n=2$の時も同様，$q_{1} = \frac{b}{\|b\|}$，また上記の式から$q_2 = -\frac{h_{11}}{h_{21}\|b\|}b + \frac{1}{h_{21}\|b\|}Ab$を代入します． $$ A q_2 = h_{12} q_1 + h_{22} q_{2} + h_{32} q_{3} \\
    A \left( -\frac{h_{11}}{h_{21}\|b\|}b + \frac{1}{h_{21}\|b\|}Ab \right) = h_{12} \frac{b}{\|b\|} + h_{22} \left( -\frac{h_{11}}{h_{21}\|b\|}b + \frac{1}{h_{21}\|b\|}Ab \right) + h_{32} q_3 \\
    -\frac{h_{11}}{h_{21} \|b\|} Ab + \frac{1}{h_{21}\|b\|} A^2 b = \frac{h_{12}}{\|b\|}b - \frac{h_{11} h_{22}}{h_{21} \|b\|} b + \frac{h_{22}}{h_{21}\|b\|} Ab + h_{32} q_3 \\
    -h_{11} Ab + A^2b = h_{12}h_{21}b - h_{11} h_{22} b + h_{22} Ab + \|b\|h_{21}h_{32} q_3\\
    \|b\| h_{21}h_{32} q_3 = \left(-h_{12}h_{32} + h_{11}h_{22}\right)b + \left(-h_{11} + h_{22}\right) Ab + A^2 b $$

    • $c_0 = -h_{12} h_{32} + h_{11} h_{22}$，$c_{1} = -h_{11} + h_{22}$，$c_2 = 1$，$C = h_{21} h_{32}$

これを手順的に行って，プロットしたものが以下の図になります．
複素平面にメッシュを作り，$p(z)$の値が$C$以下の場合の等高線を出しました．
青が濃くなるほど値が小さく，一番濃い青に囲まれた箇所に近似固有値があることから
多項式の解がHessenbergの固有値になる，という性質を満たしているが分かります．

まとめ

間違っている記述がある場合は，申し訳ありません．
勉強をするので，是非とも教えていただけると嬉しいです．

多項式の導出には，もっとスマートな方法があるはずです．
もっと理解が深められるように進めていきます．