連立一次方程式 A x = b A\boldsymbol{x}=\boldsymbol{b} A x = b
GMRES は x 0 + K k ( A , r 0 ) \boldsymbol{x}_0+\mathcal{K}_k(A,\boldsymbol{r}_0) x 0 + K k ( A , r 0 ) ∥ b − A x ∥ 2 \|\boldsymbol{b}-A\boldsymbol{x}\|_2 ∥ b − A x ∥ 2 ( k + 1 ) × k (k+1)\times k ( k + 1 ) × k Givens 回転 で 1 反復あたり O ( k ) O(k) O ( k ) 厳密な算術では、Arnoldi の過程が止まった時点で厳密解に達し、n n n
残差は r k = p k ( A ) r 0 \boldsymbol{r}_k=p_k(A)\boldsymbol{r}_0 r k = p k ( A ) r 0 p k ( 0 ) = 1 p_k(0)=1 p k ( 0 ) = 1 正規行列 では ∥ r k ∥ / ∥ r 0 ∥ ≤ min p max λ ∈ σ ( A ) ∣ p ( λ ) ∣ \|\boldsymbol{r}_k\|/\|\boldsymbol{r}_0\|\le\min_p\max_{\lambda\in\sigma(A)}|p(\lambda)| ∥ r k ∥/∥ r 0 ∥ ≤ min p max λ ∈ σ ( A ) ∣ p ( λ ) ∣ 対角化可能 な行列ではこれに κ ( V ) \kappa(V) κ ( V ) ∣ z − c ∣ ≤ ρ < ∣ c ∣ |z-c|\le\rho<|c| ∣ z − c ∣ ≤ ρ < ∣ c ∣ ( ρ / ∣ c ∣ ) k κ ( V ) (\rho/|c|)^k\kappa(V) ( ρ /∣ c ∣ ) k κ ( V ) κ ( V ) \kappa(V) κ ( V )
右前処理 A M − 1 u = b AM^{-1}\boldsymbol{u}=\boldsymbol{b} A M − 1 u = b 元の系の残差 を最小にします。対称な M 行列（Stieltjes 行列。非対角成分がすべて 0 以下の正定値な対称行列）では、非零パターンを固定した不完全 Cholesky 分解 IC(0)（Cholesky 分解の途中で新しくできる非零成分を捨てる近似）が破綻しません。証明の鍵は Schur 補行列（1 段ぶん消去したあとに残る小さな行列）が Stieltjes 行列であり続けることと、非負行列のスペクトル半径 の単調性です。
Krylov 部分空間 : r 0 \boldsymbol{r}_0 r 0 A A A Krylov 部分空間 Arnoldi の過程 : Krylov 部分空間の正規直交基底を順に作る手順。→ Arnoldi の過程 Givens 回転 : 2 つの行だけを混ぜて、その中の 1 つの成分を 0 にする直交変換。→ Givens 回転 正規行列 : ユニタリ行列で対角化できる行列。→ 正規行列 条件数 : 行列のノルムと逆行列のノルムの積で、解の揺れやすさを測る量。→ 条件数 スペクトル半径 : 固有値の絶対値の最大値。→ スペクトル半径 偏微分方程式を離散化すると、未知数が 10 5 10^5 1 0 5 10 8 10^8 1 0 8 Cholesky 分解と疎行列の直接法 ）は fill-in（分解の途中で、もとは零だった成分が非零になること）のため 3 次元の問題では記憶容量が足りなくなり、反復法が必要になります。対称正定値なら共役勾配法（Krylov 部分空間法 I ）が標準ですが、移流のある流れの問題・非対称な離散化・固有値問題のシフト反転では行列が非対称になり、共役勾配法の最適性は使えません。
GMRES は「Krylov 部分空間の上で残差を最小にする」という最も素直な一般化です。その代わり、Lanczos の三項漸化式が使えないので、反復ごとに全部の基底ベクトルと直交化する必要があり、記憶と演算が反復回数とともに増えます。そこで実務では、(1) 反復回数そのものを減らす前処理 と、(2) 一定回数で打ち切る再始動 が欠かせません。この記事では GMRES の最適性と収束の理論、前処理の 2 つの置き方の違い、代表的な前処理である不完全 Cholesky 分解がいつ安全に計算できるかを証明します。
A ∈ C n × n A\in\mathbb{C}^{n\times n} A ∈ C n × n b ∈ C n \boldsymbol{b}\in\mathbb{C}^n b ∈ C n x 0 \boldsymbol{x}_0 x 0 r 0 = b − A x 0 ≠ 0 \boldsymbol{r}_0=\boldsymbol{b}-A\boldsymbol{x}_0\ne\boldsymbol{0} r 0 = b − A x 0 = 0 β = ∥ r 0 ∥ \beta=\|\boldsymbol{r}_0\| β = ∥ r 0 ∥ P k \mathcal{P}_k P k k k k K k ( A , r 0 ) \mathcal{K}_k(A,\boldsymbol{r}_0) K k ( A , r 0 ) r 0 \boldsymbol{r}_0 r 0 d d d Q k Q_k Q k H k H_k H k H ~ k \tilde H_k H ~ k Krylov 部分空間法 I の Krylov 部分空間の定義（定義 2.1）［Krylov 部分空間法 I］ と Arnoldi の過程（定義 3.1）［Krylov 部分空間法 I］ の記号（b \boldsymbol{b} b r 0 \boldsymbol{r}_0 r 0
定義 3.1 （GMRES ） ↗
1 ≤ k ≤ d 1\le k\le d 1 ≤ k ≤ d x 0 + K k ( A , r 0 ) \boldsymbol{x}_0+\mathcal{K}_k(A,\boldsymbol{r}_0) x 0 + K k ( A , r 0 ) ∥ b − A x ∥ \|\boldsymbol{b}-A\boldsymbol{x}\| ∥ b − A x ∥ x k \boldsymbol{x}_k x k k k k r k = b − A x k \boldsymbol{r}_k=\boldsymbol{b}-A\boldsymbol{x}_k r k = b − A x k
次の定理のねらいは 2 つです。GMRES の近似解が、Arnoldi の関係式を通して小さな最小二乗問題に置き換わること、そしてその問題が Givens 回転で 1 反復ずつ安く更新できることです。Givens 回転は 2 つの行だけを混ぜて 1 つの成分を 0 にする操作で、反復が進んで新しい列が 1 本増えたときも、その列に回転をかけ直して新しい回転を 1 つ足せば済みます。
定理 3.2 （GMRES の最小残差性と実装 ） ✓
1 ≤ k < d 1\le k<d 1 ≤ k < d
H ~ k \tilde H_k H ~ k ∥ β e 1 − H ~ k y ∥ \|\beta\boldsymbol{e}_1-\tilde H_k\boldsymbol{y}\| ∥ β e 1 − H ~ k y ∥ y k ∈ C k \boldsymbol{y}_k\in\mathbb{C}^k y k ∈ C k x k \boldsymbol{x}_k x k x k = x 0 + Q k y k \boldsymbol{x}_k=\boldsymbol{x}_0+Q_k\boldsymbol{y}_k x k = x 0 + Q k y k ∥ r k ∥ = ∥ β e 1 − H ~ k y k ∥ \|\boldsymbol{r}_k\|=\|\beta\boldsymbol{e}_1-\tilde H_k\boldsymbol{y}_k\| ∥ r k ∥ = ∥ β e 1 − H ~ k y k ∥ A , b , x 0 A,\boldsymbol{b},\boldsymbol{x}_0 A , b , x 0 j j j j + 1 j+1 j + 1 G j G_j G j j = 1 , … , k j=1,\ldots,k j = 1 , … , k G k ⋯ G 1 H ~ k = ( R k 0 T ) G_k\cdots G_1\tilde H_k=\binom{R_k}{\boldsymbol{0}^{\mathsf{T}}} G k ⋯ G 1 H ~ k = ( 0 T R k ) R k R_k R k k × k k\times k k × k g = G k ⋯ G 1 ( β e 1 ) \boldsymbol{g}=G_k\cdots G_1(\beta\boldsymbol{e}_1) g = G k ⋯ G 1 ( β e 1 ) y k = R k − 1 ( g 1 , … , g k ) T \boldsymbol{y}_k=R_k^{-1}(g_1,\ldots,g_k)^{\mathsf{T}} y k = R k − 1 ( g 1 , … , g k ) T ∥ r k ∥ = ∣ g k + 1 ∣ \|\boldsymbol{r}_k\|=|g_{k+1}| ∥ r k ∥ = ∣ g k + 1 ∣ k k k k + 1 k+1 k + 1 H ~ k + 1 \tilde H_{k+1} H ~ k + 1 G 1 , … , G k G_1,\ldots,G_k G 1 , … , G k G k + 1 G_{k+1} G k + 1 証明（定理 3.2） (1) Arnoldi の関係式（定理 3.2）［Krylov 部分空間法 I］ の 1 により、x ∈ x 0 + K k \boldsymbol{x}\in\boldsymbol{x}_0+\mathcal{K}_k x ∈ x 0 + K k x = x 0 + Q k y \boldsymbol{x}=\boldsymbol{x}_0+Q_k\boldsymbol{y} x = x 0 + Q k y y ∈ C k \boldsymbol{y}\in\mathbb{C}^k y ∈ C k Q k Q_k Q k y \boldsymbol{y} y x \boldsymbol{x} x r 0 = β q 1 = Q k + 1 ( β e 1 ) \boldsymbol{r}_0=\beta\boldsymbol{q}_1=Q_{k+1}(\beta\boldsymbol{e}_1) r 0 = β q 1 = Q k + 1 ( β e 1 )
b − A x = r 0 − A Q k y = Q k + 1 ( β e 1 − H ~ k y ) \boldsymbol{b}-A\boldsymbol{x}=\boldsymbol{r}_0-AQ_k\boldsymbol{y}=Q_{k+1}\bigl(\beta\boldsymbol{e}_1-\tilde H_k\boldsymbol{y}\bigr) b − A x = r 0 − A Q k y = Q k + 1 ( β e 1 − H ~ k y ) で、Q k + 1 Q_{k+1} Q k + 1 ∥ b − A x ∥ = ∥ β e 1 − H ~ k y ∥ \|\boldsymbol{b}-A\boldsymbol{x}\|=\|\beta\boldsymbol{e}_1-\tilde H_k\boldsymbol{y}\| ∥ b − A x ∥ = ∥ β e 1 − H ~ k y ∥ H ~ k \tilde H_k H ~ k j ≤ k < d j\le k<d j ≤ k < d h j + 1 , j > 0 h_{j+1,j}>0 h j + 1 , j > 0 H ~ k \tilde H_k H ~ k k + 1 k+1 k + 1 h j + 1 , j ≠ 0 h_{j+1,j}\ne0 h j + 1 , j = 0 H ~ k \tilde H_k H ~ k 最小二乗解の特徴づけ（定理 3.2）［最小二乗問題］ の 3 により最小解 y k \boldsymbol{y}_k y k x k = x 0 + Q k y k \boldsymbol{x}_k=\boldsymbol{x}_0+Q_k\boldsymbol{y}_k x k = x 0 + Q k y k
(2) j = 1 , … , k j=1,\ldots,k j = 1 , … , k j j j ( j , j ) (j,j) ( j , j ) a a a ( j + 1 , j ) (j+1,j) ( j + 1 , j ) b b b Givens 回転の性質（命題 3.4）［QR 分解］ の回転を作って第 j j j j + 1 j+1 j + 1 ( j + 1 , j ) (j+1,j) ( j + 1 , j ) b b b h j + 1 , j > 0 h_{j+1,j}>0 h j + 1 , j > 0 j + 1 j+1 j + 1 ( a , b ) ≠ ( 0 , 0 ) (a,b)\ne(0,0) ( a , b ) = ( 0 , 0 ) j j j j j j j + 1 j+1 j + 1 H ~ k = Q ( R k 0 T ) \tilde H_k=Q\binom{R_k}{\boldsymbol{0}^{\mathsf{T}}} H ~ k = Q ( 0 T R k ) Q T = G k ⋯ G 1 Q^{\mathsf{T}}=G_k\cdots G_1 Q T = G k ⋯ G 1 QR 分解による最小二乗解（定理 5.1）［最小二乗問題］ により R k R_k R k y k = R k − 1 ( g 1 , … , g k ) T \boldsymbol{y}_k=R_k^{-1}(g_1,\ldots,g_k)^{\mathsf{T}} y k = R k − 1 ( g 1 , … , g k ) T ∣ g k + 1 ∣ |g_{k+1}| ∣ g k + 1 ∣
H ~ k + 1 \tilde H_{k+1} H ~ k + 1 H ~ k \tilde H_k H ~ k ( h 1 , k + 1 , … , h k + 2 , k + 1 ) T (h_{1,k+1},\ldots,h_{k+2,k+1})^{\mathsf{T}} ( h 1 , k + 1 , … , h k + 2 , k + 1 ) T G 1 , … , G k G_1,\ldots,G_k G 1 , … , G k 1 1 1 k + 1 k+1 k + 1 k k k ( R k 0 ) \binom{R_k}{\boldsymbol{0}} ( 0 R k ) k + 2 k+2 k + 2 G 1 , … , G k G_1,\ldots,G_k G 1 , … , G k ( k + 1 , k + 1 ) (k+1,k+1) ( k + 1 , k + 1 ) ( k + 2 , k + 1 ) (k+2,k+1) ( k + 2 , k + 1 ) G k + 1 G_{k+1} G k + 1 g \boldsymbol{g} g ( g 1 , … , g k + 1 , 0 ) T (g_1,\ldots,g_{k+1},0)^{\mathsf{T}} ( g 1 , … , g k + 1 , 0 ) T G k + 1 G_{k+1} G k + 1
∎ 実装の要点. 反復 k k k q 1 , … , q k \boldsymbol{q}_1,\ldots,\boldsymbol{q}_k q 1 , … , q k k k k O ( n k ) O(nk) O ( nk ) k k k O ( n k 2 ) O(nk^2) O ( n k 2 ) O ( n k ) O(nk) O ( nk ) ∣ g k + 1 ∣ |g_{k+1}| ∣ g k + 1 ∣ x k \boldsymbol{x}_k x k x k \boldsymbol{x}_k x k QR 分解 の 再直交化の注意（注意 6.5）［QR 分解］ ）。
命題 3.3 （Arnoldi が止まれば厳密解 ） ✓
k = d k=d k = d h d + 1 , d = 0 h_{d+1,d}=0 h d + 1 , d = 0 H d H_d H d x d = x 0 + Q d H d − 1 ( β e 1 ) \boldsymbol{x}_d=\boldsymbol{x}_0+Q_dH_d^{-1}(\beta\boldsymbol{e}_1) x d = x 0 + Q d H d − 1 ( β e 1 ) A x = b A\boldsymbol{x}=\boldsymbol{b} A x = b x 0 + K d \boldsymbol{x}_0+\mathcal{K}_d x 0 + K d r d = 0 \boldsymbol{r}_d=\boldsymbol{0} r d = 0 d ≤ n d\le n d ≤ n A A A
証明（命題 3.3） Arnoldi の関係式（定理 3.2）［Krylov 部分空間法 I］ の 2 から A Q d = Q d H d AQ_d=Q_dH_d A Q d = Q d H d H d y = 0 H_d\boldsymbol{y}=\boldsymbol{0} H d y = 0 A Q d y = 0 AQ_d\boldsymbol{y}=\boldsymbol{0} A Q d y = 0 A A A Q d y = 0 Q_d\boldsymbol{y}=\boldsymbol{0} Q d y = 0 Q d Q_d Q d y = 0 \boldsymbol{y}=\boldsymbol{0} y = 0 H d H_d H d x = x 0 + Q d y \boldsymbol{x}=\boldsymbol{x}_0+Q_d\boldsymbol{y} x = x 0 + Q d y b − A x = Q d ( β e 1 − H d y ) \boldsymbol{b}-A\boldsymbol{x}=Q_d(\beta\boldsymbol{e}_1-H_d\boldsymbol{y}) b − A x = Q d ( β e 1 − H d y ) y = H d − 1 β e 1 \boldsymbol{y}=H_d^{-1}\beta\boldsymbol{e}_1 y = H d − 1 β e 1 d ≤ n d\le n d ≤ n d d d Krylov 部分空間の基本性質（命題 2.2）［Krylov 部分空間法 I］ の 2・3 です。
∎ 止まること（breakdown）は失敗ではなく、成功です。共役勾配法と違い、GMRES の残差は単調非増加 です（K k ⊂ K k + 1 \mathcal{K}_k\subset\mathcal{K}_{k+1} K k ⊂ K k + 1
補題 4.1 （残差多項式 ） ✓
1 ≤ k ≤ d 1\le k\le d 1 ≤ k ≤ d
∥ r k ∥ = min { ∥ p ( A ) r 0 ∥ : p ∈ P k , p ( 0 ) = 1 } . \|\boldsymbol{r}_k\|=\min\bigl\{\|p(A)\boldsymbol{r}_0\|:p\in\mathcal{P}_k,\ p(0)=1\bigr\}. ∥ r k ∥ = min { ∥ p ( A ) r 0 ∥ : p ∈ P k , p ( 0 ) = 1 } . 証明（補題 4.1） Krylov 部分空間の基本性質（命題 2.2）［Krylov 部分空間法 I］ の 1 により x ∈ x 0 + K k \boldsymbol{x}\in\boldsymbol{x}_0+\mathcal{K}_k x ∈ x 0 + K k x = x 0 + q ( A ) r 0 \boldsymbol{x}=\boldsymbol{x}_0+q(A)\boldsymbol{r}_0 x = x 0 + q ( A ) r 0 q ∈ P k − 1 q\in\mathcal{P}_{k-1} q ∈ P k − 1 b − A x = r 0 − A q ( A ) r 0 = p ( A ) r 0 \boldsymbol{b}-A\boldsymbol{x}=\boldsymbol{r}_0-Aq(A)\boldsymbol{r}_0=p(A)\boldsymbol{r}_0 b − A x = r 0 − A q ( A ) r 0 = p ( A ) r 0 p ( t ) = 1 − t q ( t ) p(t)=1-tq(t) p ( t ) = 1 − tq ( t ) p ∈ P k p\in\mathcal{P}_k p ∈ P k p ( 0 ) = 1 p(0)=1 p ( 0 ) = 1 p ( t ) − 1 p(t)-1 p ( t ) − 1 p ( t ) = 1 − t q ( t ) p(t)=1-tq(t) p ( t ) = 1 − tq ( t ) q ∈ P k − 1 q\in\mathcal{P}_{k-1} q ∈ P k − 1 { b − A x } \{\boldsymbol{b}-A\boldsymbol{x}\} { b − A x } { p ( A ) r 0 } \{p(A)\boldsymbol{r}_0\} { p ( A ) r 0 } 定理 3.2 と 命題 3.3 により存在する）は等しいです。
∎ 定理 4.2 （正規行列での GMRES の収束 ） ✓
A A A 1 ≤ k ≤ d 1\le k\le d 1 ≤ k ≤ d
∥ r k ∥ ∥ r 0 ∥ ≤ min p ∈ P k , p ( 0 ) = 1 max λ ∈ σ ( A ) ∣ p ( λ ) ∣ . \frac{\|\boldsymbol{r}_k\|}{\|\boldsymbol{r}_0\|}\le\min_{p\in\mathcal{P}_k,\ p(0)=1}\ \max_{\lambda\in\sigma(A)}|p(\lambda)|. ∥ r 0 ∥ ∥ r k ∥ ≤ p ∈ P k , p ( 0 ) = 1 min λ ∈ σ ( A ) max ∣ p ( λ ) ∣. 証明（定理 4.2） 正規行列のスペクトル定理（定理 8.1）［スペクトル定理］ により A = U D U ∗ A=UDU^{*} A = U D U ∗ U U U D = diag ( λ 1 , … , λ n ) D=\operatorname{diag}(\lambda_1,\ldots,\lambda_n) D = diag ( λ 1 , … , λ n ) A j = U D j U ∗ A^j=UD^jU^{*} A j = U D j U ∗ p ( A ) = U p ( D ) U ∗ p(A)=Up(D)U^{*} p ( A ) = U p ( D ) U ∗ U U U U ∗ U^{*} U ∗ ユニタリ行列の特徴づけ（命題 3.2）［スペクトル定理］ ）から
∥ p ( A ) r 0 ∥ = ∥ p ( D ) U ∗ r 0 ∥ = ( ∑ i ∣ p ( λ i ) ∣ 2 ∣ ( U ∗ r 0 ) i ∣ 2 ) 1 / 2 ≤ max i ∣ p ( λ i ) ∣ ∥ r 0 ∥ \|p(A)\boldsymbol{r}_0\|=\|p(D)U^{*}\boldsymbol{r}_0\|=\Bigl(\sum_i|p(\lambda_i)|^2|(U^{*}\boldsymbol{r}_0)_i|^2\Bigr)^{1/2}\le\max_i|p(\lambda_i)|\,\|\boldsymbol{r}_0\| ∥ p ( A ) r 0 ∥ = ∥ p ( D ) U ∗ r 0 ∥ = ( i ∑ ∣ p ( λ i ) ∣ 2 ∣ ( U ∗ r 0 ) i ∣ 2 ) 1/2 ≤ i max ∣ p ( λ i ) ∣ ∥ r 0 ∥ です。補題 4.1 で p p p
∎ 系 4.3 （対角化可能な行列での GMRES の収束 ） ✓
A = V Λ V − 1 A=V\Lambda V^{-1} A = V Λ V − 1 Λ \Lambda Λ 1 ≤ k ≤ d 1\le k\le d 1 ≤ k ≤ d
∥ r k ∥ ∥ r 0 ∥ ≤ κ ( V ) min p ∈ P k , p ( 0 ) = 1 max λ ∈ σ ( A ) ∣ p ( λ ) ∣ , κ ( V ) = ∥ V ∥ ∥ V − 1 ∥ . \frac{\|\boldsymbol{r}_k\|}{\|\boldsymbol{r}_0\|}\le\kappa(V)\min_{p\in\mathcal{P}_k,\ p(0)=1}\ \max_{\lambda\in\sigma(A)}|p(\lambda)|,\qquad\kappa(V)=\|V\|\|V^{-1}\|. ∥ r 0 ∥ ∥ r k ∥ ≤ κ ( V ) p ∈ P k , p ( 0 ) = 1 min λ ∈ σ ( A ) max ∣ p ( λ ) ∣ , κ ( V ) = ∥ V ∥∥ V − 1 ∥. 証明（系 4.3） p ( A ) = V p ( Λ ) V − 1 p(A)=Vp(\Lambda)V^{-1} p ( A ) = V p ( Λ ) V − 1 ∥ p ( A ) r 0 ∥ ≤ ∥ V ∥ ∥ p ( Λ ) ∥ ∥ V − 1 ∥ ∥ r 0 ∥ \|p(A)\boldsymbol{r}_0\|\le\|V\|\|p(\Lambda)\|\|V^{-1}\|\|\boldsymbol{r}_0\| ∥ p ( A ) r 0 ∥ ≤ ∥ V ∥∥ p ( Λ ) ∥∥ V − 1 ∥∥ r 0 ∥ ∥ p ( Λ ) z ∥ 2 = ∑ i ∣ p ( λ i ) ∣ 2 ∣ z i ∣ 2 \|p(\Lambda)\boldsymbol{z}\|^2=\sum_i|p(\lambda_i)|^2|z_i|^2 ∥ p ( Λ ) z ∥ 2 = ∑ i ∣ p ( λ i ) ∣ 2 ∣ z i ∣ 2 ∥ p ( Λ ) ∥ = max i ∣ p ( λ i ) ∣ \|p(\Lambda)\|=\max_i|p(\lambda_i)| ∥ p ( Λ ) ∥ = max i ∣ p ( λ i ) ∣ 補題 4.1 で最小をとります。
∎ 系 4.4 （固有値が円板に入るときの評価 ） ✓
A = V Λ V − 1 A=V\Lambda V^{-1} A = V Λ V − 1 ∣ z − c ∣ ≤ ρ |z-c|\le\rho ∣ z − c ∣ ≤ ρ 0 ≤ ρ < ∣ c ∣ 0\le\rho<|c| 0 ≤ ρ < ∣ c ∣ 1 ≤ k ≤ d 1\le k\le d 1 ≤ k ≤ d ∥ r k ∥ / ∥ r 0 ∥ ≤ κ ( V ) ( ρ / ∣ c ∣ ) k \|\boldsymbol{r}_k\|/\|\boldsymbol{r}_0\|\le\kappa(V)(\rho/|c|)^k ∥ r k ∥/∥ r 0 ∥ ≤ κ ( V ) ( ρ /∣ c ∣ ) k
証明（系 4.4） p ( z ) = ( 1 − z / c ) k p(z)=(1-z/c)^k p ( z ) = ( 1 − z / c ) k P k \mathcal{P}_k P k p ( 0 ) = 1 p(0)=1 p ( 0 ) = 1 λ \lambda λ ∣ p ( λ ) ∣ = ∣ c − λ ∣ k / ∣ c ∣ k ≤ ( ρ / ∣ c ∣ ) k |p(\lambda)|=|c-\lambda|^k/|c|^k\le(\rho/|c|)^k ∣ p ( λ ) ∣ = ∣ c − λ ∣ k /∣ c ∣ k ≤ ( ρ /∣ c ∣ ) k 系 4.3 に代入します。
∎ 前処理の目標は、この評価が効く形、すなわち「固有値が 0 から離れた小さな円板に集まり、κ ( V ) \kappa(V) κ ( V ) A A A c = 1 c=1 c = 1 M ≈ A M\approx A M ≈ A M − 1 A ≈ I M^{-1}A\approx I M − 1 A ≈ I
次の例のねらいは、固有値の位置が同じでも収束の速さが大きく違いうることを、数値で見ることです。g g g κ ( V ) \kappa(V) κ ( V ) κ ( V ) \kappa(V) κ ( V )
例 4.5 （非正規行列では固有値の評価が役に立たない ）
対角に 0.5 0.5 0.5 1.5 1.5 1.5 g g g A g A_g A g g g g ∣ z − 1 ∣ ≤ 0.5 |z-1|\le0.5 ∣ z − 1∣ ≤ 0.5 系 4.4 の ( ρ / ∣ c ∣ ) k = 0.5 k (\rho/|c|)^k=0.5^k ( ρ /∣ c ∣ ) k = 0. 5 k κ ( V ) \kappa(V) κ ( V ) g = 0 g=0 g = 0 b = ( 1 , … , 1 ) T \boldsymbol{b}=(1,\ldots,1)^{\mathsf{T}} b = ( 1 , … , 1 ) T x 0 = 0 \boldsymbol{x}_0=\boldsymbol{0} x 0 = 0 ∥ r k ∥ / ∥ r 0 ∥ \|\boldsymbol{r}_k\|/\|\boldsymbol{r}_0\| ∥ r k ∥/∥ r 0 ∥ -- は収束して停止した後）。
# Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1
import numpy as np, scipy.sparse.linalg as sla
lam = np. linspace ( 0.5 , 1.5 , n ) # 固有値は円板 |z-1| <= 0.5 の中
for g in gs: # 上の副対角 g が非正規性の強さ
A = np. diag ( lam ) + g * np. eye ( n , k = 1 )
res = [] # 相対残差 |r_k|/|r_0| の履歴
sla. gmres ( A , b , rtol = 1e-13 , restart = n , maxiter = 1 ,
callback = res.append , callback_type = " pr_norm " )
fmt = lambda r , k : f " { r[k - 1 ] :.1e } " if k <= len ( r ) else " -- "
print ( " k 0.5^k g=0 g=1 g=2 " )
for k in ( 5 , 10 , 20 , 40 , 60 , 80 ):
row = " " . join ( fmt ( hist [ g ] , k ) for g in gs )
print ( f " {k :3d } { 0.5 ** k :.1e } {row} " )
5 3.1e-02 1.6e-03 6.0e-03 4.2e-02
10 9.8e-04 2.4e-06 9.3e-04 4.1e-02
20 9.5e-07 4.5e-12 3.8e-05 4.1e-02
40 9.1e-13 -- 6.2e-07 4.1e-02
60 8.7e-19 -- 1.8e-07 4.1e-02
80 8.3e-25 -- 1.7e-07 4.1e-02
固有ベクトル行列は上三角で、成分は閉じた形 v i j = g v i + 1 , j / ( λ j − λ i ) v_{ij}=g\,v_{i+1,j}/(\lambda_j-\lambda_i) v ij = g v i + 1 , j / ( λ j − λ i ) eig ではこの条件数を正しく計算できないので、mpmath の 250 桁で計算します（列は単位ベクトルに正規化）。
# （続き）固有ベクトル行列の条件数を mpmath の 250 桁で計算
import mpmath as mp # mpmath 1.3.0
lm = [ mp. mpf ( 1 ) / 2 + mp. mpf ( i ) / (n - 1 ) for i in range ( n ) ]
for j in range ( n ): # 上三角の固有ベクトル（閉じた形）
for i in range ( j - 1 , - 1 , - 1 ):
V[i, j] = g * V[i + 1 , j] / (lm[j] - lm[i])
V[i, j] /= nj # 列を単位ベクトルに
s = mp. svd_r ( V , compute_uv = False )
print ( f "g= {g} : kappa(V) =" , mp. nstr ( max ( s ) / min ( s ) , 3 ))
g=2: kappa(V) = 3.99e+101
g = 0 g=0 g = 0 0.5 k 0.5^k 0. 5 k g = 1 g=1 g = 1 k = 60 k=60 k = 60 1.7 × 10 − 7 1.7\times10^{-7} 1.7 × 1 0 − 7 g = 2 g=2 g = 2 k = 10 k=10 k = 10 4 × 10 − 2 4\times10^{-2} 4 × 1 0 − 2 g = 1 g=1 g = 1 k = 40 , 60 , 80 , 90 k=40,60,80,90 k = 40 , 60 , 80 , 90 6.2 × 10 − 7 6.2\times10^{-7} 6.2 × 1 0 − 7 1.8 × 10 − 7 1.8\times10^{-7} 1.8 × 1 0 − 7 1.73 × 10 − 7 1.73\times10^{-7} 1.73 × 1 0 − 7 1.73 × 10 − 7 1.73\times10^{-7} 1.73 × 1 0 − 7 k = d ≤ 100 k=d\le100 k = d ≤ 100 命題 3.3 ）ので、停滞はその手前の長い区間で起きていることになります。κ ( V ) ≈ 10 72 , 10 102 \kappa(V)\approx10^{72},10^{102} κ ( V ) ≈ 1 0 72 , 1 0 102
検算. (a) 既知値: g = 0 g=0 g = 0 [ 0.5 , 1.5 ] [0.5,1.5] [ 0.5 , 1.5 ] 定理 4.2 に p ( z ) = T k ( 2 ( 1 − z ) ) / T k ( 2 ) p(z)=T_k(2(1-z))/T_k(2) p ( z ) = T k ( 2 ( 1 − z )) / T k ( 2 ) T k T_k T k ∣ T k ( 2 ( 1 − z ) ) ∣ ≤ 1 |T_k(2(1-z))|\le1 ∣ T k ( 2 ( 1 − z )) ∣ ≤ 1 1 / T k ( 2 ) 1/T_k(2) 1/ T k ( 2 ) T 5 ( 2 ) = 362 T_5(2)=362 T 5 ( 2 ) = 362 T 10 ( 2 ) = 262087 T_{10}(2)=262087 T 10 ( 2 ) = 262087 T 20 ( 2 ) = 2 T 10 ( 2 ) 2 − 1 ≈ 1.37 × 10 11 T_{20}(2)=2T_{10}(2)^2-1\approx1.37\times10^{11} T 20 ( 2 ) = 2 T 10 ( 2 ) 2 − 1 ≈ 1.37 × 1 0 11 2.8 × 10 − 3 2.8\times10^{-3} 2.8 × 1 0 − 3 3.8 × 10 − 6 3.8\times10^{-6} 3.8 × 1 0 − 6 7.3 × 10 − 12 7.3\times10^{-12} 7.3 × 1 0 − 12 1.6 × 10 − 3 1.6\times10^{-3} 1.6 × 1 0 − 3 2.4 × 10 − 6 2.4\times10^{-6} 2.4 × 1 0 − 6 4.5 × 10 − 12 4.5\times10^{-12} 4.5 × 1 0 − 12 κ ( V ) \kappa(V) κ ( V ) g = 0 g=0 g = 0 ∥ r 20 ∥ / ∥ r 10 ∥ ≈ 1.9 × 10 − 6 \|\boldsymbol{r}_{20}\|/\|\boldsymbol{r}_{10}\|\approx1.9\times10^{-6} ∥ r 20 ∥/∥ r 10 ∥ ≈ 1.9 × 1 0 − 6 T 10 ( 2 ) / T 20 ( 2 ) ≈ 1.9 × 10 − 6 T_{10}(2)/T_{20}(2)\approx1.9\times10^{-6} T 10 ( 2 ) / T 20 ( 2 ) ≈ 1.9 × 1 0 − 6
定義 5.1 （左前処理と右前処理 ） ↗
正則な M ∈ C n × n M\in\mathbb{C}^{n\times n} M ∈ C n × n 前処理行列 。M ≈ A M\approx A M ≈ A M z = c M\boldsymbol{z}=\boldsymbol{c} M z = c M − 1 A x = M − 1 b M^{-1}A\boldsymbol{x}=M^{-1}\boldsymbol{b} M − 1 A x = M − 1 b 左前処理 、A M − 1 u = b AM^{-1}\boldsymbol{u}=\boldsymbol{b} A M − 1 u = b x = M − 1 u \boldsymbol{x}=M^{-1}\boldsymbol{u} x = M − 1 u 右前処理 という。
命題 5.2 （右前処理の GMRES は元の残差を最小にする ） ✓
A M − 1 AM^{-1} A M − 1 u 0 = M x 0 \boldsymbol{u}_0=M\boldsymbol{x}_0 u 0 = M x 0 k k k u k \boldsymbol{u}_k u k x k = M − 1 u k \boldsymbol{x}_k=M^{-1}\boldsymbol{u}_k x k = M − 1 u k x k \boldsymbol{x}_k x k
x 0 + K k ( M − 1 A , M − 1 r 0 ) \boldsymbol{x}_0+\mathcal{K}_k(M^{-1}A,\ M^{-1}\boldsymbol{r}_0) x 0 + K k ( M − 1 A , M − 1 r 0 ) の中で ∥ b − A x ∥ \|\boldsymbol{b}-A\boldsymbol{x}\| ∥ b − A x ∥ M − 1 A M^{-1}A M − 1 A M − 1 b M^{-1}\boldsymbol{b} M − 1 b x 0 \boldsymbol{x}_0 x 0 k k k ∥ M − 1 ( b − A x ) ∥ \|M^{-1}(\boldsymbol{b}-A\boldsymbol{x})\| ∥ M − 1 ( b − A x ) ∥ k k k
証明（命題 5.2） 右前処理: 系 A M − 1 u = b AM^{-1}\boldsymbol{u}=\boldsymbol{b} A M − 1 u = b b − A M − 1 u 0 = b − A x 0 = r 0 \boldsymbol{b}-AM^{-1}\boldsymbol{u}_0=\boldsymbol{b}-A\boldsymbol{x}_0=\boldsymbol{r}_0 b − A M − 1 u 0 = b − A x 0 = r 0 u \boldsymbol{u} u b − A M − 1 u = b − A ( M − 1 u ) \boldsymbol{b}-AM^{-1}\boldsymbol{u}=\boldsymbol{b}-A(M^{-1}\boldsymbol{u}) b − A M − 1 u = b − A ( M − 1 u ) u ∈ u 0 + K k ( A M − 1 , r 0 ) \boldsymbol{u}\in\boldsymbol{u}_0+\mathcal{K}_k(AM^{-1},\boldsymbol{r}_0) u ∈ u 0 + K k ( A M − 1 , r 0 ) u ↦ M − 1 u \boldsymbol{u}\mapsto M^{-1}\boldsymbol{u} u ↦ M − 1 u j j j M − 1 ( A M − 1 ) j = ( M − 1 A ) j M − 1 M^{-1}(AM^{-1})^j=(M^{-1}A)^jM^{-1} M − 1 ( A M − 1 ) j = ( M − 1 A ) j M − 1
M − 1 ( u 0 + K k ( A M − 1 , r 0 ) ) = x 0 + span { ( M − 1 A ) j M − 1 r 0 : j < k } = x 0 + K k ( M − 1 A , M − 1 r 0 ) M^{-1}\bigl(\boldsymbol{u}_0+\mathcal{K}_k(AM^{-1},\boldsymbol{r}_0)\bigr)=\boldsymbol{x}_0+\operatorname{span}\{(M^{-1}A)^jM^{-1}\boldsymbol{r}_0:j<k\}=\boldsymbol{x}_0+\mathcal{K}_k(M^{-1}A,M^{-1}\boldsymbol{r}_0) M − 1 ( u 0 + K k ( A M − 1 , r 0 ) ) = x 0 + span {( M − 1 A ) j M − 1 r 0 : j < k } = x 0 + K k ( M − 1 A , M − 1 r 0 ) です。よって x k = M − 1 u k \boldsymbol{x}_k=M^{-1}\boldsymbol{u}_k x k = M − 1 u k ∥ b − A x ∥ \|\boldsymbol{b}-A\boldsymbol{x}\| ∥ b − A x ∥
左前処理: 系 M − 1 A x = M − 1 b M^{-1}A\boldsymbol{x}=M^{-1}\boldsymbol{b} M − 1 A x = M − 1 b M − 1 r 0 M^{-1}\boldsymbol{r}_0 M − 1 r 0 M − 1 ( b − A x ) M^{-1}(\boldsymbol{b}-A\boldsymbol{x}) M − 1 ( b − A x ) 定義 3.1 ）どおりこれを x 0 + K k ( M − 1 A , M − 1 r 0 ) \boldsymbol{x}_0+\mathcal{K}_k(M^{-1}A,M^{-1}\boldsymbol{r}_0) x 0 + K k ( M − 1 A , M − 1 r 0 )
∎ 2 つは同じ部分空間を探しますが、最小にする量が違います。左前処理では、停止判定に使う残差 ∥ M − 1 r k ∥ \|M^{-1}\boldsymbol{r}_k\| ∥ M − 1 r k ∥ M M M ∥ r k ∥ \|\boldsymbol{r}_k\| ∥ r k ∥ gmres に M= を渡すと左前処理になる（前処理つき残差を最小にし、最後に元の残差で収束を確かめる）ので、右前処理にしたいときは A M − 1 AM^{-1} A M − 1 LinearOperator として渡します。
例 5.3 （移流拡散方程式：再始動と ILU 前処理 ）
単位正方形で − Δ u + β ( u x + u y ) = 1 -\Delta u+\beta(u_x+u_y)=1 − Δ u + β ( u x + u y ) = 1 u = 0 u=0 u = 0 h = 1 / ( N + 1 ) h=1/(N+1) h = 1/ ( N + 1 ) β h / 2 \beta h/2 β h /2 10 − 8 10^{-8} 1 0 − 8 spilu（しきい値つき不完全 LU、drop_tol=1e-3、fill_factor=5）を右前処理にした場合で比べます。参照解は spsolve の直接法です。
# Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1
import numpy as np, scipy.sparse as sp, scipy.sparse.linalg as sla
def conv_diff ( N , beta ) : # -Δu + beta (u_x + u_y)、中心差分
h = 1.0 / (N + 1 ); I = sp. eye ( N )
T = sp. diags ( [ - 1.0 , 2.0 , - 1.0 ], [ - 1 , 0 , 1 ], shape = (N, N) ) / h ** 2
C = sp. diags ( [ - 1.0 , 1.0 ], [ - 1 , 1 ], shape = (N, N) ) * beta / ( 2 * h)
return (sp. kron ( I , T + C ) + sp. kron ( T + C , I )). tocsc ()
def run ( A , b , **kw ) : # 内側の反復回数を数える
cb = lambda r : it. __setitem__ ( 0 , it [ 0 ] + 1 )
x, info = sla. gmres ( A , b , rtol = 1e-8 , callback = cb ,
callback_type = " pr_norm " , ** kw )
print ( " N beta full GMRES(20) ILU-right fill err(right) " )
for N, beta in (( 32 , 10.0 ), ( 32 , 100.0 ), ( 64 , 10.0 ), ( 64 , 100.0 )):
A = conv_diff ( N , beta ); n = N * N
b = np. ones ( n ); xs = sla. spsolve ( A , b ) # 参照解
k1 = run ( A , b , restart = n , maxiter = 1 ) [ 0 ]
k2 = run ( A , b , restart = 20 , maxiter = 500 ) [ 0 ]
ilu = sla. spilu ( A , drop_tol = 1e-3 , fill_factor = 5 )
AM = sla. LinearOperator ( (n, n) , matvec =lambda v : A @ ilu. solve ( v ))
k3, u = run ( AM , b , restart = n , maxiter = 1 ) # 右前処理
x = ilu. solve ( u ) # x = M^{-1} u
e = np.linalg. norm ( x - xs ) / np.linalg. norm ( xs )
fill = (ilu.L.nnz + ilu.U.nnz) / A.nnz
print ( f " {N :3d } {beta :5.0f } {k1 :4d } {k2 :9d } {k3 :9d } "
N beta full GMRES(20) ILU-right fill err(right)
32 10 91 206 35 3.5 4.3e-10
32 100 71 217 20 3.3 5.9e-10
64 10 181 299 76 3.5 6.2e-10
64 100 131 224 40 3.2 1.4e-10
再始動なしの GMRES の反復回数は N N N h − 1 h^{-1} h − 1
検算. (a) 既知値: 右辺と境界条件を固定した離散問題の厳密解（直接法の解）を参照値とし、右前処理の解との相対差は 10 − 10 10^{-10} 1 0 − 10 10 − 8 10^{-8} 1 0 − 8 spsolve）で計算しました。(c) 次数: full GMRES の反復回数の比は β = 10 \beta=10 β = 10 181 / 91 ≈ 1.99 181/91\approx1.99 181/91 ≈ 1.99 β = 100 \beta=100 β = 100 131 / 71 ≈ 1.85 131/71\approx1.85 131/71 ≈ 1.85 h − 1 h^{-1} h − 1
対称正定値の系には、Cholesky 分解 A = L L T A=LL^{\mathsf{T}} A = L L T Cholesky 分解の存在（定理 2.1）［Cholesky 分解と帯行列・疎行列の直接法］ ）の fill-in を捨てた近似 M = L L T M=LL^{\mathsf{T}} M = L L T
定義 6.1 （Z 行列と Stieltjes 行列 ） ↗
非対角成分がすべて 0 0 0 Z 行列 という。正定値な対称 Z 行列を Stieltjes 行列 （対称な M 行列）という。実行列 N , N ′ N,N' N , N ′ N ≥ 0 N\ge0 N ≥ 0 N ′ ≤ N N'\le N N ′ ≤ N
離散 Laplace 行列（2 次元・3 次元を含む）は Stieltjes 行列の代表例です。
次の補題は、成分を小さくした非負行列のスペクトル半径は大きくならない、という事実です。IC の算法で成分を捨てたあとも正定値のままか、を確かめる場面で使います。証明の方針は、N / t N/t N / t N ′ / t N'/t N ′ / t
補題 6.2 （非負行列のスペクトル半径の単調性 ） ✓
N , N ′ ∈ R n × n N,N'\in\mathbb{R}^{n\times n} N , N ′ ∈ R n × n 0 ≤ N ′ ≤ N 0\le N'\le N 0 ≤ N ′ ≤ N ρ ( N ′ ) ≤ ρ ( N ) \rho(N')\le\rho(N) ρ ( N ′ ) ≤ ρ ( N )
証明（補題 6.2） t > ρ ( N ) t>\rho(N) t > ρ ( N ) k k k 0 ≤ ( N ′ / t ) k ≤ ( N / t ) k 0\le(N'/t)^k\le(N/t)^k 0 ≤ ( N ′ / t ) k ≤ ( N / t ) k ρ ( N / t ) < 1 \rho(N/t)<1 ρ ( N / t ) < 1 行列の冪の収束（系 5.3）［行列のノルムと線形方程式の条件数］ により ( N / t ) k → 0 (N/t)^k\to0 ( N / t ) k → 0 ( N ′ / t ) k → 0 (N'/t)^k\to0 ( N ′ / t ) k → 0 ρ ( N ′ / t ) < 1 \rho(N'/t)<1 ρ ( N ′ / t ) < 1 ρ ( N ′ ) < t \rho(N')<t ρ ( N ′ ) < t t > ρ ( N ) t>\rho(N) t > ρ ( N ) ρ ( N ′ ) ≤ ρ ( N ) \rho(N')\le\rho(N) ρ ( N ′ ) ≤ ρ ( N )
∎ 補題 6.3 （Stieltjes 行列の特徴づけ ） ✓
A A A s ≥ max i a i i s\ge\max_ia_{ii} s ≥ max i a ii N = s I − A N=sI-A N = s I − A N ≥ 0 N\ge0 N ≥ 0 A A A ⟺ ρ ( N ) < s \iff\rho(N)<s ⟺ ρ ( N ) < s
証明（補題 6.3） N N N s − a i i ≥ 0 s-a_{ii}\ge0 s − a ii ≥ 0 − a i j ≥ 0 -a_{ij}\ge0 − a ij ≥ 0 N N N 主軸定理（定理 6.3）［スペクトル定理］ により実の正規直交固有ベクトルの基底があり、固有値 μ 1 ≥ ⋯ ≥ μ n \mu_1\ge\cdots\ge\mu_n μ 1 ≥ ⋯ ≥ μ n A A A s − μ i s-\mu_i s − μ i A A A ⟺ \iff ⟺ i i i s − μ i > 0 ⟺ μ 1 < s s-\mu_i>0\iff\mu_1<s s − μ i > 0 ⟺ μ 1 < s μ 1 = ρ ( N ) \mu_1=\rho(N) μ 1 = ρ ( N ) N x = μ x N\boldsymbol{x}=\mu\boldsymbol{x} N x = μ x x \boldsymbol{x} x N ≥ 0 N\ge0 N ≥ 0 Rayleigh 商の最大・最小（定理 7.1）［スペクトル定理］ から
∣ μ ∣ = ∣ x T N x ∣ ≤ ∣ x ∣ T N ∣ x ∣ ≤ μ 1 |\mu|=|\boldsymbol{x}^{\mathsf{T}}N\boldsymbol{x}|\le|\boldsymbol{x}|^{\mathsf{T}}N|\boldsymbol{x}|\le\mu_1 ∣ μ ∣ = ∣ x T N x ∣ ≤ ∣ x ∣ T N ∣ x ∣ ≤ μ 1 です（∣ x ∣ |\boldsymbol{x}| ∣ x ∣ ρ ( N ) ≤ μ 1 ≤ ρ ( N ) \rho(N)\le\mu_1\le\rho(N) ρ ( N ) ≤ μ 1 ≤ ρ ( N )
∎ IC の算法（ねらい）. Cholesky 分解を 1 列ずつ進めると、各段で「まだ分解していない部分」を表す小さな行列（Schur 補行列）が残ります。IC は、その行列の成分のうち、あらかじめ決めた非零パターンの外にあるものを 0 に捨てて次の段へ進みます。捨てると正定値性が壊れて平方根の中が負になりかねないので、次の定理は、Stieltjes 行列ならそれが起きないことを段ごとの帰納法で示します。
IC の算法. 非零パターン S ⊂ { 1 , … , n } 2 S\subset\{1,\ldots,n\}^2 S ⊂ { 1 , … , n } 2 ( i , j ) ∈ S ⇒ ( j , i ) ∈ S (i,j)\in S\Rightarrow(j,i)\in S ( i , j ) ∈ S ⇒ ( j , i ) ∈ S A A A S S S A A A B ( 1 ) = A B^{(1)}=A B ( 1 ) = A k = 1 , … , n k=1,\ldots,n k = 1 , … , n n − k + 1 n-k+1 n − k + 1 B ( k ) B^{(k)} B ( k ) k , … , n k,\ldots,n k , … , n
B ( k ) = ( d k w k T w k C k ) , l k k = d k , l i k = ( w k ) i d k ( i > k ) B^{(k)}=\begin{pmatrix}d_k&\boldsymbol{w}_k^{\mathsf{T}}\\\boldsymbol{w}_k&C_k\end{pmatrix},\qquad l_{kk}=\sqrt{d_k},\quad l_{ik}=\frac{(\boldsymbol{w}_k)_i}{\sqrt{d_k}}\ (i>k) B ( k ) = ( d k w k w k T C k ) , l k k = d k , l ik = d k ( w k ) i ( i > k ) と分け、Schur 補行列 S k = C k − w k w k T / d k S_k=C_k-\boldsymbol{w}_k\boldsymbol{w}_k^{\mathsf{T}}/d_k S k = C k − w k w k T / d k ( i , j ) ∉ S (i,j)\notin S ( i , j ) ∈ / S 0 0 0 B ( k + 1 ) B^{(k+1)} B ( k + 1 ) F k F_k F k B ( k + 1 ) = S k + F k B^{(k+1)}=S_k+F_k B ( k + 1 ) = S k + F k
定理 6.4 （Stieltjes 行列の不完全 Cholesky 分解は破綻しない ） ✓
A A A d k > 0 d_k>0 d k > 0 B ( k ) B^{(k)} B ( k ) L = ( l i j ) L=(l_{ij}) L = ( l ij ) ( i , j ) ∉ S (i,j)\notin S ( i , j ) ∈ / S l i j = 0 l_{ij}=0 l ij = 0
L L T = A + R , R = ∑ k F k ( 位置を埋め込んだ和 ) ≥ 0 , ( i , j ) ∈ S で r i j = 0 LL^{\mathsf{T}}=A+R,\qquad R=\sum_kF_k\ (\text{位置を埋め込んだ和})\ \ge0,\qquad(i,j)\in S\text{ で }r_{ij}=0 L L T = A + R , R = k ∑ F k ( 位置を埋め込んだ和 ) ≥ 0 , ( i , j ) ∈ S で r ij = 0 が成り立つ。とくに M = L L T M=LL^{\mathsf{T}} M = L L T
証明（定理 6.4） 帰納法の主張. B ( k ) B^{(k)} B ( k ) ( i , j ) ∉ S (i,j)\notin S ( i , j ) ∈ / S k = 1 k=1 k = 1 S S S A A A
d k > 0 d_k>0 d k > 0 e 1 T B ( k ) e 1 > 0 \boldsymbol{e}_1^{\mathsf{T}}B^{(k)}\boldsymbol{e}_1>0 e 1 T B ( k ) e 1 > 0
Schur 補行列. 行と列の順序を入れ替えた ( C k w k w k T d k ) \begin{pmatrix}C_k&\boldsymbol{w}_k\\\boldsymbol{w}_k^{\mathsf{T}}&d_k\end{pmatrix} ( C k w k T w k d k ) B ( k ) B^{(k)} B ( k ) x T P T B P x = ( P x ) T B ( P x ) \boldsymbol{x}^{\mathsf{T}}P^{\mathsf{T}}BP\boldsymbol{x}=(P\boldsymbol{x})^{\mathsf{T}}B(P\boldsymbol{x}) x T P T B P x = ( P x ) T B ( P x ) Schur 補行列（定理 4.2）［行列の固有値評価と不等式］ の 3（右下の d k > 0 d_k>0 d k > 0 S k = C k − w k d k − 1 w k T S_k=C_k-\boldsymbol{w}_kd_k^{-1}\boldsymbol{w}_k^{\mathsf{T}} S k = C k − w k d k − 1 w k T ( w k ) i , ( w k ) j ≤ 0 (\boldsymbol{w}_k)_i,(\boldsymbol{w}_k)_j\le0 ( w k ) i , ( w k ) j ≤ 0 ( w k ) i ( w k ) j ≥ 0 (\boldsymbol{w}_k)_i(\boldsymbol{w}_k)_j\ge0 ( w k ) i ( w k ) j ≥ 0 ( S k ) i j ≤ ( C k ) i j ≤ 0 (S_k)_{ij}\le(C_k)_{ij}\le0 ( S k ) ij ≤ ( C k ) ij ≤ 0 S k S_k S k
捨てる操作. s = max i ( S k ) i i s=\max_i(S_k)_{ii} s = max i ( S k ) ii S k = s I − N S_k=sI-N S k = s I − N 補題 6.3 により N ≥ 0 N\ge0 N ≥ 0 ρ ( N ) < s \rho(N)<s ρ ( N ) < s B ( k + 1 ) B^{(k+1)} B ( k + 1 ) S k S_k S k B ( k + 1 ) = s I − N ′ B^{(k+1)}=sI-N' B ( k + 1 ) = s I − N ′ 0 ≤ N ′ ≤ N 0\le N'\le N 0 ≤ N ′ ≤ N 補題 6.2 で ρ ( N ′ ) ≤ ρ ( N ) < s \rho(N')\le\rho(N)<s ρ ( N ′ ) ≤ ρ ( N ) < s 補題 6.3 で B ( k + 1 ) B^{(k+1)} B ( k + 1 ) ( i , j ) ∉ S (i,j)\notin S ( i , j ) ∈ / S
分解の式. ( w k ) i (\boldsymbol{w}_k)_i ( w k ) i B ( k ) B^{(k)} B ( k ) ( i , k ) (i,k) ( i , k ) ( i , k ) ∉ S (i,k)\notin S ( i , k ) ∈ / S 0 0 0 l i k = 0 l_{ik}=0 l ik = 0 l k = ( l k k , … , l n k ) T \boldsymbol{l}_k=(l_{kk},\ldots,l_{nk})^{\mathsf{T}} l k = ( l k k , … , l nk ) T B ( k ) − l k l k T = ( 0 0 T 0 S k ) = ( 0 0 T 0 B ( k + 1 ) − F k ) B^{(k)}-\boldsymbol{l}_k\boldsymbol{l}_k^{\mathsf{T}}=\begin{pmatrix}0&\boldsymbol{0}^{\mathsf{T}}\\\boldsymbol{0}&S_k\end{pmatrix}=\begin{pmatrix}0&\boldsymbol{0}^{\mathsf{T}}\\\boldsymbol{0}&B^{(k+1)}-F_k\end{pmatrix} B ( k ) − l k l k T = ( 0 0 0 T S k ) = ( 0 0 0 T B ( k + 1 ) − F k ) n × n n\times n n × n k = 1 , … , n k=1,\ldots,n k = 1 , … , n A = ∑ k l k l k T − ∑ k F k = L L T − R A=\sum_k\boldsymbol{l}_k\boldsymbol{l}_k^{\mathsf{T}}-\sum_kF_k=LL^{\mathsf{T}}-R A = ∑ k l k l k T − ∑ k F k = L L T − R F k F_k F k ( i , j ) ∉ S (i,j)\notin S ( i , j ) ∈ / S − ( S k ) i j ≥ 0 -(S_k)_{ij}\ge0 − ( S k ) ij ≥ 0 L L L x T L L T x = ∥ L T x ∥ 2 > 0 \boldsymbol{x}^{\mathsf{T}}LL^{\mathsf{T}}\boldsymbol{x}=\|L^{\mathsf{T}}\boldsymbol{x}\|^2>0 x T L L T x = ∥ L T x ∥ 2 > 0 x ≠ 0 \boldsymbol{x}\ne\boldsymbol{0} x = 0
∎ S S S
例 6.5 （2 次元 Laplace 方程式での IC(0) 前処理 ）
2 次元の離散 Laplace 行列（tridiag ( − 1 , 2 , − 1 ) \operatorname{tridiag}(-1,2,-1) tridiag ( − 1 , 2 , − 1 ) 10 − 8 10^{-8} 1 0 − 8 spsolve です。
行ごとの版が上の算法と同じ L L L ( i , j ) ∈ S (i,j)\in S ( i , j ) ∈ S i ≥ j i\ge j i ≥ j k < j k<j k < j ( w k ) i ( w k ) j / d k = l i k l j k (\boldsymbol{w}_k)_i(\boldsymbol{w}_k)_j/d_k=l_{ik}l_{jk} ( w k ) i ( w k ) j / d k = l ik l j k S S S j j j B ( j ) B^{(j)} B ( j ) a i j − ∑ k < j l i k l j k a_{ij}-\sum_{k<j}l_{ik}l_{jk} a ij − ∑ k < j l ik l j k d j = l j j \sqrt{d_j}=l_{jj} d j = l j j l i j l_{ij} l ij
l j j = ( a j j − ∑ k < j l j k 2 ) 1 / 2 , l i j = 1 l j j ( a i j − ∑ k < j l i k l j k ) ( ( i , j ) ∈ S , i > j ) l_{jj}=\Bigl(a_{jj}-\sum_{k<j}l_{jk}^2\Bigr)^{1/2},\qquad l_{ij}=\frac1{l_{jj}}\Bigl(a_{ij}-\sum_{k<j}l_{ik}l_{jk}\Bigr)\quad((i,j)\in S,\ i>j) l j j = ( a j j − k < j ∑ l j k 2 ) 1/2 , l ij = l j j 1 ( a ij − k < j ∑ l ik l j k ) (( i , j ) ∈ S , i > j ) で、( i , j ) ∉ S (i,j)\notin S ( i , j ) ∈ / S l i j = 0 l_{ij}=0 l ij = 0 定理 6.4 の証明の「分解の式」）。下のコードはこの式を i i i j j j l l l
# Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1
import numpy as np, scipy.sparse as sp, scipy.sparse.linalg as sla
def ic0 ( A ) : # 下三角部分の非零パターンだけで分解
T = sp. tril ( A , format = " csr " ); n = A.shape[ 0 ]
L = [ dict () for _ in range ( n ) ] # L[i] = {列: 値}
lo, hi = T.indptr[i], T.indptr[i + 1 ]
for j, a in sorted ( zip ( T.indices [ lo:hi ] , T.data [ lo:hi ])):
s = a - sum ( v * L [ j ] . get ( k , 0.0 )
for k , v in L [ i ] . items () if k < j )
L[i] [ j ] = np. sqrt ( s ) if j == i else s / L[j] [ j ]
t = [ (i, k, x) for i in range ( n ) for k, x in L[i]. items () ]
return sp. csr_matrix ( (v, (r, c)) , shape = (n, n) )
T1 = lambda N : sp. diags ( [ - 1.0 , 2.0 , - 1.0 ], [ - 1 , 0 , 1 ], shape = (N, N) )
print ( " N CG CG+IC(0) ratio err(IC) " )
for N in ( 16 , 32 , 64 , 128 ):
I = sp. eye ( N ); A = (sp. kron ( I , T1 ( N )) + sp. kron ( T1 ( N ) , I )). tocsr ()
b = np. ones ( N * N ); xs = sla. spsolve ( A. tocsc () , b ) # 参照解
L = ic0 ( A ); Lt = L.T. tocsr ()
sol = lambda r : sla. spsolve_triangular (
Lt , sla. spsolve_triangular ( L , r , lower = True ) , lower = False )
M = sla. LinearOperator ( A.shape , matvec = sol ); cnt = []
for kw in ({}, { " M " : M}):
it = [ 0 ] ; cb = lambda x : it. __setitem__ ( 0 , it [ 0 ] + 1 )
x, _ = sla. cg ( A , b , rtol = 1e-8 , callback = cb , ** kw )
e = np.linalg. norm ( x - xs ) / np.linalg. norm ( xs )
print ( f " {N :4d } {cnt [ 0 ] :4d } {cnt [ 1 ] :7d } {cnt [ 0 ] / cnt [ 1 ] :5.2f } "
N CG CG+IC(0) ratio err(IC)
前処理なしの反復回数は N N N κ ∝ h − 1 \sqrt{\kappa}\propto h^{-1} κ ∝ h − 1 h − 1 h^{-1} h − 1 1.7 → 1.8 → 1.9 1.7\to1.8\to1.9 1.7 → 1.8 → 1.9 h − 2 h^{-2} h − 2
検算. (a) 既知値: 定理 6.4 のとおり、どの格子でも平方根の中は正で分解は最後まで進みました。前処理つきの解は spsolve の参照解と相対 10 − 10 10^{-10} 1 0 − 10 59 / 28 , 119 / 59 , 239 / 119 59/28,119/59,239/119 59/28 , 119/59 , 239/119 2.1 , 2.0 , 2.0 2.1,2.0,2.0 2.1 , 2.0 , 2.0 ∝ κ \propto\sqrt{\kappa} ∝ κ ノルムの章 の κ ≈ 4 / ( π 2 h 2 ) \kappa\approx4/(\pi^2h^2) κ ≈ 4/ ( π 2 h 2 ) h − 1 h^{-1} h − 1
反復法が行列に要求するのは行列ベクトル積だけです。格納形式と演算量は Cholesky 分解と疎行列の直接法 の 疎行列の格納形式（定義 6.1）［Cholesky 分解と帯行列・疎行列の直接法］ と 疎行列ベクトル積の演算量（命題 6.2）［Cholesky 分解と帯行列・疎行列の直接法］ で扱いました。反復法の観点でまとめると次のとおりです（n n z \mathrm{nnz} nnz k k k
形式・道具 記憶 行列ベクトル積 反復法での使いどころ COO（coo_array） 3 n n z 3\,\mathrm{nnz} 3 nnz 2 n n z 2\,\mathrm{nnz} 2 nnz 組み立て。計算の前に CSR・CSC に変換する CSR（csr_array） 2 n n z + n + 1 2\,\mathrm{nnz}+n+1 2 nnz + n + 1 2 n n z 2\,\mathrm{nnz} 2 nnz A x A\boldsymbol{x} A x CSC（csc_array） 2 n n z + n + 1 2\,\mathrm{nnz}+n+1 2 nnz + n + 1 2 n n z 2\,\mathrm{nnz} 2 nnz splu・spilu・spsolve の入力LinearOperator関数だけ 関数の手間 行列を作らない演算子、A M − 1 AM^{-1} A M − 1 M − 1 M^{-1} M − 1 GMRES（再始動なし） O ( n k ) O(nk) O ( nk ) k k k 直交化 O ( n k 2 ) O(nk^2) O ( n k 2 ) k k k 共役勾配法 O ( n ) O(n) O ( n ) k k k 正定値。短い漸化式で直交化は不要
LinearOperator((n, n), matvec=f) は、ベクトルを受け取って A v A\boldsymbol{v} A v f だけで行列の代わりになります。上の 例 5.3 では A M − 1 AM^{-1} A M − 1 例 6.5 では ( L L T ) − 1 (LL^{\mathsf{T}})^{-1} ( L L T ) − 1 J v J\boldsymbol{v} J v 2 n n z 2\,\mathrm{nnz} 2 nnz
演習 8.1 易
A A A λ 1 , λ 2 \lambda_1,\lambda_2 λ 1 , λ 2
解答 p ( z ) = ( 1 − z / λ 1 ) ( 1 − z / λ 2 ) p(z)=(1-z/\lambda_1)(1-z/\lambda_2) p ( z ) = ( 1 − z / λ 1 ) ( 1 − z / λ 2 ) P 2 \mathcal{P}_2 P 2 p ( 0 ) = 1 p(0)=1 p ( 0 ) = 1 p ( λ 1 ) = p ( λ 2 ) = 0 p(\lambda_1)=p(\lambda_2)=0 p ( λ 1 ) = p ( λ 2 ) = 0 系 4.3 から ∥ r 2 ∥ ≤ κ ( V ) ⋅ 0 = 0 \|\boldsymbol{r}_2\|\le\kappa(V)\cdot0=0 ∥ r 2 ∥ ≤ κ ( V ) ⋅ 0 = 0 d ≥ 2 d\ge2 d ≥ 2 d = 1 d=1 d = 1 命題 3.3 により 1 回で厳密解です。
演習 8.2 標準
A A A A e i = e i + 1 A\boldsymbol{e}_i=\boldsymbol{e}_{i+1} A e i = e i + 1 A e n = e 1 A\boldsymbol{e}_n=\boldsymbol{e}_1 A e n = e 1 b = e 1 \boldsymbol{b}=\boldsymbol{e}_1 b = e 1 x 0 = 0 \boldsymbol{x}_0=\boldsymbol{0} x 0 = 0 k < n k<n k < n ∥ r k ∥ = 1 \|\boldsymbol{r}_k\|=1 ∥ r k ∥ = 1 k = n k=n k = n r n = 0 \boldsymbol{r}_n=\boldsymbol{0} r n = 0 A A A 定理 4.2 と矛盾しないことを説明せよ。
解答 K k ( A , e 1 ) = span { e 1 , … , e k } \mathcal{K}_k(A,\boldsymbol{e}_1)=\operatorname{span}\{\boldsymbol{e}_1,\ldots,\boldsymbol{e}_k\} K k ( A , e 1 ) = span { e 1 , … , e k } A K k = span { e 2 , … , e k + 1 } A\mathcal{K}_k=\operatorname{span}\{\boldsymbol{e}_2,\ldots,\boldsymbol{e}_{k+1}\} A K k = span { e 2 , … , e k + 1 } k < n k<n k < n e 1 \boldsymbol{e}_1 e 1 x ∈ K k \boldsymbol{x}\in\mathcal{K}_k x ∈ K k ∥ e 1 − A x ∥ 2 = 1 + ∥ A x ∥ 2 ≥ 1 \|\boldsymbol{e}_1-A\boldsymbol{x}\|^2=1+\|A\boldsymbol{x}\|^2\ge1 ∥ e 1 − A x ∥ 2 = 1 + ∥ A x ∥ 2 ≥ 1 x = 0 \boldsymbol{x}=\boldsymbol{0} x = 0 ∥ r k ∥ = 1 \|\boldsymbol{r}_k\|=1 ∥ r k ∥ = 1 k = n k=n k = n x = e n \boldsymbol{x}=\boldsymbol{e}_n x = e n A x = e 1 A\boldsymbol{x}=\boldsymbol{e}_1 A x = e 1 A A A n n n ω j \omega^j ω j p ( 0 ) = 1 p(0)=1 p ( 0 ) = 1 deg p ≤ k < n \deg p\le k<n deg p ≤ k < n 1 n ∑ j p ( ω j ) = p ( 0 ) = 1 \frac1n\sum_j p(\omega^j)=p(0)=1 n 1 ∑ j p ( ω j ) = p ( 0 ) = 1 ω j \omega^j ω j max j ∣ p ( ω j ) ∣ ≥ 1 \max_j|p(\omega^j)|\ge1 max j ∣ p ( ω j ) ∣ ≥ 1
演習 8.3 標準
(1) A = ( 1 − 2 − 2 1 ) A=\begin{pmatrix}1&-2\\-2&1\end{pmatrix} A = ( 1 − 2 − 2 1 ) S S S d 2 < 0 d_2<0 d 2 < 0 定理 6.4 のどの仮定が欠けているか答えよ。(2) 同じ A A A 補題 6.3 の N N N ρ ( N ) \rho(N) ρ ( N )
解答 (1) d 1 = 1 d_1=1 d 1 = 1 w 1 = ( − 2 ) \boldsymbol{w}_1=(-2) w 1 = ( − 2 ) d 2 = 1 − ( − 2 ) 2 / 1 = − 3 < 0 d_2=1-(-2)^2/1=-3<0 d 2 = 1 − ( − 2 ) 2 /1 = − 3 < 0 A A A 1 ± 2 1\pm2 1 ± 2 3 , − 1 3,-1 3 , − 1 s = 1 s=1 s = 1 N = I − A = ( 0 2 2 0 ) N=I-A=\begin{pmatrix}0&2\\2&0\end{pmatrix} N = I − A = ( 0 2 2 0 ) ± 2 \pm2 ± 2 ρ ( N ) = 2 > s = 1 \rho(N)=2>s=1 ρ ( N ) = 2 > s = 1
Y. Saad, Iterative Methods for Sparse Linear Systems , 2nd ed., SIAM, 2003 — Ch. 6（Arnoldi・GMRES）、Ch. 9–10（前処理、不完全分解）。
A. Greenbaum, Iterative Methods for Solving Linear Systems , SIAM, 1997 — Ch. 2–3（GMRES・MINRES と非正規行列での収束）。
Y. Saad and M. H. Schultz, “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems”, SIAM J. Sci. Stat. Comput. 7 (1986), 856–869.
J. A. Meijerink and H. A. van der Vorst, “An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix”, Math. Comp. 31 (1977), 148–162.
A. Greenbaum, V. Pták, and Z. Strakoš, “Any nonincreasing convergence curve is possible for GMRES”, SIAM J. Matrix Anal. Appl. 17 (1996), 465–469.
U. Trottenberg, C. W. Oosterlee, and A. Schüller, Multigrid , Academic Press, 2001.
執筆: Claude Code（Opus 5.5）、Claude Code（Sonnet 5.5）・監査 2 回（Claude Code）・最終 2026-10-12