# Krylov 部分空間法 II：GMRES・前処理・疎行列の反復解法 > GMRES の最小残差性と Givens 回転による実装、厳密な停止、正規行列・対角化可能な行列での収束評価を証明し、右前処理が真の残差を最小にすること、Stieltjes 行列の不完全 Cholesky 分解 IC(0) が破綻しないことを示す。移流拡散方程式と非正規行列で数値的に確かめる。 > https://rikai.mugen-giken.com/mathematics/numerical-linear-algebra/gmres-and-preconditioning ## 0. この記事の要点 連立一次方程式 $A\boldsymbol{x}=\boldsymbol{b}$ の未知数が何百万もあるとき、行列を丸ごと分解するのは重すぎます。そこで、行列とベクトルの積だけを繰り返して、答えに少しずつ近づく方法を使います。この記事では、そのときまでに作れる候補の中で誤差（残差）が一番小さいものを選ぶ GMRES という方法と、近づく速さを上げる前処理を扱います。読み終えると、速く収束する場合と遅い場合の違い、前処理を置いても計算が途中で壊れない条件を、証明つきで説明できるようになります。 - **GMRES** は $\boldsymbol{x}_0+\mathcal{K}_k(A,\boldsymbol{r}_0)$ の上で残差 $\|\boldsymbol{b}-A\boldsymbol{x}\|_2$ を最小にします。Arnoldi の関係式により、これは $(k+1)\times k$ の Hessenberg 行列（対角の 1 つ下までしか成分をもたない行列）の最小二乗問題になり、[Givens 回転](/mathematics/numerical-linear-algebra/qr-factorization#def-givens)で 1 反復あたり $O(k)$ の追加計算で逐次に解けます。 - 厳密な算術では、Arnoldi の過程が止まった時点で厳密解に達し、$n$ 回以内に終わります。 - 残差は $\boldsymbol{r}_k=p_k(A)\boldsymbol{r}_0$（$p_k(0)=1$）の形で、[正規行列](/mathematics/linear-algebra/spectral-theorem#def-hermitian-unitary)では $\|\boldsymbol{r}_k\|/\|\boldsymbol{r}_0\|\le\min_p\max_{\lambda\in\sigma(A)}|p(\lambda)|$、[対角化可能](/mathematics/linear-algebra/diagonalization-and-jordan-form#def-diagonalizable)な行列ではこれに $\kappa(V)$ が掛かります。固有値が円板 $|z-c|\le\rho<|c|$ に入れば $(\rho/|c|)^k\kappa(V)$ です。非正規性が強いと $\kappa(V)$ が天文学的になり、固有値だけでは収束を予言できません。 - **右前処理** $AM^{-1}\boldsymbol{u}=\boldsymbol{b}$ の GMRES は、左前処理と同じ部分空間の上で**元の系の残差**を最小にします。 - 対称な M 行列（Stieltjes 行列。非対角成分がすべて 0 以下の正定値な対称行列）では、非零パターンを固定した**不完全 Cholesky 分解** IC(0)（Cholesky 分解の途中で新しくできる非零成分を捨てる近似）が破綻しません。証明の鍵は Schur 補行列（1 段ぶん消去したあとに残る小さな行列）が Stieltjes 行列であり続けることと、非負行列の[スペクトル半径](/mathematics/numerical-linear-algebra/norms-and-conditioning#def-spectral-radius)の単調性です。 ## 1. 動機 偏微分方程式を離散化すると、未知数が $10^5$〜$10^8$ の疎な線形方程式が現れます。直接法（[Cholesky 分解と疎行列の直接法](/mathematics/numerical-linear-algebra/cholesky-and-sparse-direct)）は fill-in（分解の途中で、もとは零だった成分が非零になること）のため 3 次元の問題では記憶容量が足りなくなり、反復法が必要になります。対称正定値なら共役勾配法（[Krylov 部分空間法 I](/mathematics/numerical-linear-algebra/lanczos-and-conjugate-gradient)）が標準ですが、移流のある流れの問題・非対称な離散化・固有値問題のシフト反転では行列が非対称になり、共役勾配法の最適性は使えません。 GMRES は「Krylov 部分空間の上で残差を最小にする」という最も素直な一般化です。その代わり、Lanczos の三項漸化式が使えないので、反復ごとに全部の基底ベクトルと直交化する必要があり、記憶と演算が反復回数とともに増えます。そこで実務では、(1) 反復回数そのものを減らす**前処理**と、(2) 一定回数で打ち切る**再始動**が欠かせません。この記事では GMRES の最適性と収束の理論、前処理の 2 つの置き方の違い、代表的な前処理である不完全 Cholesky 分解がいつ安全に計算できるかを証明します。 ## 2. 準備 $A\in\mathbb{C}^{n\times n}$ は正則、$\boldsymbol{b}\in\mathbb{C}^n$、初期値 $\boldsymbol{x}_0$、$\boldsymbol{r}_0=\boldsymbol{b}-A\boldsymbol{x}_0\ne\boldsymbol{0}$、$\beta=\|\boldsymbol{r}_0\|$ とします。ノルムは 2 ノルム、$\mathcal{P}_k$ は次数 $k$ 以下の多項式の全体です。Krylov 部分空間 $\mathcal{K}_k(A,\boldsymbol{r}_0)$、$\boldsymbol{r}_0$ の次数 $d$、Arnoldi の過程の $Q_k$、$H_k$、$\tilde H_k$ は [Krylov 部分空間法 I](/mathematics/numerical-linear-algebra/lanczos-and-conjugate-gradient) の と の記号（$\boldsymbol{b}$ を $\boldsymbol{r}_0$ に置き換える）をそのまま使います。中心になるのは次の事実です。 - : $\mathcal{K}_k=\{q(A)\boldsymbol{r}_0:q\in\mathcal{P}_{k-1}\}$、$d\le n$。 - : $k0$、$Q_k$ の列は $\mathcal{K}_k$ の正規直交基底、$AQ_k=Q_{k+1}\tilde H_k$。$k=d$ では $AQ_d=Q_dH_d$。 ## 3. GMRES $1\le k\le d$ について、$\boldsymbol{x}_0+\mathcal{K}_k(A,\boldsymbol{r}_0)$ の中で $\|\boldsymbol{b}-A\boldsymbol{x}\|$ を最小にする $\boldsymbol{x}_k$ を GMRES の第 $k$ 近似解、$\boldsymbol{r}_k=\boldsymbol{b}-A\boldsymbol{x}_k$ をその残差という（存在と一意性は下の定理で示す）。 次の定理のねらいは 2 つです。GMRES の近似解が、Arnoldi の関係式を通して小さな最小二乗問題に置き換わること、そしてその問題が Givens 回転で 1 反復ずつ安く更新できることです。Givens 回転は 2 つの行だけを混ぜて 1 つの成分を 0 にする操作で、反復が進んで新しい列が 1 本増えたときも、その列に回転をかけ直して新しい回転を 1 つ足せば済みます。 $1\le k (1) の 1 により、$\boldsymbol{x}\in\boldsymbol{x}_0+\mathcal{K}_k$ は $\boldsymbol{x}=\boldsymbol{x}_0+Q_k\boldsymbol{y}$（$\boldsymbol{y}\in\mathbb{C}^k$）と書け、$Q_k$ の列は一次独立なので $\boldsymbol{y}$ は $\boldsymbol{x}$ から一意に決まります。同じ定理の 2 と $\boldsymbol{r}_0=\beta\boldsymbol{q}_1=Q_{k+1}(\beta\boldsymbol{e}_1)$ から $$ \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) $$ で、$Q_{k+1}$ の列は正規直交なので $\|\boldsymbol{b}-A\boldsymbol{x}\|=\|\beta\boldsymbol{e}_1-\tilde H_k\boldsymbol{y}\|$ です。$\tilde H_k$ は上 Hessenberg で、$j\le k0$ です。$\tilde H_k$ の第 2〜$k+1$ 行は対角成分 $h_{j+1,j}\ne0$ の上三角行列をなすので、$\tilde H_k$ の列は一次独立です。 の 3 により最小解 $\boldsymbol{y}_k$ は一意で、$\boldsymbol{x}_k=\boldsymbol{x}_0+Q_k\boldsymbol{y}_k$ も一意です。 (2) $j=1,\ldots,k$ の順に、その時点の第 $j$ 列の $(j,j)$ 成分 $a$ と $(j+1,j)$ 成分 $b$ から の回転を作って第 $j$・$j+1$ 行に作用させると、$(j+1,j)$ 成分が 0 になります。$b$ は元の $h_{j+1,j}>0$ のままです（それより前の回転は第 $j+1$ 行を動かさない）から $(a,b)\ne(0,0)$ です。第 $j$ 列より左の列の第 $j$・$j+1$ 行はすでに 0 なので壊れません。こうして $\tilde H_k=Q\binom{R_k}{\boldsymbol{0}^{\mathsf{T}}}$、$Q^{\mathsf{T}}=G_k\cdots G_1$ は QR 分解で、 により $R_k$ は正則、$\boldsymbol{y}_k=R_k^{-1}(g_1,\ldots,g_k)^{\mathsf{T}}$、最小残差は $|g_{k+1}|$ です。 $\tilde H_{k+1}$ は $\tilde H_k$ の下に零の行を足し、右に新しい列 $(h_{1,k+1},\ldots,h_{k+2,k+1})^{\mathsf{T}}$ を足したものです。$G_1,\ldots,G_k$ は第 $1$〜$k+1$ 行にしか作用しないので、最初の $k$ 列は $\binom{R_k}{\boldsymbol{0}}$ のまま（第 $k+2$ 行は 0）です。新しい列に $G_1,\ldots,G_k$ を作用させ、その $(k+1,k+1)$・$(k+2,k+1)$ 成分から $G_{k+1}$ を作れば QR 分解が更新され、$\boldsymbol{g}$ も $(g_1,\ldots,g_{k+1},0)^{\mathsf{T}}$ に $G_{k+1}$ を作用させるだけです。 **実装の要点.** 反復 $k$ ごとの主な手間は、行列ベクトル積 1 回と、$\boldsymbol{q}_1,\ldots,\boldsymbol{q}_k$ との直交化（内積 $k$ 回、$O(nk)$）です。$k$ 回の合計は $O(nk^2)$、記憶は $O(nk)$ です。残差のノルム $|g_{k+1}|$ は $\boldsymbol{x}_k$ を作らずに分かるので、停止の判定は安く済み、$\boldsymbol{x}_k$ は最後に 1 回だけ作ります。浮動小数の Arnoldi では修正 Gram–Schmidt か再直交化を使います（[QR 分解](/mathematics/numerical-linear-algebra/qr-factorization) の ）。 $k=d$（Arnoldi の過程が $h_{d+1,d}=0$ で止まる）のとき、$H_d$ は正則で $\boldsymbol{x}_d=\boldsymbol{x}_0+Q_dH_d^{-1}(\beta\boldsymbol{e}_1)$ は $A\boldsymbol{x}=\boldsymbol{b}$ の厳密解であり、$\boldsymbol{x}_0+\mathcal{K}_d$ の上で残差を最小にする（$\boldsymbol{r}_d=\boldsymbol{0}$）。とくに厳密な算術では GMRES は $d\le n$ 回以内に厳密解に達し、$A$ が対角化可能なら相異なる固有値の個数以内に達する。 の 2 から $AQ_d=Q_dH_d$ です。$H_d\boldsymbol{y}=\boldsymbol{0}$ なら $AQ_d\boldsymbol{y}=\boldsymbol{0}$、$A$ は正則なので $Q_d\boldsymbol{y}=\boldsymbol{0}$、$Q_d$ の列は一次独立なので $\boldsymbol{y}=\boldsymbol{0}$ です。よって $H_d$ は正則です。$\boldsymbol{x}=\boldsymbol{x}_0+Q_d\boldsymbol{y}$ について $\boldsymbol{b}-A\boldsymbol{x}=Q_d(\beta\boldsymbol{e}_1-H_d\boldsymbol{y})$ で、$\boldsymbol{y}=H_d^{-1}\beta\boldsymbol{e}_1$ とすれば 0 になります。残差のノルムは非負なのでこれが最小です。$d\le n$ と、対角化可能なとき $d$ が相異なる固有値の個数以下であることは の 2・3 です。 止まること（breakdown）は失敗ではなく、成功です。共役勾配法と違い、GMRES の残差は**単調非増加**です（$\mathcal{K}_k\subset\mathcal{K}_{k+1}$ の上での最小化だから）。 ## 4. 収束の評価 $1\le k\le d$ について $$ \|\boldsymbol{r}_k\|=\min\bigl\{\|p(A)\boldsymbol{r}_0\|:p\in\mathcal{P}_k,\ p(0)=1\bigr\}. $$ の 1 により $\boldsymbol{x}\in\boldsymbol{x}_0+\mathcal{K}_k$ は $\boldsymbol{x}=\boldsymbol{x}_0+q(A)\boldsymbol{r}_0$（$q\in\mathcal{P}_{k-1}$）と書け、$\boldsymbol{b}-A\boldsymbol{x}=\boldsymbol{r}_0-Aq(A)\boldsymbol{r}_0=p(A)\boldsymbol{r}_0$、$p(t)=1-tq(t)$ です。逆に $p\in\mathcal{P}_k$、$p(0)=1$ なら $p(t)-1$ は定数項が 0 なので $p(t)=1-tq(t)$（$q\in\mathcal{P}_{k-1}$）と書けます。したがって $\{\boldsymbol{b}-A\boldsymbol{x}\}$ と $\{p(A)\boldsymbol{r}_0\}$ は同じ集合で、最小値（ と により存在する）は等しいです。 $A$ が正規行列なら、$1\le k\le d$ について $$ \frac{\|\boldsymbol{r}_k\|}{\|\boldsymbol{r}_0\|}\le\min_{p\in\mathcal{P}_k,\ p(0)=1}\ \max_{\lambda\in\sigma(A)}|p(\lambda)|. $$ により $A=UDU^{*}$（$U$ はユニタリ、$D=\operatorname{diag}(\lambda_1,\ldots,\lambda_n)$）です。$A^j=UD^jU^{*}$ なので $p(A)=Up(D)U^{*}$ で、$U$ と $U^{*}$ はノルムを保つ（）から $$ \|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=V\Lambda V^{-1}$（$\Lambda$ は対角）と対角化できるなら、$1\le k\le d$ について $$ \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}\|. $$ $p(A)=Vp(\Lambda)V^{-1}$ なので $\|p(A)\boldsymbol{r}_0\|\le\|V\|\|p(\Lambda)\|\|V^{-1}\|\|\boldsymbol{r}_0\|$ で、対角行列では $\|p(\Lambda)\boldsymbol{z}\|^2=\sum_i|p(\lambda_i)|^2|z_i|^2$ から $\|p(\Lambda)\|=\max_i|p(\lambda_i)|$ です。 で最小をとります。 $A=V\Lambda V^{-1}$ と対角化でき、すべての固有値が円板 $|z-c|\le\rho$（$0\le\rho<|c|$）に入るなら、$1\le k\le d$ について $\|\boldsymbol{r}_k\|/\|\boldsymbol{r}_0\|\le\kappa(V)(\rho/|c|)^k$。 $p(z)=(1-z/c)^k$ は $\mathcal{P}_k$ に属し $p(0)=1$ で、固有値 $\lambda$ について $|p(\lambda)|=|c-\lambda|^k/|c|^k\le(\rho/|c|)^k$ です。 に代入します。 前処理の目標は、この評価が効く形、すなわち「固有値が 0 から離れた小さな円板に集まり、$\kappa(V)$ も小さい」形に $A$ を変えることです。$c=1$ の近くに固有値を集めるのが典型で、$M\approx A$ なら $M^{-1}A\approx I$ です。 次の例のねらいは、固有値の位置が同じでも収束の速さが大きく違いうることを、数値で見ることです。$g$ は行列の「正規からのずれ」の強さを決める数で、$\kappa(V)$ は固有ベクトルを並べた行列の条件数です。上の収束の評価が $\kappa(V)$ を含んでいたことが、ここで効いてきます。 対角に $0.5$ から $1.5$ まで等間隔の 100 個の固有値を並べ、すぐ上の副対角に $g$ を置いた上三角行列 $A_g$ を考えます。固有値は $g$ によらず円板 $|z-1|\le0.5$ に入り、 の $(\rho/|c|)^k=0.5^k$ に $\kappa(V)$ を掛けたものが評価です。$g=0$ は対角（正規）行列です。$\boldsymbol{b}=(1,\ldots,1)^{\mathsf{T}}$、$\boldsymbol{x}_0=\boldsymbol{0}$ で、再始動なしの GMRES の相対残差 $\|\boldsymbol{r}_k\|/\|\boldsymbol{r}_0\|$ を記録します（`--` は収束して停止した後）。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 import numpy as np, scipy.sparse.linalg as sla n = 100 lam = np.linspace(0.5, 1.5, n) # 固有値は円板 |z-1| <= 0.5 の中 b = np.ones(n) gs = (0.0, 1.0, 2.0) hist = {} 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") hist[g] = res 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}") ``` ```text k 0.5^k g=0 g=1 g=2 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_{ij}=g\,v_{i+1,j}/(\lambda_j-\lambda_i)$ で書けます。浮動小数の `eig` ではこの条件数を正しく計算できないので、mpmath の 250 桁で計算します（列は単位ベクトルに正規化）。 ```python # （続き）固有ベクトル行列の条件数を mpmath の 250 桁で計算 import mpmath as mp # mpmath 1.3.0 mp.mp.dps = 250 n = 100 lm = [mp.mpf(1) / 2 + mp.mpf(i) / (n - 1) for i in range(n)] for g in (1, 2): V = mp.zeros(n, n) for j in range(n): # 上三角の固有ベクトル（閉じた形） V[j, j] = 1 for i in range(j - 1, -1, -1): V[i, j] = g * V[i + 1, j] / (lm[j] - lm[i]) nj = mp.norm(V[:, j]) for i in range(j + 1): 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)) ``` ```text g=1: kappa(V) = 7.35e+71 g=2: kappa(V) = 3.99e+101 ``` $g=0$ では残差は $0.5^k$ よりずっと速く減り、$g=1$ では遅くなって $k=60$ 以降 $1.7\times10^{-7}$ の付近で止まり、$g=2$ では $k=10$ 以降 $4\times10^{-2}$ で停滞します。$g=1$ の停滞は丸め誤差のせいではなく、厳密な GMRES の性質です。GMRES を mpmath の 40 桁で計算しても $k=40,60,80,90$ の相対残差は $6.2\times10^{-7}$、$1.8\times10^{-7}$、$1.73\times10^{-7}$、$1.73\times10^{-7}$ で、倍精度の表と一致しました。厳密な算術では残差は単調非増加で、$k=d\le100$ で 0 になる（）ので、停滞はその手前の長い区間で起きていることになります。$\kappa(V)\approx10^{72},10^{102}$ を掛けた評価は意味を持ちません。固有値が同じでも、非正規性の強さで収束はまったく違います。 **検算.** (a) 既知値: $g=0$ では固有値が実区間 $[0.5,1.5]$ にあり、 に $p(z)=T_k(2(1-z))/T_k(2)$（$T_k$ は Chebyshev 多項式、区間上で $|T_k(2(1-z))|\le1$）を入れると評価は $1/T_k(2)$ です。$T_5(2)=362$、$T_{10}(2)=262087$、$T_{20}(2)=2T_{10}(2)^2-1\approx1.37\times10^{11}$ から評価は $2.8\times10^{-3}$、$3.8\times10^{-6}$、$7.3\times10^{-12}$ で、表の $1.6\times10^{-3}$、$2.4\times10^{-6}$、$4.5\times10^{-12}$ はその内側にあります。(b) 別精度: $\kappa(V)$ は mpmath の 250 桁で計算し、350 桁でも表示桁が一致することを確かめました。(c) 次数: $g=0$ の残差の比 $\|\boldsymbol{r}_{20}\|/\|\boldsymbol{r}_{10}\|\approx1.9\times10^{-6}$ は、Chebyshev の評価の比 $T_{10}(2)/T_{20}(2)\approx1.9\times10^{-6}$ とほぼ一致します。 **任意の収束曲線（主張のみ、定理環境にしない）.** 任意の非増加な列 $1=f_0\ge f_1\ge\cdots\ge f_{n-1}>0$ と、任意の $n$ 個の 0 でない複素数に対し、それらを固有値にもち、GMRES の相対残差がちょうど $f_k$ になる行列と右辺が存在します（Greenbaum–Pták–Strakoš 1996）。固有値だけから GMRES の収束について何も言えない、という上の例の極端な形です。非正規行列の収束は、擬スペクトルや数値域で論じます（Greenbaum の教科書 Ch. 3）。この記事の証明では使いません。 **再始動と他の Krylov 法.** 記憶と直交化の手間を抑えるため、$m$ 回ごとに $\boldsymbol{x}_m$ を新しい初期値として最初からやり直すのが GMRES($m$) です。残差は単調非増加のままですが、部分空間を捨てるので停滞しうる（たとえば $A\boldsymbol{r}_0\perp\boldsymbol{r}_0$ 型の行列で 1 歩も進まない）ことが知られています。BiCGSTAB は短い漸化式で記憶を一定に保つ非対称向けの方法、MINRES は対称（不定値でもよい）な行列で GMRES と同じ最小残差の解を三項漸化式で計算する方法です（Saad の教科書 Ch. 6–7、Greenbaum の教科書 Ch. 2）。この記事では定理として扱いません。 ## 5. 前処理 正則な $M\in\mathbb{C}^{n\times n}$（**前処理行列**。$M\approx A$ で、$M\boldsymbol{z}=\boldsymbol{c}$ が安く解けるもの）について、$M^{-1}A\boldsymbol{x}=M^{-1}\boldsymbol{b}$ に反復法を使うことを**左前処理**、$AM^{-1}\boldsymbol{u}=\boldsymbol{b}$ を解いて $\boldsymbol{x}=M^{-1}\boldsymbol{u}$ とすることを**右前処理**という。 $AM^{-1}$ と初期値 $\boldsymbol{u}_0=M\boldsymbol{x}_0$ に GMRES を使った第 $k$ 近似解を $\boldsymbol{u}_k$ とし、$\boldsymbol{x}_k=M^{-1}\boldsymbol{u}_k$ とおく。このとき $\boldsymbol{x}_k$ は $$ \boldsymbol{x}_0+\mathcal{K}_k(M^{-1}A,\ M^{-1}\boldsymbol{r}_0) $$ の中で $\|\boldsymbol{b}-A\boldsymbol{x}\|$ を最小にする。一方、$M^{-1}A$、$M^{-1}\boldsymbol{b}$、$\boldsymbol{x}_0$ に GMRES を使う左前処理の第 $k$ 近似解は、同じ集合の中で $\|M^{-1}(\boldsymbol{b}-A\boldsymbol{x})\|$ を最小にする（$k$ はそれぞれの Krylov 部分空間の次数以下とする）。 右前処理: 系 $AM^{-1}\boldsymbol{u}=\boldsymbol{b}$ の初期残差は $\boldsymbol{b}-AM^{-1}\boldsymbol{u}_0=\boldsymbol{b}-A\boldsymbol{x}_0=\boldsymbol{r}_0$ で、$\boldsymbol{u}$ の残差は $\boldsymbol{b}-AM^{-1}\boldsymbol{u}=\boldsymbol{b}-A(M^{-1}\boldsymbol{u})$ です。GMRES は $\boldsymbol{u}\in\boldsymbol{u}_0+\mathcal{K}_k(AM^{-1},\boldsymbol{r}_0)$ の中でこれを最小にします。$\boldsymbol{u}\mapsto M^{-1}\boldsymbol{u}$ は全単射で、$j$ についての帰納法で $M^{-1}(AM^{-1})^j=(M^{-1}A)^jM^{-1}$ なので $$ 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）どおりこれを $\boldsymbol{x}_0+\mathcal{K}_k(M^{-1}A,M^{-1}\boldsymbol{r}_0)$ の上で最小にします。 2 つは同じ部分空間を探しますが、最小にする量が違います。左前処理では、停止判定に使う残差 $\|M^{-1}\boldsymbol{r}_k\|$ が小さくても、$M$ が悪条件なら元の残差 $\|\boldsymbol{r}_k\|$ は小さいとは限りません。右前処理では停止判定の残差がそのまま元の系の残差です。SciPy の `gmres` に `M=` を渡すと左前処理になる（前処理つき残差を最小にし、最後に元の残差で収束を確かめる）ので、右前処理にしたいときは $AM^{-1}$ を `LinearOperator` として渡します。 単位正方形で $-\Delta u+\beta(u_x+u_y)=1$、境界で $u=0$ を、格子幅 $h=1/(N+1)$ の中心差分で離散化します。行列は非対称で、$\beta h/2$（格子 Péclet 数）が大きいほど非正規性が増します。相対残差 $10^{-8}$ までの反復回数を、再始動なし（full）、GMRES(20)、`spilu`（しきい値つき不完全 LU、`drop_tol=1e-3`、`fill_factor=5`）を右前処理にした場合で比べます。参照解は `spsolve` の直接法です。 ```python # 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): # 内側の反復回数を数える it = [0] 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) return it[0], x 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}" f" {fill:4.1f} {e:.1e}") ``` ```text 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$ を 2 倍にするとほぼ 2 倍（$h^{-1}$ に比例）で、GMRES(20) は 1.6〜3.1 倍の反復を要します（1 反復の直交化は安い）。ILU の右前処理は、元の行列の 3〜3.5 倍の非零を使って反復回数をおよそ 0.3〜0.4 倍に減らします。 **検算.** (a) 既知値: 右辺と境界条件を固定した離散問題の厳密解（直接法の解）を参照値とし、右前処理の解との相対差は $10^{-10}$ の程度で、相対残差 $10^{-8}$ で止めた解として妥当な大きさです。(b) 別手法: 参照値は反復法と独立な疎 LU 分解（`spsolve`）で計算しました。(c) 次数: full GMRES の反復回数の比は $\beta=10$ で $181/91\approx1.99$、$\beta=100$ で $131/71\approx1.85$ と、$h^{-1}$ の増え方に近い値です。 ## 6. 不完全 Cholesky 分解 対称正定値の系には、Cholesky 分解 $A=LL^{\mathsf{T}}$（）の fill-in を捨てた近似 $M=LL^{\mathsf{T}}$ を共役勾配法の前処理に使うのが標準です。fill-in を捨てると正定値性が保たれる保証がなくなり、途中で平方根の中が負になって破綻しえます。ここでは、それが起きない行列のクラスを特定します。 非対角成分がすべて $0$ 以下の実正方行列を **Z 行列**という。正定値な対称 Z 行列を **Stieltjes 行列**（対称な M 行列）という。実行列 $N,N'$ について $N\ge0$ や $N'\le N$ は成分ごとの不等式とする。 離散 Laplace 行列（2 次元・3 次元を含む）は Stieltjes 行列の代表例です。 次の補題は、成分を小さくした非負行列のスペクトル半径は大きくならない、という事実です。IC の算法で成分を捨てたあとも正定値のままか、を確かめる場面で使います。証明の方針は、$N/t$ の冪が 0 に近づくなら、それ以下の成分をもつ $N'/t$ の冪も 0 に近づく、という比較です。 $N,N'\in\mathbb{R}^{n\times n}$ が $0\le N'\le N$ を満たすなら $\rho(N')\le\rho(N)$。 $t>\rho(N)$ を任意にとります。成分が非負の行列の積は各因子について成分ごとに単調なので、$k$ についての帰納法で $0\le(N'/t)^k\le(N/t)^k$ です。$\rho(N/t)<1$ なので により $(N/t)^k\to0$、はさみうちで $(N'/t)^k\to0$、同じ系から $\rho(N'/t)<1$、すなわち $\rho(N')\rho(N)$ は任意なので $\rho(N')\le\rho(N)$ です。 $A$ を対称 Z 行列、$s\ge\max_ia_{ii}$、$N=sI-A$ とすると $N\ge0$ で、$A$ が正定値 $\iff\rho(N) $N$ の対角は $s-a_{ii}\ge0$、非対角は $-a_{ij}\ge0$ です。$N$ は実対称なので、 により実の正規直交固有ベクトルの基底があり、固有値 $\mu_1\ge\cdots\ge\mu_n$ は実数、$A$ の固有値は $s-\mu_i$ です。同じ定理の 2 により $A$ が正定値 $\iff$ 全 $i$ で $s-\mu_i>0\iff\mu_1 から $$ |\mu|=|\boldsymbol{x}^{\mathsf{T}}N\boldsymbol{x}|\le|\boldsymbol{x}|^{\mathsf{T}}N|\boldsymbol{x}|\le\mu_1 $$ です（$|\boldsymbol{x}|$ も単位ベクトル）。よって $\rho(N)\le\mu_1\le\rho(N)$ です。 **IC の算法（ねらい）.** Cholesky 分解を 1 列ずつ進めると、各段で「まだ分解していない部分」を表す小さな行列（Schur 補行列）が残ります。IC は、その行列の成分のうち、あらかじめ決めた非零パターンの外にあるものを 0 に捨てて次の段へ進みます。捨てると正定値性が壊れて平方根の中が負になりかねないので、次の定理は、Stieltjes 行列ならそれが起きないことを段ごとの帰納法で示します。 **IC の算法.** 非零パターン $S\subset\{1,\ldots,n\}^2$ を、対称（$(i,j)\in S\Rightarrow(j,i)\in S$）で対角を含み、$A$ の非零成分の位置をすべて含むものとします（$S$ を $A$ の非零パターンそのものにとるのが IC(0)）。$B^{(1)}=A$ とし、$k=1,\ldots,n$ について、$n-k+1$ 次の対称行列 $B^{(k)}$（添字は $k,\ldots,n$）を $$ 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) $$ と分け、Schur 補行列 $S_k=C_k-\boldsymbol{w}_k\boldsymbol{w}_k^{\mathsf{T}}/d_k$ の $(i,j)\notin S$ の成分を $0$ に置き換えたものを $B^{(k+1)}$ とします。置き換えの差を $F_k$（$B^{(k+1)}=S_k+F_k$）とします。 $A$ が Stieltjes 行列なら、上の算法のすべての段で $d_k>0$ であり、各 $B^{(k)}$ は Stieltjes 行列である。得られる下三角行列 $L=(l_{ij})$ は対角成分が正で、$(i,j)\notin S$ では $l_{ij}=0$ であり、 $$ LL^{\mathsf{T}}=A+R,\qquad R=\sum_kF_k\ (\text{位置を埋め込んだ和})\ \ge0,\qquad(i,j)\in S\text{ で }r_{ij}=0 $$ が成り立つ。とくに $M=LL^{\mathsf{T}}$ は正定値である（Meijerink–van der Vorst の対称な場合）。 **帰納法の主張.** $B^{(k)}$ は Stieltjes 行列で、$(i,j)\notin S$ の成分は 0 である。$k=1$ では仮定どおりです（$S$ は $A$ の非零パターンを含む）。 **$d_k>0$.** 正定値行列の対角成分は $\boldsymbol{e}_1^{\mathsf{T}}B^{(k)}\boldsymbol{e}_1>0$ です。 **Schur 補行列.** 行と列の順序を入れ替えた $\begin{pmatrix}C_k&\boldsymbol{w}_k\\\boldsymbol{w}_k^{\mathsf{T}}&d_k\end{pmatrix}$ は $B^{(k)}$ と置換行列で合同なので正定値です（$\boldsymbol{x}^{\mathsf{T}}P^{\mathsf{T}}BP\boldsymbol{x}=(P\boldsymbol{x})^{\mathsf{T}}B(P\boldsymbol{x})$）。 の 3（右下の $d_k>0$ は正定値）により $S_k=C_k-\boldsymbol{w}_kd_k^{-1}\boldsymbol{w}_k^{\mathsf{T}}$ は正定値です。非対角成分は、$(\boldsymbol{w}_k)_i,(\boldsymbol{w}_k)_j\le0$ から $(\boldsymbol{w}_k)_i(\boldsymbol{w}_k)_j\ge0$ なので $(S_k)_{ij}\le(C_k)_{ij}\le0$ です。よって $S_k$ は Stieltjes 行列です。 **捨てる操作.** $s=\max_i(S_k)_{ii}$ とし、$S_k=sI-N$ と書くと、 により $N\ge0$、$\rho(N) で $\rho(N')\le\rho(N) で $B^{(k+1)}$ は正定値で、対称 Z 行列なので Stieltjes 行列です。$(i,j)\notin S$ の成分が 0 であることは作り方のとおりで、帰納法が進みます。 **分解の式.** $(\boldsymbol{w}_k)_i$ は $B^{(k)}$ の $(i,k)$ 成分なので、$(i,k)\notin S$ なら $0$ で、$l_{ik}=0$ です。$\boldsymbol{l}_k=(l_{kk},\ldots,l_{nk})^{\mathsf{T}}$ とおくと $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}$ です。これを $n\times n$ に埋め込んで $k=1,\ldots,n$ について足すと $A=\sum_k\boldsymbol{l}_k\boldsymbol{l}_k^{\mathsf{T}}-\sum_kF_k=LL^{\mathsf{T}}-R$ です。$F_k$ の 0 でない成分は $(i,j)\notin S$ の位置にあり、値は $-(S_k)_{ij}\ge0$ です。$L$ は対角が正の下三角なので正則で、$\boldsymbol{x}^{\mathsf{T}}LL^{\mathsf{T}}\boldsymbol{x}=\|L^{\mathsf{T}}\boldsymbol{x}\|^2>0$（$\boldsymbol{x}\ne\boldsymbol{0}$）です。 $S$ を全成分にとれば何も捨てないので、これは Cholesky 分解そのものです。正定値でも Z 行列でない行列では、IC(0) が破綻する例が知られています（Saad の教科書 Ch. 10）。 2 次元の離散 Laplace 行列（$\operatorname{tridiag}(-1,2,-1)$ の Kronecker 和、Stieltjes 行列）について、前処理なしの共役勾配法と、自作の IC(0)（上の算法を行ごとに並べ替えたもの）を前処理にした共役勾配法の反復回数を、相対残差 $10^{-8}$ まで比べます。参照解は `spsolve` です。 行ごとの版が上の算法と同じ $L$ を与えることは、成分ごとに見ると分かります。$(i,j)\in S$（$i\ge j$）の成分は、第 $kj) $$ で、$(i,j)\notin S$ では $l_{ij}=0$ です（ の証明の「分解の式」）。下のコードはこの式を $i$ の小さい順、各行の中で $j$ の小さい順に計算するもので、右辺に現れる $l$ はすべて計算済みです。 ```python # 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] = {列: 値} for i in range(n): 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()] r, c, v = zip(*t) 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) cnt.append(it[0]) 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}" f" {e:.1e}") ``` ```text N CG CG+IC(0) ratio err(IC) 16 28 17 1.65 2.1e-10 32 59 29 2.03 2.5e-10 64 119 52 2.29 2.7e-10 128 239 100 2.39 1.0e-10 ``` 前処理なしの反復回数は $N$ を 2 倍にするとほぼ 2 倍（$\sqrt{\kappa}\propto h^{-1}$）で、IC(0) でも増え方は $h^{-1}$ に近いまま（比 $1.7\to1.8\to1.9$）、反復回数が 0.4〜0.6 倍になるだけです。IC(0) は条件数の次数 $h^{-2}$ を変えず、定数を改善する前処理だということが読み取れます。 **検算.** (a) 既知値: のとおり、どの格子でも平方根の中は正で分解は最後まで進みました。前処理つきの解は `spsolve` の参照解と相対 $10^{-10}$ の程度で一致します。(b) 別手法: 前処理の有無で同じ解に収束しています。(c) 次数: 前処理なしの反復回数の比 $59/28,119/59,239/119$ は $2.1,2.0,2.0$ で、共役勾配法の評価（反復回数 $\propto\sqrt{\kappa}$）と [ノルムの章](/mathematics/numerical-linear-algebra/norms-and-conditioning) の $\kappa\approx4/(\pi^2h^2)$ から予測される $h^{-1}$ の増え方に一致します。 **修正不完全 Cholesky と格子に依らない前処理.** 捨てる成分を対角に足し戻す修正不完全 Cholesky 分解（MIC）は、モデル問題で $\kappa(M^{-1}A)$ を $O(h^{-1})$ に下げ、反復回数を $h^{-1/2}$ に比例させます（Gustafsson 1978。Saad の教科書 Ch. 10）。多重格子法と領域分割法は、適切に構成すれば反復回数が格子幅に依らない前処理（あるいは解法）になります（Trottenberg–Oosterlee–Schüller）。どちらも理論が長いので、この記事では証明しません。 ## 7. 疎行列の反復解法の実務 反復法が行列に要求するのは行列ベクトル積だけです。格納形式と演算量は [Cholesky 分解と疎行列の直接法](/mathematics/numerical-linear-algebra/cholesky-and-sparse-direct) の と で扱いました。反復法の観点でまとめると次のとおりです（$\mathrm{nnz}$ は非零の個数、$k$ は反復回数）。 | 形式・道具 | 記憶 | 行列ベクトル積 | 反復法での使いどころ | |---|---|---|---| | COO（`coo_array`） | $3\,\mathrm{nnz}$ | $2\,\mathrm{nnz}$ | 組み立て。計算の前に CSR・CSC に変換する | | CSR（`csr_array`） | $2\,\mathrm{nnz}+n+1$ | $2\,\mathrm{nnz}$ | $A\boldsymbol{x}$ の標準。行ごとの処理（IC・ILU の行版） | | CSC（`csc_array`） | $2\,\mathrm{nnz}+n+1$ | $2\,\mathrm{nnz}$ | `splu`・`spilu`・`spsolve` の入力 | | `LinearOperator` | 関数だけ | 関数の手間 | 行列を作らない演算子、$AM^{-1}$ の合成、前処理 $M^{-1}$ | | GMRES（再始動なし） | $O(nk)$ | $k$ 回 | 直交化 $O(nk^2)$。$k$ が大きいと再始動か前処理が要る | | 共役勾配法 | $O(n)$ | $k$ 回 | 正定値。短い漸化式で直交化は不要 | `LinearOperator((n, n), matvec=f)` は、ベクトルを受け取って $A\boldsymbol{v}$ を返す関数 `f` だけで行列の代わりになります。上の では $AM^{-1}$ を、 では $(LL^{\mathsf{T}})^{-1}$ をこの形で渡しました。Jacobian を作らずに差分で $J\boldsymbol{v}$ を近似する Newton–Krylov 法や、FFT で畳み込みを計算する演算子も同じ形で反復法に渡せます。時間の計測値は環境に依存するので、比べるときは反復回数と、1 反復の演算量（$2\,\mathrm{nnz}$ と前処理の適用の手間）を掛けて見積もるのがよいと思います。 ## 8. 演習 $A$ が対角化可能で、相異なる固有値が 2 つ $\lambda_1,\lambda_2$（ともに 0 でない）しかないとき、GMRES は高々 2 回で厳密解に達することを、残差多項式を具体的に与えて示せ。 $p(z)=(1-z/\lambda_1)(1-z/\lambda_2)$ は $\mathcal{P}_2$ に属し $p(0)=1$、$p(\lambda_1)=p(\lambda_2)=0$ です。 から $\|\boldsymbol{r}_2\|\le\kappa(V)\cdot0=0$（$d\ge2$ のとき）。$d=1$ なら により 1 回で厳密解です。 $A$ を巡回置換行列（$A\boldsymbol{e}_i=\boldsymbol{e}_{i+1}$、$A\boldsymbol{e}_n=\boldsymbol{e}_1$）、$\boldsymbol{b}=\boldsymbol{e}_1$、$\boldsymbol{x}_0=\boldsymbol{0}$ とする。$k と矛盾しないことを説明せよ。 $\mathcal{K}_k(A,\boldsymbol{e}_1)=\operatorname{span}\{\boldsymbol{e}_1,\ldots,\boldsymbol{e}_k\}$ で、$A\mathcal{K}_k=\operatorname{span}\{\boldsymbol{e}_2,\ldots,\boldsymbol{e}_{k+1}\}$ は $k (1) $A=\begin{pmatrix}1&-2\\-2&1\end{pmatrix}$ は対称 Z 行列だが、Cholesky 分解（$S$ を全成分にとった IC）の第 2 段で $d_2<0$ となることを確かめ、 のどの仮定が欠けているか答えよ。(2) 同じ $A$ について の $N$ と $\rho(N)$ を求めよ。 (1) $d_1=1$、$\boldsymbol{w}_1=(-2)$ で、$d_2=1-(-2)^2/1=-3<0$ です。$A$ の固有値は $1\pm2$、すなわち $3,-1$ で、正定値でありません。Stieltjes 行列の仮定のうち正定値性が欠けています。(2) $s=1$ とすると $N=I-A=\begin{pmatrix}0&2\\2&0\end{pmatrix}$、固有値 $\pm2$ で $\rho(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.