# Krylov 部分空間法 I：Lanczos 法と共役勾配法の最適性と Chebyshev 多項式による収束評価

> Krylov 部分空間と Arnoldi の関係式、Lanczos の三項漸化式と Ritz 値の最適性・交互性、Kaniel–Paige の評価、共役勾配法の最適性と √κ に比例する収束評価、最急降下法、前処理つき共役勾配法の同値性を証明し、離散 Laplace 方程式で確かめる。
> https://rikai.mugen-giken.com/mathematics/numerical-linear-algebra/lanczos-and-conjugate-gradient

## 0. この記事の要点

未知数が何万もある連立一次方程式 $A\boldsymbol{x} = \boldsymbol{b}$ や、大きな行列の固有値を、行列とベクトルの掛け算だけで効率よく求める方法の話です。ベクトルに行列を繰り返し掛けると新しいベクトルが次々にできますが、基本の発想は、それらを捨てずに全部使って最もよい近似を選ぶことです。読み終えると、共役勾配法がなぜ少ない回数で解けるのか、反復回数が行列の「性質の悪さ」を表す数の平方根に比例するとはどういうことかが分かります。

- $\boldsymbol{b}, A\boldsymbol{b}, \ldots, A^{k-1}\boldsymbol{b}$ の張る **Krylov 部分空間** $\mathcal{K}_k(A, \boldsymbol{b})$ は「$\deg p < k$ の多項式 $p$ による $p(A)\boldsymbol{b}$ の全体」です。[冪乗法](/mathematics/numerical-linear-algebra/power-and-inverse-iteration#def-power-method)は最後のベクトル $A^{k-1}\boldsymbol{b}$ しか使いませんが、Krylov 部分空間法はこの空間全体から最良のものを選びます。
- [正規直交基底](/mathematics/linear-algebra/inner-product-spaces#def-orthonormal)を作る Arnoldi の過程（$\mathcal{K}_k$ の基底を 1 本ずつ直交化しながら作る手順）は $AQ_k = Q_{k+1}\tilde{H}_k$ を満たし、$A$ が Hermite 行列なら $H_k$ は三重対角になって、3 項の漸化式（**Lanczos 法**）で済みます。
- $H_k$ の固有値（**Ritz 値**。$\mathcal{K}_k$ の中で見つかる固有値の近似値）の最大は $\mathcal{K}_k$ 上の Rayleigh 商の最大で、$k$ とともに単調に増えて $\lambda_1$ に近づきます。誤差は $(\lambda_1 - \lambda_n)\tan^2\angle(\boldsymbol{b}, \boldsymbol{u}_1)/T_{k-1}(1 + 2\gamma)^2$ 以下です（**Kaniel–Paige の評価**）。
- 正定値行列の**共役勾配法**は、$\boldsymbol{x}_0 + \mathcal{K}_k$ の上で誤差の $A$ ノルム（誤差 $\boldsymbol{e}$ の大きさを $(\boldsymbol{e}^{*}A\boldsymbol{e})^{1/2}$ で測ったもの）を最小にします。そのため $\|\boldsymbol{e}_k\|_A \le 2\bigl((\sqrt{\kappa} - 1)/(\sqrt{\kappa} + 1)\bigr)^k\|\boldsymbol{e}_0\|_A$ が成り立ち、反復回数は $\sqrt{\kappa}$ に比例します。最急降下法の $(\kappa - 1)/(\kappa + 1)$ より桁違いに速く、厳密な算術では相異なる固有値の個数以下の回数で終わります。
- 共役勾配法の残差は Lanczos ベクトルの定数倍で、共役勾配法は Lanczos 法で作った三重対角の系を解くことと同じです。前処理つき共役勾配法（前処理は、方程式を解きやすい形に変える下ごしらえ）は $M^{-1/2}AM^{-1/2}$ への共役勾配法と同じ反復で、収束は $\kappa(M^{-1}A)$ で決まります。

<Aside type="tip" title="この記事で使う事前知識">

- **正定値**: 零でないどんなベクトル $\boldsymbol{x}$ についても $\boldsymbol{x}^{*}A\boldsymbol{x} > 0$ となる Hermite 行列の性質。→ [正定値](/mathematics/linear-algebra/spectral-theorem#def-positive-definite)
- **条件数**: 解が入力の誤差にどれだけ敏感かを表す数。この記事では正定値行列の最大固有値と最小固有値の比 $\kappa$ として使います。→ [条件数](/mathematics/numerical-linear-algebra/norms-and-conditioning#def-condition-number)
- **冪乗法**: ベクトルに $A$ を繰り返し掛けて、絶対値最大の固有値の固有ベクトルに近づける方法。→ [冪乗法](/mathematics/numerical-linear-algebra/power-and-inverse-iteration#def-power-method)
- **Rayleigh 商**: $R_A(\boldsymbol{x}) = \boldsymbol{x}^{*}A\boldsymbol{x}/\boldsymbol{x}^{*}\boldsymbol{x}$。固有ベクトルでは固有値に一致する量。→ [固有値の計算 I](/mathematics/numerical-linear-algebra/power-and-inverse-iteration)
- **正規直交基底**: 長さが 1 で互いに直交するベクトルからなる基底。→ [正規直交基底](/mathematics/linear-algebra/inner-product-spaces#def-orthonormal)
- **定常反復法**: 同じ行列 $G$ を誤差に繰り返し掛けて解に近づける反復（Jacobi 法など）。→ [定常反復法](/mathematics/numerical-linear-algebra/stationary-iterative-methods)

</Aside>

## 1. 動機

[定常反復法の章](/mathematics/numerical-linear-algebra/stationary-iterative-methods) では、2 次元の離散 Laplace 方程式で Jacobi 法の反復回数が $h^{-2}$ に比例し、最適な $\omega$ を知っていても SOR 法は $h^{-1}$ が限界であることを見ました。しかも $\omega$ の最適値は Jacobi 法のスペクトル半径という、ふつうは分からない量に依存していました。

定常反復の $k$ 回目の誤差は $G^k\boldsymbol{e}_0$ です。$G = I - M^{-1}A$ なので、これは「$p(0) = 1$ を満たす特定の $k$ 次多項式 $p(t) = (1 - t)^k$ に $M^{-1}A$ を代入して $\boldsymbol{e}_0$ に掛けたもの」です。それなら、$p(0) = 1$ の $k$ 次多項式の中から**最良の**ものを選べばもっと速くなるはずです。この発想を、パラメータを一切知らずに実行するのが共役勾配法です。

固有値の計算でも事情は同じです。[冪乗法](/mathematics/numerical-linear-algebra/power-and-inverse-iteration) は $A^{k-1}\boldsymbol{b}$ だけを使い、それまでに作った $\boldsymbol{b}, A\boldsymbol{b}, \ldots$ を捨てていました。これらが張る空間の中で Rayleigh 商を最大にすれば、同じ行列ベクトル積の回数ではるかによい近似が得られます。これが Lanczos 法です。

2 つの算法は同じ Krylov 部分空間の上で動き、収束の評価はどちらも「区間の外で急に大きくなる多項式」、すなわち Chebyshev 多項式で与えられます。

## 2. Krylov 部分空間

この章でも $\mathbb{C}^n$ の標準内積 $\boldsymbol{y}^{*}\boldsymbol{x}$、ノルム $\|\cdot\|$ と、[冪乗法の章](/mathematics/numerical-linear-algebra/power-and-inverse-iteration) の角（<Ref to="mathematics/numerical-linear-algebra/power-and-inverse-iteration#def-angle" />）と Rayleigh 商 $R_A(\boldsymbol{x})$ を使います。$\mathcal{P}_k$ は次数 $k$ 以下の多項式全体です。

<Definition id="def-krylov" title="Krylov 部分空間">
$A \in \mathbb{C}^{n\times n}$、$\boldsymbol{b} \in \mathbb{C}^n$、$\boldsymbol{b} \ne \boldsymbol{0}$ とする。$k \ge 1$ について
$$
\mathcal{K}_k(A, \boldsymbol{b}) = \operatorname{span}\{\boldsymbol{b}, A\boldsymbol{b}, \ldots, A^{k-1}\boldsymbol{b}\}
$$
を $k$ 次の **Krylov 部分空間**という。$A^{d}\boldsymbol{b} \in \mathcal{K}_d(A, \boldsymbol{b})$ となる最小の $d \ge 1$ を $\boldsymbol{b}$ の**次数**（grade）という。
</Definition>

<Proposition id="prop-krylov-polynomial" title="Krylov 部分空間の基本性質">
<Ref to="def-krylov" /> の記号で次が成り立つ。

1. $\mathcal{K}_k(A, \boldsymbol{b}) = \{p(A)\boldsymbol{b} : p \in \mathcal{P}_{k-1}\}$。
2. $\dim\mathcal{K}_k(A, \boldsymbol{b}) = \min(k, d)$ であり、$k \ge d$ なら $\mathcal{K}_k = \mathcal{K}_d$ で、$\mathcal{K}_d$ は $A$ で不変である。
3. $A$ が対角化可能で相異なる固有値が $s$ 個なら $d \le s$ である。
</Proposition>

<Proof of="prop-krylov-polynomial">
1. $p(t) = \sum_{j<k}c_jt^j$ なら $p(A)\boldsymbol{b} = \sum_jc_jA^j\boldsymbol{b}$ なので、定義の言い換えです。

2. $k < d$ では $d$ の最小性から $A^{k}\boldsymbol{b} \notin \mathcal{K}_k$ なので、帰納法で $\boldsymbol{b}, \ldots, A^{k}\boldsymbol{b}$ は一次独立、つまり $\dim\mathcal{K}_{k+1} = k + 1$ です（$k + 1 \le d$）。$A^d\boldsymbol{b} \in \mathcal{K}_d$ なら $A\mathcal{K}_d = \operatorname{span}\{A\boldsymbol{b}, \ldots, A^d\boldsymbol{b}\} \subset \mathcal{K}_d$ で、$\mathcal{K}_d$ は $A$ で不変です。よって $A^j\boldsymbol{b} \in \mathcal{K}_d$（$j \ge d$）で、$k \ge d$ なら $\mathcal{K}_k = \mathcal{K}_d$ です。

3. 対角化可能なので固有ベクトルの基底があり（<Ref to="mathematics/linear-algebra/diagonalization-and-jordan-form#thm-diagonalizable-basis" />）、相異なる固有値を $\mu_1, \ldots, \mu_s$ として $q(t) = \prod_{i=1}^{s}(t - \mu_i)$ とおくと、$q(A)$ はすべての固有ベクトルを $0$ に写すので $q(A) = O$ です。$q$ はモニックな $s$ 次式なので $A^s\boldsymbol{b} = A^s\boldsymbol{b} - q(A)\boldsymbol{b} \in \mathcal{K}_s$ で、$d \le s$ です。
</Proof>

## 3. Arnoldi の関係式と Lanczos の三項漸化式

$\boldsymbol{b}, A\boldsymbol{b}, A^2\boldsymbol{b}, \ldots$ をそのまま基底に使うと、冪乗法と同じ理由で後のベクトルほど $\boldsymbol{v}_1$ の方向にそろい、数値的に一次従属に近づきます。そこで 1 本ずつ $A$ を掛けては直交化します。

<Definition id="def-arnoldi" title="Arnoldi の過程">
$\boldsymbol{q}_1 = \boldsymbol{b}/\|\boldsymbol{b}\|$ とし、$j = 1, 2, \ldots$ について
$$
h_{ij} = \boldsymbol{q}_i^{*}A\boldsymbol{q}_j\ (i \le j),\qquad \boldsymbol{w}_j = A\boldsymbol{q}_j - \sum_{i=1}^{j}h_{ij}\boldsymbol{q}_i,\qquad h_{j+1,j} = \|\boldsymbol{w}_j\|
$$
を計算し、$h_{j+1,j} \ne 0$ なら $\boldsymbol{q}_{j+1} = \boldsymbol{w}_j/h_{j+1,j}$ とする。$h_{j+1,j} = 0$ なら停止する。$Q_k = (\boldsymbol{q}_1\ \cdots\ \boldsymbol{q}_k)$、$H_k = (h_{ij})_{1\le i,j\le k}$（$i > j + 1$ では $h_{ij} = 0$）、$\tilde{H}_k$ を $H_k$ の下に行 $(0, \ldots, 0, h_{k+1,k})$ を足した $(k+1)\times k$ 行列とする。
</Definition>

<Theorem id="thm-arnoldi-relation" title="Arnoldi の関係式">
$d$ を $\boldsymbol{b}$ の次数とする。Arnoldi の過程は $j = 1, \ldots, d - 1$ では停止せず（$h_{j+1,j} > 0$）、$j = d$ で停止する（$h_{d+1,d} = 0$）。$1 \le k \le d$ について次が成り立つ。

1. $\boldsymbol{q}_1, \ldots, \boldsymbol{q}_k$ は正規直交で $\mathcal{K}_k(A, \boldsymbol{b})$ を張る。
2. $k < d$ なら $AQ_k = Q_{k+1}\tilde{H}_k$、$k = d$ なら $AQ_d = Q_dH_d$。
3. $H_k = Q_k^{*}AQ_k$。
</Theorem>

証明は $j$ についての帰納法です。各段で「$A\boldsymbol{q}_j$ はまだ $\mathcal{K}_j$ に入っていないので、$\mathcal{K}_j$ への成分を引くと零でない新しい方向が残り、それを正規化したものが $\boldsymbol{q}_{j+1}$ になる」ことを確かめ、最後に $j = d$ で新しい方向が残らなくなることを見ます。

<Proof of="thm-arnoldi-relation">
$j$ についての帰納法で、$j < d$ のとき「$\boldsymbol{q}_1, \ldots, \boldsymbol{q}_j$ は正規直交で $\mathcal{K}_j$ を張り、$h_{j+1,j} > 0$ で、$\boldsymbol{q}_1, \ldots, \boldsymbol{q}_{j+1}$ は正規直交で $\mathcal{K}_{j+1}$ を張る」ことを示します。$\mathcal{K}_1 = \operatorname{span}\{\boldsymbol{q}_1\}$ は明らかです。

$\boldsymbol{q}_1, \ldots, \boldsymbol{q}_j$ が正規直交で $\mathcal{K}_j$ を張るとします。$\boldsymbol{q}_j \in \mathcal{K}_j$ ですが、$j \ge 2$ では $\boldsymbol{q}_j \perp \mathcal{K}_{j-1}$ なので $\boldsymbol{q}_j \notin \mathcal{K}_{j-1}$ です。よって $\boldsymbol{q}_j = p(A)\boldsymbol{b}$ の $p$ はちょうど $j - 1$ 次で（<Ref to="prop-krylov-polynomial" /> の 1）、$A\boldsymbol{q}_j = Ap(A)\boldsymbol{b}$ を $\boldsymbol{b}, A\boldsymbol{b}, \ldots$ で展開したときの $A^j\boldsymbol{b}$ の係数は $0$ でありません。$j < d$ では $\boldsymbol{b}, \ldots, A^j\boldsymbol{b}$ は一次独立（<Ref to="prop-krylov-polynomial" /> の 2）なので $A\boldsymbol{q}_j \notin \mathcal{K}_j$ です。$\boldsymbol{w}_j$ は $A\boldsymbol{q}_j$ から $\mathcal{K}_j$ への直交射影を引いたもの（<Ref to="mathematics/linear-algebra/inner-product-spaces#thm-orthogonal-decomposition" />）なので $\boldsymbol{w}_j \ne \boldsymbol{0}$、$\boldsymbol{w}_j \perp \mathcal{K}_j$ で、$\boldsymbol{q}_{j+1}$ は $\boldsymbol{q}_1, \ldots, \boldsymbol{q}_j$ と直交する単位ベクトルです。$\boldsymbol{q}_{j+1} \in \operatorname{span}\{A\boldsymbol{q}_j, \mathcal{K}_j\} \subset \mathcal{K}_{j+1}$ で、次元を比べて $\operatorname{span}\{\boldsymbol{q}_1, \ldots, \boldsymbol{q}_{j+1}\} = \mathcal{K}_{j+1}$ です。

$j = d$ では $\mathcal{K}_d$ が $A$ で不変（<Ref to="prop-krylov-polynomial" /> の 2）なので $A\boldsymbol{q}_d \in \mathcal{K}_d$ で、$\boldsymbol{w}_d = \boldsymbol{0}$、つまり $h_{d+1,d} = 0$ です。

定義を書き直すと $A\boldsymbol{q}_j = \sum_{i=1}^{j+1}h_{ij}\boldsymbol{q}_i$（$j = d$ では和は $i \le d$ まで）で、これを列に並べたものが 2 です。2 の左から $Q_k^{*}$ を掛けると、$Q_k^{*}Q_{k+1} = (I_k\ \ \boldsymbol{0})$ から $Q_k^{*}AQ_k = H_k$ です。
</Proof>

<Corollary id="cor-lanczos-three-term" title="Lanczos の三項漸化式">
$A$ が Hermite 行列なら、$H_k$ は実の三重対角行列
$$
H_k = T_k = \begin{pmatrix} \alpha_1 & \beta_1 & & \\ \beta_1 & \alpha_2 & \ddots & \\ & \ddots & \ddots & \beta_{k-1} \\ & & \beta_{k-1} & \alpha_k \end{pmatrix},\qquad \alpha_j = \boldsymbol{q}_j^{*}A\boldsymbol{q}_j \in \mathbb{R},\quad \beta_j = h_{j+1,j} > 0
$$
であり、Arnoldi の過程は
$$
\beta_j\boldsymbol{q}_{j+1} = A\boldsymbol{q}_j - \alpha_j\boldsymbol{q}_j - \beta_{j-1}\boldsymbol{q}_{j-1}\qquad(\beta_0\boldsymbol{q}_0 = \boldsymbol{0})
$$
に簡約される（**Lanczos 法**）。
</Corollary>

証明の中心は、上 Hessenberg 行列（対角線のすぐ下より下がすべて $0$ の行列）が Hermite でもあれば三重対角になる、という点です。

<Proof of="cor-lanczos-three-term">
<Ref to="thm-arnoldi-relation" /> の 3 から $H_k^{*} = Q_k^{*}A^{*}Q_k = H_k$ で、$H_k$ は Hermite 行列です。$H_k$ は上 Hessenberg（$i > j + 1$ で $0$）なので、Hermite 性から $j > i + 1$ でも $0$、すなわち三重対角です。対角成分は Hermite 行列の対角なので実数、副対角は $h_{ij} = \overline{h_{ji}}$ と $h_{j+1,j} = \|\boldsymbol{w}_j\| > 0$ から $h_{j,j+1} = h_{j+1,j}$ で実数です。漸化式は Arnoldi の $\boldsymbol{w}_j$ の和の中で $h_{ij} = 0$（$i < j - 1$）の項を落としたものです。
</Proof>

1 回の Lanczos の反復は、行列ベクトル積 1 回とベクトル演算 $O(n)$ だけで済み、保持するベクトルも 3 本です。Arnoldi の過程は $j$ 本すべてとの内積が要るので、反復が進むほど重くなります。この差が、Hermite 行列の Krylov 部分空間法（Lanczos 法・共役勾配法）が一般の行列の方法（Arnoldi 法・GMRES）より軽い理由です。

## 4. Ritz 値と Lanczos 法の収束

以下この節では $A$ を $n$ 次 Hermite 行列とし、固有値を $\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_n$、対応する正規直交な固有ベクトルを $\boldsymbol{u}_1, \ldots, \boldsymbol{u}_n$ とします（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-spectral" />）。

<Definition id="def-ritz-value" title="Ritz 値">
$1 \le k \le d$ について、三重対角行列 $T_k = Q_k^{*}AQ_k$ の固有値 $\theta_1^{(k)} \ge \cdots \ge \theta_k^{(k)}$ を $\mathcal{K}_k(A, \boldsymbol{b})$ に関する **Ritz 値**という。$T_k\boldsymbol{s} = \theta\boldsymbol{s}$ のとき $Q_k\boldsymbol{s}$ を **Ritz ベクトル**という。
</Definition>

<Theorem id="thm-ritz-optimality" title="Ritz 値の最適性と交互性">
$1 \le k \le d$ とする。

1. $\theta_1^{(k)} = \max\{R_A(\boldsymbol{x}) : \boldsymbol{x} \in \mathcal{K}_k,\ \boldsymbol{x} \ne \boldsymbol{0}\}$、$\theta_k^{(k)} = \min\{R_A(\boldsymbol{x}) : \boldsymbol{x} \in \mathcal{K}_k,\ \boldsymbol{x} \ne \boldsymbol{0}\}$。
2. $\lambda_i \ge \theta_i^{(k)} \ge \lambda_{i+n-k}$（$1 \le i \le k$）。
3. $k + 1 \le d$ なら $\theta_i^{(k+1)} \ge \theta_i^{(k)} \ge \theta_{i+1}^{(k+1)}$（$1 \le i \le k$）。とくに $\theta_1^{(k)}$ は $k$ について単調に増加して $\lambda_1$ を超えず、$\theta_k^{(k)}$ は単調に減少して $\lambda_n$ を下回らない。
4. $k = d$ なら Ritz 値はすべて $A$ の固有値である。
</Theorem>

<Proof of="thm-ritz-optimality">
1. $\boldsymbol{x} \in \mathcal{K}_k$ は $\boldsymbol{x} = Q_k\boldsymbol{y}$ と一意に書け、$Q_k$ の列は正規直交なので $\|\boldsymbol{x}\| = \|\boldsymbol{y}\|$、$\boldsymbol{x}^{*}A\boldsymbol{x} = \boldsymbol{y}^{*}T_k\boldsymbol{y}$ です。よって $R_A(\boldsymbol{x}) = R_{T_k}(\boldsymbol{y})$ で、$\boldsymbol{y}$ は $\mathbb{C}^k \setminus \{\boldsymbol{0}\}$ 全体を動きます。Hermite 行列 $T_k$ の Rayleigh 商の最大・最小は最大・最小の固有値です（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-rayleigh" /> の 1）。

2. Cauchy の交互定理（<Ref to="mathematics/linear-algebra/matrix-inequalities#thm-interlacing" />）を $P = Q_k$（$Q_k^{*}Q_k = I_k$）に使えばよい。

3. $T_k$ は $T_{k+1}$ の左上の $k$ 次主部分行列なので、交互定理を $T_{k+1}$ と $P = (I_k\ \ \boldsymbol{0})^{\mathsf{T}}$ に使うと $\theta_i^{(k+1)} \ge \theta_i^{(k)} \ge \theta_{i+1}^{(k+1)}$ です。後半はこれと 2 から従います。

4. $AQ_d = Q_dT_d$（<Ref to="thm-arnoldi-relation" /> の 2）なので、$T_d\boldsymbol{s} = \theta\boldsymbol{s}$、$\boldsymbol{s} \ne \boldsymbol{0}$ なら $A(Q_d\boldsymbol{s}) = Q_dT_d\boldsymbol{s} = \theta Q_d\boldsymbol{s}$ で、$Q_d\boldsymbol{s} \ne \boldsymbol{0}$ です。
</Proof>

冪乗法の Rayleigh 商 $R_A(A^{k-1}\boldsymbol{b})$ は $\mathcal{K}_k$ の 1 点での値なので、1 から $\theta_1^{(k)} \ge R_A(A^{k-1}\boldsymbol{b})$ です。同じ $k - 1$ 回の行列ベクトル積で、Lanczos 法はつねに冪乗法以上の近似を与えます。どれだけ上回るかを定量化したのが次の定理で、ここで Chebyshev 多項式 $T_j$（<Ref to="mathematics/computer-assisted-proofs/rigorous-quadrature#def-chebyshev" />）が登場します。使う性質は 2 つ、$[-1, 1]$ で $|T_j(x)| \le 1$ であること（$T_j(\cos\phi) = \cos j\phi$。<Ref to="mathematics/numerical-analysis/polynomial-interpolation#lem-chebyshev-cosine" />）と、$x \ge 1$ で $T_j(x) \ge \tfrac{1}{2}(x + \sqrt{x^2 - 1})^j$ と急に大きくなること（<Ref to="mathematics/numerical-analysis/polynomial-interpolation#lem-chebyshev-outside" />）です。

次の定理の $\gamma$ は、最大固有値 $\lambda_1$ が残りの固有値の区間 $[\lambda_n, \lambda_2]$ からどれだけ離れているかを、区間の幅との比で測った量で、離れているほど Lanczos 法は速く収束します。Chebyshev 多項式を使うのは、区間 $[\lambda_n, \lambda_2]$ の上では絶対値が 1 以下におさまり、その外の $\lambda_1$ では急に大きくなる多項式だからです。証明では、$\mathcal{K}_k$ の元 $p(A)\boldsymbol{b}$ の中で $\boldsymbol{u}_1$ の成分だけが大きく他の成分が小さいものとして、この多項式を選びます。

<Theorem id="thm-kaniel-paige" title="最大の Ritz 値の収束（Kaniel–Paige の評価）">
$\lambda_1 > \lambda_2 > \lambda_n$ とし、$\boldsymbol{u}_1^{*}\boldsymbol{b} \ne 0$、$1 \le k \le d$ とする。
$$
\gamma = \frac{\lambda_1 - \lambda_2}{\lambda_2 - \lambda_n}
$$
とおくと
$$
0 \le \lambda_1 - \theta_1^{(k)} \le (\lambda_1 - \lambda_n)\left(\frac{\tan\angle(\boldsymbol{b}, \boldsymbol{u}_1)}{T_{k-1}(1 + 2\gamma)}\right)^2
$$
が成り立つ。
</Theorem>

<Proof of="thm-kaniel-paige">
左の不等号は <Ref to="thm-ritz-optimality" /> の 2 です。$p \in \mathcal{P}_{k-1}$ で $p(\lambda_1) \ne 0$ なものを 1 つとり、$\boldsymbol{x} = p(A)\boldsymbol{b} \in \mathcal{K}_k$ とします（<Ref to="prop-krylov-polynomial" /> の 1）。$c_i = \boldsymbol{u}_i^{*}\boldsymbol{b}$ とおくと $\boldsymbol{x} = \sum_ic_ip(\lambda_i)\boldsymbol{u}_i \ne \boldsymbol{0}$ で、
$$
\lambda_1 - R_A(\boldsymbol{x}) = \frac{\sum_{i\ge2}|c_i|^2|p(\lambda_i)|^2(\lambda_1 - \lambda_i)}{\sum_{i}|c_i|^2|p(\lambda_i)|^2} \le (\lambda_1 - \lambda_n)\max_{i\ge2}\frac{|p(\lambda_i)|^2}{|p(\lambda_1)|^2}\cdot\frac{\sum_{i\ge2}|c_i|^2}{|c_1|^2}
$$
です（分母を $|c_1|^2|p(\lambda_1)|^2$ に減らした）。最後の比は $\tan^2\angle(\boldsymbol{b}, \boldsymbol{u}_1)$ です（<Ref to="mathematics/numerical-linear-algebra/power-and-inverse-iteration#lem-angle" /> の 4）。<Ref to="thm-ritz-optimality" /> の 1 から $\theta_1^{(k)} \ge R_A(\boldsymbol{x})$ なので、$\lambda_1 - \theta_1^{(k)}$ も同じ式で抑えられます。

ここで $\ell(t) = (2t - \lambda_2 - \lambda_n)/(\lambda_2 - \lambda_n)$ とおき、$p(t) = T_{k-1}(\ell(t))$ とします（$\ell$ は 1 次式なので $p \in \mathcal{P}_{k-1}$）。$\ell$ は $[\lambda_n, \lambda_2]$ を $[-1, 1]$ に写し、$\ell(\lambda_1) = 1 + 2\gamma > 1$ です。$i \ge 2$ では $\ell(\lambda_i) \in [-1, 1]$ なので、$\ell(\lambda_i) = \cos\phi$ と書いて $|p(\lambda_i)| = |\cos((k-1)\phi)| \le 1$（<Ref to="mathematics/numerical-analysis/polynomial-interpolation#lem-chebyshev-cosine" />）です。一方 $p(\lambda_1) = T_{k-1}(1 + 2\gamma) \ge \tfrac12(1 + 2\gamma + \sqrt{(1 + 2\gamma)^2 - 1})^{k-1} > 0$（<Ref to="mathematics/numerical-analysis/polynomial-interpolation#lem-chebyshev-outside" />）です。これを代入すれば主張を得ます。
</Proof>

冪乗法（<Ref to="mathematics/numerical-linear-algebra/power-and-inverse-iteration#cor-power-rayleigh" />。$A$ が正定値なら $|\lambda_2/\lambda_1|$ が比）と比べます。$\gamma$ が小さいとき $T_{k-1}(1 + 2\gamma) \approx \tfrac12 e^{2(k-1)\sqrt{\gamma}}$ なので、Lanczos 法の誤差は 1 回あたり約 $e^{-4\sqrt{\gamma}}$ 倍、冪乗法は $(\lambda_2/\lambda_1)^2 \approx e^{-2(\lambda_1 - \lambda_2)/\lambda_1}$ 倍です。たとえば固有値が $[0, 1]$ に広がり $\lambda_1 - \lambda_2 = 10^{-2}$ なら、$\gamma \approx 10^{-2}$ で Lanczos 法は 1 回あたり $e^{-0.4}$、冪乗法は $e^{-0.02}$ で、Lanczos 法は約 20 倍少ない反復で済みます。差は隔たり $\gamma$ の平方根から生じます。

<Remark id="rem-finite-precision-lanczos" title="有限精度の Lanczos 法：直交性の喪失とゴースト">
三項漸化式は $\boldsymbol{q}_{j+1}$ を直前の 2 本とだけ直交化するので、丸め誤差があると $\boldsymbol{q}_j$ たちの直交性は失われます。Paige は、直交性の喪失がちょうど Ritz 値が固有値に収束したときに起こり、その結果すでに収束した固有値が Ritz 値として何度も現れる（**ゴースト**、複製）ことを示しました。Ritz 値の位置そのものは正しく、失われるのは重複度の情報です。この解析は研究水準の議論で、本章では使いません。Paige の 1971 年の学位論文と 1976・1980 年の論文、Parlett『The Symmetric Eigenvalue Problem』第 13 章、Meurant–Strakoš, "The Lanczos and conjugate gradient algorithms in finite precision arithmetic", *Acta Numerica* 15 (2006) を参照してください。実用では、収束した Ritz ベクトルに対してだけ直交化し直す選択的再直交化などが使われます。
</Remark>

再直交化とは、新しいベクトルを直前の 2 本だけでなく、それまでに作った全部のベクトルと直交させ直す手間のことです。これを省くと丸め誤差で直交性が崩れ、同じ固有値が Ritz 値として何度も現れます（ゴースト）。

次のコードは、固有値が上の端で疎になるように並べた対角行列（Strakoš が有限精度の解析で使った分布）で、Lanczos 法を再直交化なしと完全再直交化で行い、最大固有値 $\lambda_{\max}$ に $10^{-8}$ 以内まで近い Ritz 値の個数を数えたものです。

```python
# Python 3.12.10 / NumPy 2.5.2（x86-64, OpenBLAS）
import numpy as np

n = 100                     # 固有値の分布（Strakoš の例、rho = 0.9）
i = np.arange(n)
lam = 0.1 + i / (n - 1) * 99.9 * 0.9 ** (n - 1 - i)  # A = diag(lam)

def lanczos(k, reorth):
    Q = np.zeros((n, k + 1))
    Q[:, 0] = np.ones(n) / np.sqrt(n)
    al, be = np.zeros(k), np.zeros(k)
    for j in range(k):      # j = 0 では be[-1] = 0 なので第 2 項は 0
        w = lam * Q[:, j] - be[j - 1] * Q[:, j - 1]
        al[j] = Q[:, j] @ w
        w -= al[j] * Q[:, j]
        for _ in range(2 * reorth):   # 完全再直交化（2 回）
            w -= Q[:, :j + 1] @ (Q[:, :j + 1].T @ w)
        be[j] = np.linalg.norm(w)
        Q[:, j + 1] = w / be[j]
    T = np.diag(al) + np.diag(be[:-1], 1) + np.diag(be[:-1], -1)
    return np.linalg.eigvalsh(T), Q[:, :k]

print("k    再直交化  lam_max に 1e-8 以内の Ritz 値の個数")
for k, r in [(k, r) for k in range(20, 101, 20) for r in (0, 1)]:
    th, Q = lanczos(k, r)
    c = np.sum(abs(th - lam[-1]) < 1e-8)
    e = np.linalg.norm(Q.T @ Q - np.eye(k))
    es = "< 1e-14" if e < 1e-14 else f"{e:.0e}"
    print(f"{k:3d}  {'あり' if r else 'なし'}  {c:2d}"
          f"   ||Q^T Q - I|| = {es}")
```

```text
k    再直交化  lam_max に 1e-8 以内の Ritz 値の個数
 20  なし   1   ||Q^T Q - I|| = 3e-08
 20  あり   1   ||Q^T Q - I|| = < 1e-14
 40  なし   2   ||Q^T Q - I|| = 3e+00
 40  あり   1   ||Q^T Q - I|| = < 1e-14
 60  なし   3   ||Q^T Q - I|| = 6e+00
 60  あり   1   ||Q^T Q - I|| = < 1e-14
 80  なし   4   ||Q^T Q - I|| = 9e+00
 80  あり   1   ||Q^T Q - I|| = < 1e-14
100  なし   5   ||Q^T Q - I|| = 1e+01
100  あり   1   ||Q^T Q - I|| = < 1e-14
```

**検算.** (a) 固有値はすべて相異なり（構成した値）、$\lambda_{\max} = 100$ は単純なので、厳密な算術では <Ref to="thm-ritz-optimality" /> の 3 と 2 から $\lambda_{\max}$ の近くの Ritz 値は高々 1 個です。(b) 完全再直交化した計算ではつねに 1 個で、$Q$ の直交性も単位丸めの程度に保たれます。再直交化なしでは、$k = 20$ ですでに直交性が $10^{-8}$ の程度まで崩れ始め、$k = 40$ 以降は約 20 回ごとに $\lambda_{\max}$ の複製が 1 個ずつ増えます。(c) ゴーストが現れる周期は、Ritz 値が収束して直交性が失われる周期に対応しています（<Ref to="rem-finite-precision-lanczos" /> の Paige の理論）。

## 5. 共役勾配法

以下 $A$ は[正定値](/mathematics/linear-algebra/spectral-theorem#def-positive-definite) Hermite 行列、$\boldsymbol{x}^{*} = A^{-1}\boldsymbol{b}$ とし、誤差を $A$ ノルム（エネルギーノルム）$\|\boldsymbol{e}\|_A = (\boldsymbol{e}^{*}A\boldsymbol{e})^{1/2}$ で測ります。$\boldsymbol{x}$ の誤差 $\boldsymbol{x} - \boldsymbol{x}^{*}$ と残差 $\boldsymbol{r} = \boldsymbol{b} - A\boldsymbol{x}$ の間には $\boldsymbol{r} = -A(\boldsymbol{x} - \boldsymbol{x}^{*})$ の関係があります。

<Definition id="def-cg" title="共役勾配法">
$\boldsymbol{x}_0 \in \mathbb{C}^n$ から $\boldsymbol{r}_0 = \boldsymbol{b} - A\boldsymbol{x}_0$、$\boldsymbol{p}_0 = \boldsymbol{r}_0$ とし、$\boldsymbol{r}_k \ne \boldsymbol{0}$ である限り
$$
\begin{aligned}
\alpha_k &= \frac{\boldsymbol{r}_k^{*}\boldsymbol{r}_k}{\boldsymbol{p}_k^{*}A\boldsymbol{p}_k}, & \boldsymbol{x}_{k+1} &= \boldsymbol{x}_k + \alpha_k\boldsymbol{p}_k, & \boldsymbol{r}_{k+1} &= \boldsymbol{r}_k - \alpha_kA\boldsymbol{p}_k,\\
\beta_k &= \frac{\boldsymbol{r}_{k+1}^{*}\boldsymbol{r}_{k+1}}{\boldsymbol{r}_k^{*}\boldsymbol{r}_k}, & \boldsymbol{p}_{k+1} &= \boldsymbol{r}_{k+1} + \beta_k\boldsymbol{p}_k &&
\end{aligned}
$$
を繰り返す。$\boldsymbol{r}_k = \boldsymbol{0}$ になったら止める。これを**共役勾配法**（CG 法）という。
</Definition>

1 回の反復は行列ベクトル積 1 回と内積 2 回・ベクトルの更新 3 回で、保持するベクトルは 4 本です。この単純な算法が、次の強い性質をもちます。

<Theorem id="thm-cg-optimality" title="共役勾配法の直交性と最適性">
$\boldsymbol{r}_0, \ldots, \boldsymbol{r}_{k-1}$ がすべて $\boldsymbol{0}$ でないとする。このとき $\boldsymbol{x}_k$, $\boldsymbol{r}_k$ まで定義され（分母は正）、次が成り立つ。

1. $\boldsymbol{r}_j = \boldsymbol{b} - A\boldsymbol{x}_j$（$j \le k$）。
2. $\boldsymbol{r}_i^{*}\boldsymbol{r}_j = 0$（$0 \le i < j \le k$）、$\boldsymbol{p}_i^{*}A\boldsymbol{p}_j = 0$（$0 \le i < j \le k - 1$）。
3. $\operatorname{span}\{\boldsymbol{r}_0, \ldots, \boldsymbol{r}_{k-1}\} = \operatorname{span}\{\boldsymbol{p}_0, \ldots, \boldsymbol{p}_{k-1}\} = \mathcal{K}_k(A, \boldsymbol{r}_0)$ で、その次元は $k$ である。
4. $\boldsymbol{x}_k \in \boldsymbol{x}_0 + \mathcal{K}_k(A, \boldsymbol{r}_0)$ であり、$\boldsymbol{x}_k$ は $\boldsymbol{x}_0 + \mathcal{K}_k(A, \boldsymbol{r}_0)$ の中で $\|\boldsymbol{x} - \boldsymbol{x}^{*}\|_A$ を最小にするただ 1 つの点である。同じことだが、$\boldsymbol{r}_k \perp \mathcal{K}_k(A, \boldsymbol{r}_0)$ である。
</Theorem>

証明は $k$ についての帰納法で、残差 $\boldsymbol{r}_j$ どうしが直交し、探索方向 $\boldsymbol{p}_j$ どうしが $A$ について直交する、という 2 つの性質を同時に保ちながら進めます。まず算法の分母が正で反復が続くことを確かめ、最後の 4 でこの直交性から $\boldsymbol{x}_k$ が最良点であることを導きます。

<Proof of="thm-cg-optimality">
$k$ についての帰納法で 1〜3 を示します。$k = 0$ では 1 は定義、2・3 は主張がありません。$k$ で成り立ち、$\boldsymbol{r}_k \ne \boldsymbol{0}$ とします。

**分母が正であること.** $k \ge 1$ では 3 から $\boldsymbol{p}_{k-1} \in \operatorname{span}\{\boldsymbol{r}_0, \ldots, \boldsymbol{r}_{k-1}\}$ なので 2 から $\boldsymbol{r}_k^{*}\boldsymbol{p}_{k-1} = 0$ で、$\boldsymbol{p}_k = \boldsymbol{r}_k + \beta_{k-1}\boldsymbol{p}_{k-1}$ から
$$
\boldsymbol{r}_k^{*}\boldsymbol{p}_k = \boldsymbol{r}_k^{*}\boldsymbol{r}_k > 0
$$
です（$k = 0$ でも $\boldsymbol{p}_0 = \boldsymbol{r}_0$ なので同じ）。とくに $\boldsymbol{p}_k \ne \boldsymbol{0}$ で、$A$ は正定値なので $\boldsymbol{p}_k^{*}A\boldsymbol{p}_k > 0$、$\alpha_k > 0$ です。$\beta_k$ の分母も $\boldsymbol{r}_k \ne \boldsymbol{0}$ から正です。$\alpha_k$, $\beta_k$ は正の実数です。

**1.** $\boldsymbol{r}_{k+1} = \boldsymbol{b} - A\boldsymbol{x}_k - \alpha_kA\boldsymbol{p}_k = \boldsymbol{b} - A\boldsymbol{x}_{k+1}$ です。

**$\boldsymbol{p}$ の共役性.** $i < k$ について $\boldsymbol{p}_k^{*}A\boldsymbol{p}_i = 0$ を示します（$k \ge 1$ の場合）。$\alpha_i > 0$ なので $A\boldsymbol{p}_i = (\boldsymbol{r}_i - \boldsymbol{r}_{i+1})/\alpha_i$ です。$i < k - 1$ なら、$\boldsymbol{r}_k^{*}A\boldsymbol{p}_i = 0$（2 で $i, i+1 < k$）と帰納法の仮定 $\boldsymbol{p}_{k-1}^{*}A\boldsymbol{p}_i = 0$ から $\boldsymbol{p}_k^{*}A\boldsymbol{p}_i = 0$ です。$i = k - 1$ なら
$$
\boldsymbol{p}_k^{*}A\boldsymbol{p}_{k-1} = \boldsymbol{r}_k^{*}A\boldsymbol{p}_{k-1} + \beta_{k-1}\boldsymbol{p}_{k-1}^{*}A\boldsymbol{p}_{k-1} = -\frac{\boldsymbol{r}_k^{*}\boldsymbol{r}_k}{\alpha_{k-1}} + \frac{\boldsymbol{r}_k^{*}\boldsymbol{r}_k}{\boldsymbol{r}_{k-1}^{*}\boldsymbol{r}_{k-1}}\,\boldsymbol{p}_{k-1}^{*}A\boldsymbol{p}_{k-1} = 0
$$
です（第 1 項は $\boldsymbol{r}_k^{*}\boldsymbol{r}_{k-1} = 0$ を使い、最後は $\alpha_{k-1}$ の定義）。

**$\boldsymbol{r}$ の直交性.** $i < k$ なら、$\boldsymbol{r}_i \in \operatorname{span}\{\boldsymbol{p}_0, \ldots, \boldsymbol{p}_i\}$（3）と前段落から $\boldsymbol{r}_i^{*}A\boldsymbol{p}_k = 0$ なので、$\boldsymbol{r}_i^{*}\boldsymbol{r}_{k+1} = \boldsymbol{r}_i^{*}\boldsymbol{r}_k - \alpha_k\boldsymbol{r}_i^{*}A\boldsymbol{p}_k = 0$ です。$i = k$ では、$\boldsymbol{r}_k = \boldsymbol{p}_k - \beta_{k-1}\boldsymbol{p}_{k-1}$ と共役性から $\boldsymbol{r}_k^{*}A\boldsymbol{p}_k = \boldsymbol{p}_k^{*}A\boldsymbol{p}_k$ なので（$k = 0$ では $\boldsymbol{r}_0 = \boldsymbol{p}_0$）、$\boldsymbol{r}_k^{*}\boldsymbol{r}_{k+1} = \boldsymbol{r}_k^{*}\boldsymbol{r}_k - \alpha_k\boldsymbol{p}_k^{*}A\boldsymbol{p}_k = 0$ です。これで 2 が $k + 1$ で成り立ちます。

**3.** 帰納法で、$\boldsymbol{r}_j, \boldsymbol{p}_j \in \mathcal{K}_{j+1}(A, \boldsymbol{r}_0)$ です（$\boldsymbol{r}_{j+1} = \boldsymbol{r}_j - \alpha_jA\boldsymbol{p}_j$ で $A\mathcal{K}_{j+1} \subset \mathcal{K}_{j+2}$、$\boldsymbol{p}_{j+1} = \boldsymbol{r}_{j+1} + \beta_j\boldsymbol{p}_j$）。$\boldsymbol{r}_0, \ldots, \boldsymbol{r}_k$ は零でなく互いに直交するので一次独立で、$k + 1$ 次元の空間を張ります。それは高々 $k + 1$ 次元の $\mathcal{K}_{k+1}$ に含まれるので一致します。$\boldsymbol{p}$ の張る空間は、$\boldsymbol{p}_j = \boldsymbol{r}_j + \beta_{j-1}\boldsymbol{p}_{j-1}$ から帰納的に $\operatorname{span}\{\boldsymbol{r}_0, \ldots, \boldsymbol{r}_j\} = \operatorname{span}\{\boldsymbol{p}_0, \ldots, \boldsymbol{p}_j\}$ です。

**4.** $\boldsymbol{x}_k = \boldsymbol{x}_0 + \sum_{j<k}\alpha_j\boldsymbol{p}_j \in \boldsymbol{x}_0 + \mathcal{K}_k$ で、2 と 3 から $\boldsymbol{r}_k \perp \mathcal{K}_k$ です。$\boldsymbol{x} = \boldsymbol{x}_k + \boldsymbol{z}$（$\boldsymbol{z} \in \mathcal{K}_k$）とすると、$\boldsymbol{e}_k = \boldsymbol{x}_k - \boldsymbol{x}^{*}$、$A\boldsymbol{e}_k = -\boldsymbol{r}_k$ から
$$
\|\boldsymbol{x} - \boldsymbol{x}^{*}\|_A^2 = \|\boldsymbol{e}_k\|_A^2 + 2\operatorname{Re}(\boldsymbol{z}^{*}A\boldsymbol{e}_k) + \|\boldsymbol{z}\|_A^2 = \|\boldsymbol{e}_k\|_A^2 - 2\operatorname{Re}(\boldsymbol{z}^{*}\boldsymbol{r}_k) + \|\boldsymbol{z}\|_A^2 = \|\boldsymbol{e}_k\|_A^2 + \|\boldsymbol{z}\|_A^2
$$
で、最小は $\boldsymbol{z} = \boldsymbol{0}$ でだけとられます。逆に $\boldsymbol{x}_0 + \mathcal{K}_k$ の点 $\boldsymbol{x}$ の残差が $\mathcal{K}_k$ に直交するなら、同じ計算を $\boldsymbol{x}$ を基準に行えば $\boldsymbol{x}$ が最小点で、一意性から $\boldsymbol{x} = \boldsymbol{x}_k$ です。
</Proof>

名前の「共役」は 2 の後半（探索方向が $A$ に関して直交すること）、「勾配」は $\boldsymbol{r} = -\nabla\bigl(\tfrac12\boldsymbol{x}^{*}A\boldsymbol{x} - \operatorname{Re}\boldsymbol{b}^{*}\boldsymbol{x}\bigr)$ が汎関数の最急降下の方向であることに由来します。

<Corollary id="cor-cg-finite-termination" title="有限回での停止">
厳密な算術では、共役勾配法は $\boldsymbol{r}_0$ の次数 $d$ 以下の回数で $\boldsymbol{r}_k = \boldsymbol{0}$（$\boldsymbol{x}_k = \boldsymbol{x}^{*}$）となって止まる。とくに反復回数は $A$ の相異なる固有値の個数以下である。
</Corollary>

<Proof of="cor-cg-finite-termination">
$\boldsymbol{r}_0, \ldots, \boldsymbol{r}_{k-1}$ がすべて零でなければ、<Ref to="thm-cg-optimality" /> の 3 から $\dim\mathcal{K}_k(A, \boldsymbol{r}_0) = k$ なので、<Ref to="prop-krylov-polynomial" /> の 2 から $k \le d$ です。よって $\boldsymbol{r}_0, \ldots, \boldsymbol{r}_{d-1}$ がすべて零でないとしても、$\boldsymbol{r}_d \in \mathcal{K}_{d+1} = \mathcal{K}_d$（<Ref to="prop-krylov-polynomial" /> の 2）かつ $\boldsymbol{r}_d \perp \mathcal{K}_d$（<Ref to="thm-cg-optimality" /> の 4）から $\boldsymbol{r}_d = \boldsymbol{0}$ です。Hermite 行列は対角化可能なので（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-spectral" />）、<Ref to="prop-krylov-polynomial" /> の 3 から $d$ は相異なる固有値の個数以下です。
</Proof>

### 5.1. Chebyshev 多項式による収束の評価

共役勾配法は $\boldsymbol{x}_0 + \mathcal{K}_k$ の中の最良点を選ぶので、その誤差は、$p(0) = 1$ の $k$ 次多項式 $p$ を使った $p(A)\boldsymbol{e}_0$ の最小値になります。こう言い換えると、固有値の上で小さな値をとる多項式を 1 つ見つけるだけで、誤差の上界が得られます。最適性を多項式の言葉に翻訳します。$\boldsymbol{x} \in \boldsymbol{x}_0 + \mathcal{K}_k(A, \boldsymbol{r}_0)$ は $\boldsymbol{x} = \boldsymbol{x}_0 + q(A)\boldsymbol{r}_0$（$q \in \mathcal{P}_{k-1}$）と書け、$\boldsymbol{r}_0 = -A\boldsymbol{e}_0$ なので
$$
\boldsymbol{x} - \boldsymbol{x}^{*} = \boldsymbol{e}_0 - q(A)A\boldsymbol{e}_0 = p(A)\boldsymbol{e}_0,\qquad p(t) = 1 - tq(t)
$$
です。$q$ が $\mathcal{P}_{k-1}$ を動くとき、$p$ は $p(0) = 1$ を満たす $\mathcal{P}_k$ の元全体を動きます。

<Theorem id="thm-cg-chebyshev-bound" title="共役勾配法の収束の評価">
$A$ の最大・最小の固有値を $\lambda_{\max}$, $\lambda_{\min}$、$\kappa = \lambda_{\max}/\lambda_{\min}$ とする。共役勾配法の誤差 $\boldsymbol{e}_k = \boldsymbol{x}_k - \boldsymbol{x}^{*}$ は（停止した後は $\boldsymbol{e}_k = \boldsymbol{0}$ とみなして）
$$
\|\boldsymbol{e}_k\|_A = \min_{p\in\mathcal{P}_k,\ p(0)=1}\|p(A)\boldsymbol{e}_0\|_A \le \min_{p\in\mathcal{P}_k,\ p(0)=1}\ \max_{\lambda\in\sigma(A)}|p(\lambda)|\ \|\boldsymbol{e}_0\|_A \le 2\left(\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}\right)^{k}\|\boldsymbol{e}_0\|_A
$$
を満たす。ここで $\sigma(A)$ は $A$ の固有値の集合である。
</Theorem>

<Proof of="thm-cg-chebyshev-bound">
停止した後は $\boldsymbol{e}_k = \boldsymbol{0}$ なので、停止していない $k$ を考えます。等号は <Ref to="thm-cg-optimality" /> の 4 と上の翻訳から従います。$\boldsymbol{e}_0 = \sum_ic_i\boldsymbol{u}_i$ と正規直交な固有ベクトルで展開すると（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-spectral" />）
$$
\|p(A)\boldsymbol{e}_0\|_A^2 = \sum_i\lambda_i|p(\lambda_i)|^2|c_i|^2 \le \max_i|p(\lambda_i)|^2\sum_i\lambda_i|c_i|^2 = \max_i|p(\lambda_i)|^2\,\|\boldsymbol{e}_0\|_A^2
$$
で、2 つ目の不等号を得ます。

$\kappa = 1$ なら $A = \lambda_{\max}I$ で、$p(t) = 1 - t/\lambda_{\max}$ が $\sigma(A)$ 上で $0$ になるので主張は明らかです。$\kappa > 1$ とし、
$$
p(t) = \frac{T_k(\ell(t))}{T_k(\ell(0))},\qquad \ell(t) = \frac{\lambda_{\max} + \lambda_{\min} - 2t}{\lambda_{\max} - \lambda_{\min}}
$$
とおきます。$\ell$ は $[\lambda_{\min}, \lambda_{\max}]$ を $[-1, 1]$ に写し、$\ell(0) = (\kappa + 1)/(\kappa - 1) > 1$ です。$T_k(\ell(0)) > 0$（<Ref to="mathematics/numerical-analysis/polynomial-interpolation#lem-chebyshev-outside" />）なので $p$ は定義でき、$p \in \mathcal{P}_k$、$p(0) = 1$ です。固有値 $\lambda$ では $\ell(\lambda) = \cos\phi$ と書けて $|T_k(\ell(\lambda))| = |\cos k\phi| \le 1$（<Ref to="mathematics/numerical-analysis/polynomial-interpolation#lem-chebyshev-cosine" />）です。分母は
$$
\ell(0) + \sqrt{\ell(0)^2 - 1} = \frac{\kappa + 1 + 2\sqrt{\kappa}}{\kappa - 1} = \frac{(\sqrt{\kappa} + 1)^2}{(\sqrt{\kappa} - 1)(\sqrt{\kappa} + 1)} = \frac{\sqrt{\kappa} + 1}{\sqrt{\kappa} - 1}
$$
と <Ref to="mathematics/numerical-analysis/polynomial-interpolation#lem-chebyshev-outside" /> から $T_k(\ell(0)) \ge \tfrac12\bigl((\sqrt{\kappa} + 1)/(\sqrt{\kappa} - 1)\bigr)^k$ です。よって $\max_{\lambda\in\sigma(A)}|p(\lambda)| \le 2\bigl((\sqrt{\kappa} - 1)/(\sqrt{\kappa} + 1)\bigr)^k$ です。
</Proof>

<Figure caption="κ = 10、k = 5 のときの評価に使う多項式 p(t) = T₅(ℓ(t))/T₅(ℓ(0))。p(0) = 1 を保ったまま、固有値の区間 [1, 10] の上では ±1/T₅(11/9) ≈ ±0.0756 の間（網かけの帯）を等しい振幅で振動する。この振幅は評価式の 2((√10 − 1)/(√10 + 1))⁵ ≈ 0.0757 とほぼ等しい。">
<svg viewBox="0 0 720 360" width="100%" role="img" aria-label="区間 0 から 10 の上の 5 次の Chebyshev 型多項式のグラフ。t = 0 で 1、区間 1 から 10 の上では小さな振幅で 5 回振動する。">
  <rect x="130" y="261" width="540" height="38" fill="var(--sl-color-accent)" opacity="0.12" />
  <g stroke="currentColor" stroke-width="1" opacity="0.4">
    <line x1="70" y1="280" x2="680" y2="280" />
    <line x1="70" y1="20" x2="70" y2="335" />
    <line x1="130" y1="275" x2="130" y2="285" />
    <line x1="670" y1="275" x2="670" y2="285" />
  </g>
  <g fill="currentColor" font-size="13" opacity="0.8">
    <text x="62" y="34" text-anchor="end">1</text>
    <text x="62" y="284" text-anchor="end">0</text>
    <text x="70" y="352" text-anchor="middle">0</text>
    <text x="130" y="352" text-anchor="middle">1 = λ_min</text>
    <text x="400" y="352" text-anchor="middle">5.5</text>
    <text x="670" y="352" text-anchor="middle">10 = λ_max</text>
    <text x="690" y="284">t</text>
  </g>
  <polyline fill="none" stroke="var(--sl-color-accent)" stroke-width="2.5" stroke-linejoin="round" points="70,30 85,116 100,181 115,229 130,261 145,282 160,294 175,298 190,298 205,295 220,289 235,283 250,276 265,271 280,266 295,263 310,261 325,261 340,263 355,266 370,270 385,275 400,280 415,285 430,290 445,294 460,297 475,299 490,299 505,297 520,294 535,289 550,284 565,277 580,271 595,265 610,262 625,262 640,266 655,278 670,299" />
  <circle cx="70" cy="30" r="4" fill="var(--sl-color-accent)" />
  <text x="82" y="30" font-size="13" fill="currentColor">p(0) = 1</text>
</svg>
</Figure>

$(\sqrt{\kappa} - 1)/(\sqrt{\kappa} + 1) \approx 1 - 2/\sqrt{\kappa}$ なので、誤差を $\varepsilon$ 倍にする反復回数は約 $\tfrac12\sqrt{\kappa}\ln(2/\varepsilon)$ です。2 次元の離散 Laplace 行列では $\kappa \approx 4/(\pi^2h^2)$ なので回数は $h^{-1}$ に比例し、最適 SOR 法と同じ次数を、パラメータの推定なしで実現します。評価は固有値の**区間**だけを使うので、固有値がいくつかの塊に集まっていれば、実際の収束はずっと速くなります（<Ref to="ex-clustered-eigenvalues" />）。

<Proposition id="prop-steepest-descent" title="最急降下法の収束">
最急降下法 $\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \alpha_k\boldsymbol{r}_k$、$\alpha_k = \boldsymbol{r}_k^{*}\boldsymbol{r}_k/\boldsymbol{r}_k^{*}A\boldsymbol{r}_k$（$\boldsymbol{r}_k = \boldsymbol{b} - A\boldsymbol{x}_k \ne \boldsymbol{0}$）の誤差は
$$
\|\boldsymbol{e}_{k+1}\|_A \le \frac{\kappa - 1}{\kappa + 1}\|\boldsymbol{e}_k\|_A
$$
を満たす。
</Proposition>

<Proof of="prop-steepest-descent">
$\boldsymbol{r} = \boldsymbol{r}_k$、$\boldsymbol{e} = \boldsymbol{e}_k$ と書くと $A\boldsymbol{e} = -\boldsymbol{r}$、$\boldsymbol{e}_{k+1} = \boldsymbol{e} + \alpha_k\boldsymbol{r}$ で、
$$
\|\boldsymbol{e}_{k+1}\|_A^2 = \boldsymbol{e}^{*}A\boldsymbol{e} + 2\alpha_k\operatorname{Re}(\boldsymbol{r}^{*}A\boldsymbol{e}) + \alpha_k^2\boldsymbol{r}^{*}A\boldsymbol{r} = \|\boldsymbol{e}\|_A^2 - 2\alpha_k\boldsymbol{r}^{*}\boldsymbol{r} + \alpha_k^2\boldsymbol{r}^{*}A\boldsymbol{r} = \|\boldsymbol{e}\|_A^2 - \frac{(\boldsymbol{r}^{*}\boldsymbol{r})^2}{\boldsymbol{r}^{*}A\boldsymbol{r}}
$$
です。$\|\boldsymbol{e}\|_A^2 = \boldsymbol{e}^{*}A\boldsymbol{e} = \boldsymbol{r}^{*}A^{-1}\boldsymbol{r}$ なので
$$
\frac{\|\boldsymbol{e}_{k+1}\|_A^2}{\|\boldsymbol{e}\|_A^2} = 1 - \frac{(\boldsymbol{r}^{*}\boldsymbol{r})^2}{(\boldsymbol{r}^{*}A\boldsymbol{r})(\boldsymbol{r}^{*}A^{-1}\boldsymbol{r})} \le 1 - \frac{4\lambda_{\max}\lambda_{\min}}{(\lambda_{\max} + \lambda_{\min})^2} = \left(\frac{\lambda_{\max} - \lambda_{\min}}{\lambda_{\max} + \lambda_{\min}}\right)^2
$$
です。不等号は Kantorovich の不等式（<Ref to="mathematics/linear-algebra/matrix-inequalities#thm-kantorovich" />）を $m = \lambda_{\min}$, $M = \lambda_{\max}$ で使いました。
</Proof>

最急降下法は 1 回ごとに最適な歩幅をとりますが、過去の方向を忘れるので、反復回数は $\kappa$ に比例します。共役勾配法は過去の方向すべてと $A$ 共役な方向をとることで、$\mathcal{K}_k$ 全体での最適化を 1 本の漸化式で実現し、回数を $\sqrt{\kappa}$ に下げています。

<Example id="ex-clustered-eigenvalues" title="固有値が塊に集まると速く終わる">
固有値が $1$ の近くの $s$ 個の値と、$[100, 101]$ の中の値だけからなるとします。<Ref to="thm-cg-chebyshev-bound" /> の最後の評価は $\kappa \approx 101$ で、1 回あたり $(\sqrt{101} - 1)/(\sqrt{101} + 1) \approx 0.82$ 倍しか保証しません。しかし真ん中の評価で $p(t) = \prod_{i=1}^{s}(1 - t/\mu_i)\cdot\tilde{p}(t)$（$\mu_i$ は小さい固有値、$\tilde{p}$ は $[100, 101]$ に合わせた Chebyshev 多項式で $\tilde{p}(0) = 1$）を選ぶと、$p$ は小さい固有値で $0$、$[100, 101]$ の上では $|\tilde{p}| \cdot\prod_i|1 - t/\mu_i|$ 以下です。$[100,101]$ の条件数は $1.01$ なので $\tilde{p}$ は 1 回あたり約 $(\sqrt{1.01} - 1)/(\sqrt{1.01} + 1) \approx 0.0025$ 倍で減り、$s$ 回の「準備」のあとは非常に速く収束します（$\mu_i \approx 1$ なら $[100, 101]$ の上で $\prod_i|1 - t/\mu_i|$ は約 $100^s$ なので、その分を取り返す数回の遅れは生じます）。前処理の目的は、条件数を下げることに加えて、固有値をこのような塊に集めることです。
</Example>

次のコードは、2 次元の離散 Laplace 方程式（右辺は $f = 1$、$\boldsymbol{x}_0 = \boldsymbol{0}$）を `scipy.sparse.linalg.cg` で相対残差 $10^{-8}$ まで解いた反復回数です。$\sqrt{\kappa}$ は閉じた形の固有値 $8\sin^2(\pi h/2)$ と $8\cos^2(\pi h/2)$ から計算しました。「予測」は <Ref to="thm-cg-chebyshev-bound" /> から相対残差 $10^{-8}$ を保証する回数 $\ln(2\sqrt{\kappa}\cdot10^{8})/\ln\bigl((\sqrt{\kappa}+1)/(\sqrt{\kappa}-1)\bigr)$ です（$\|\boldsymbol{r}\| \le \sqrt{\lambda_{\max}}\|\boldsymbol{e}\|_A$、$\|\boldsymbol{e}\|_A \le \|\boldsymbol{r}\|/\sqrt{\lambda_{\min}}$ で残差に換算した分 $\sqrt{\kappa}$ がかかります）。

```python
# Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1（x86-64, OpenBLAS）
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import cg

def laplace2d(m):
    T = sp.diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(m, m),
                 dtype=float)
    I = sp.identity(m, dtype=float)
    return (sp.kron(I, T) + sp.kron(T, I)).tocsr()

print(" 1/h   反復   sqrt(kappa)  反復/sqrt(kappa)  予測")
for m in (31, 63, 127, 255):
    h = 1 / (m + 1)
    A = laplace2d(m)
    b = np.ones(m * m) * h**2
    it = [0]
    def cb(xk):
        it[0] += 1
    x, info = cg(A, b, rtol=1e-8, callback=cb)
    lmin = 8 * np.sin(np.pi * h / 2) ** 2      # 閉じた形
    lmax = 8 * np.cos(np.pi * h / 2) ** 2
    sk = np.sqrt(lmax / lmin)
    q = (sk - 1) / (sk + 1)
    pred = np.log(2e8 * sk) / -np.log(q)
    res = np.linalg.norm(b - A @ x) / np.linalg.norm(b)
    assert info == 0 and res <= 1e-8
    print(f"{m+1:4d}  {it[0]:5d}  {sk:9.2f}  {it[0]/sk:12.3f}"
          f"  {pred:8.0f}")
```

```text
 1/h   反復   sqrt(kappa)  反復/sqrt(kappa)  予測
  32     58      20.36         2.849       225
  64    118      40.74         2.897       465
 128    237      81.48         2.909       958
 256    468     162.97         2.872      1973
```

**検算.** (a) $\sqrt{\kappa}$ は閉じた形の固有値（[ノルムの章](/mathematics/numerical-linear-algebra/norms-and-conditioning) の三重対角 Toeplitz 行列の固有値の和）から計算し、$\sqrt{\kappa} \approx 2/(\pi h)$（$h = 1/256$ で $162.97$ 対 $162.97$）とも合います。(b) 各行で `assert` により、返された解の相対残差が実際に $10^{-8}$ 以下であることを確かめています。(c) 格子幅を半分にするごとに反復回数は $2.03, 2.01, 1.97$ 倍で、反復回数 $/\sqrt{\kappa}$ はほぼ一定（約 $2.9$）です。<Ref to="thm-cg-chebyshev-bound" /> の予言どおり $\sqrt{\kappa} \propto h^{-1}$ に比例します。実際の回数は評価から保証される回数の 4 分の 1 程度で、評価は上界として成り立っています。[定常反復法の章](/mathematics/numerical-linear-algebra/stationary-iterative-methods) の最適 SOR 法（$h = 1/64$ で 257 回）と比べても半分以下です。

## 6. 共役勾配法と Lanczos 法、前処理

<Theorem id="thm-cg-lanczos" title="共役勾配法と Lanczos 法の関係">
共役勾配法が $k$ 回目まで停止しないとし、$\boldsymbol{b}$ の代わりに $\boldsymbol{r}_0$ から始めた Lanczos 法（<Ref to="cor-lanczos-three-term" />）の $Q_k$, $T_k$ を考える。このとき

1. 各 $j < k$ について、$\boldsymbol{r}_j$ は $\boldsymbol{q}_{j+1}$ の定数倍で、$\|\boldsymbol{r}_j\| = |\boldsymbol{q}_{j+1}^{*}\boldsymbol{r}_j|$ である。
2. $T_k$ は正定値で、$\boldsymbol{x}_k = \boldsymbol{x}_0 + Q_k\boldsymbol{y}_k$、$T_k\boldsymbol{y}_k = \|\boldsymbol{r}_0\|\boldsymbol{e}_1$ である（$\boldsymbol{e}_1$ は第 1 単位ベクトル）。
</Theorem>

<Proof of="thm-cg-lanczos">
1. <Ref to="thm-cg-optimality" /> の 3 から $\mathcal{K}_k(A, \boldsymbol{r}_0)$ の次元は $k$ なので、$\boldsymbol{r}_0$ の次数は $k$ 以上で、Lanczos 法は $k$ 本の $\boldsymbol{q}_j$ を作ります（<Ref to="thm-arnoldi-relation" />）。$\boldsymbol{r}_j \in \mathcal{K}_{j+1}$ で、$\boldsymbol{r}_j \perp \mathcal{K}_j$（<Ref to="thm-cg-optimality" /> の 2・3）です。$\mathcal{K}_{j+1}$ の中で $\mathcal{K}_j = \operatorname{span}\{\boldsymbol{q}_1, \ldots, \boldsymbol{q}_j\}$ に直交するベクトルは $\boldsymbol{q}_{j+1}$ の定数倍だけなので（$\boldsymbol{q}_1, \ldots, \boldsymbol{q}_{j+1}$ は $\mathcal{K}_{j+1}$ の正規直交基底）、$\boldsymbol{r}_j = c\,\boldsymbol{q}_{j+1}$ で、$|c| = \|\boldsymbol{r}_j\|$ です。

2. $T_k = Q_k^{*}AQ_k$（<Ref to="thm-arnoldi-relation" /> の 3）で、$\boldsymbol{y} \ne \boldsymbol{0}$ なら $\boldsymbol{y}^{*}T_k\boldsymbol{y} = (Q_k\boldsymbol{y})^{*}A(Q_k\boldsymbol{y}) > 0$ なので正定値、とくに正則です。$\boldsymbol{x}_k \in \boldsymbol{x}_0 + \mathcal{K}_k$ なので $\boldsymbol{x}_k = \boldsymbol{x}_0 + Q_k\boldsymbol{y}_k$ と書け、$\boldsymbol{r}_k = \boldsymbol{r}_0 - AQ_k\boldsymbol{y}_k \perp \mathcal{K}_k$（<Ref to="thm-cg-optimality" /> の 4）から $Q_k^{*}\boldsymbol{r}_0 - T_k\boldsymbol{y}_k = \boldsymbol{0}$ です。$\boldsymbol{q}_1 = \boldsymbol{r}_0/\|\boldsymbol{r}_0\|$ と直交性から $Q_k^{*}\boldsymbol{r}_0 = \|\boldsymbol{r}_0\|\boldsymbol{e}_1$ です。
</Proof>

したがって、共役勾配法を走らせながら $\alpha_k$, $\beta_k$ から $T_k$ を組み立てれば、Ritz 値（$A$ の端の固有値の近似）が副産物として得られます。条件数の推定や停止判定に使われる技法です。逆に、Lanczos 法で作った $T_k$ の小さな系 $T_k\boldsymbol{y} = \|\boldsymbol{r}_0\|\boldsymbol{e}_1$ を解けば共役勾配法の近似解が得られます。

前処理は $M^{-1}A$ の条件数を小さくするための変形ですが、$M^{-1}A$ はふつう Hermite にならず、そのままでは共役勾配法を使いにくくなります。そこで $M$ の平方根 $M^{1/2}$ で両側から挟んだ $M^{-1/2}AM^{-1/2}$ を考えます。これは正定値 Hermite のままで、$M^{-1}A$ と相似なので固有値も同じです。次の定理は、この変形を実際には行わず $M^{-1}$ を掛けるだけの反復が、変形した行列への共役勾配法と一致することを示します。

<Theorem id="thm-pcg-equivalence" title="前処理つき共役勾配法">
$M$ を正定値 Hermite 行列とし、$M^{1/2}$ をその正定値な平方根とする。**前処理つき共役勾配法**を、$\boldsymbol{r}_0 = \boldsymbol{b} - A\boldsymbol{x}_0$、$\boldsymbol{z}_0 = M^{-1}\boldsymbol{r}_0$、$\boldsymbol{p}_0 = \boldsymbol{z}_0$ とし、
$$
\begin{aligned}
\alpha_k &= \frac{\boldsymbol{r}_k^{*}\boldsymbol{z}_k}{\boldsymbol{p}_k^{*}A\boldsymbol{p}_k}, & \boldsymbol{x}_{k+1} &= \boldsymbol{x}_k + \alpha_k\boldsymbol{p}_k, & \boldsymbol{r}_{k+1} &= \boldsymbol{r}_k - \alpha_kA\boldsymbol{p}_k,\\
\boldsymbol{z}_{k+1} &= M^{-1}\boldsymbol{r}_{k+1}, & \beta_k &= \frac{\boldsymbol{r}_{k+1}^{*}\boldsymbol{z}_{k+1}}{\boldsymbol{r}_k^{*}\boldsymbol{z}_k}, & \boldsymbol{p}_{k+1} &= \boldsymbol{z}_{k+1} + \beta_k\boldsymbol{p}_k
\end{aligned}
$$
で定める。$\hat{A} = M^{-1/2}AM^{-1/2}$、$\hat{\boldsymbol{b}} = M^{-1/2}\boldsymbol{b}$ とおくと、$\hat{\boldsymbol{x}}_k = M^{1/2}\boldsymbol{x}_k$ は $\hat{A}\hat{\boldsymbol{x}} = \hat{\boldsymbol{b}}$ に $\hat{\boldsymbol{x}}_0 = M^{1/2}\boldsymbol{x}_0$ から共役勾配法を使った反復と一致する。とくに、$M^{-1}A$ の固有値（正の実数）の最大と最小の比を $\hat{\kappa}$ として
$$
\|\boldsymbol{x}_k - \boldsymbol{x}^{*}\|_A \le 2\left(\frac{\sqrt{\hat{\kappa}} - 1}{\sqrt{\hat{\kappa}} + 1}\right)^k\|\boldsymbol{x}_0 - \boldsymbol{x}^{*}\|_A .
$$
</Theorem>

<Proof of="thm-pcg-equivalence">
**平方根.** スペクトル定理（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-spectral" />）で $M = U\Lambda U^{*}$ と書くと、正定値性から $\Lambda$ の対角成分は正です（<Ref to="mathematics/linear-algebra/spectral-theorem#thm-principal-axis" /> の 2）。$M^{1/2} = U\Lambda^{1/2}U^{*}$ は正定値 Hermite で $(M^{1/2})^2 = M$、逆行列は $M^{-1/2} = U\Lambda^{-1/2}U^{*}$ です。$\hat{A}$ は Hermite で、$\boldsymbol{y} \ne \boldsymbol{0}$ なら $\boldsymbol{y}^{*}\hat{A}\boldsymbol{y} = (M^{-1/2}\boldsymbol{y})^{*}A(M^{-1/2}\boldsymbol{y}) > 0$ なので正定値です。

**反復の対応.** $\hat{\boldsymbol{r}}_k = M^{-1/2}\boldsymbol{r}_k$、$\hat{\boldsymbol{p}}_k = M^{1/2}\boldsymbol{p}_k$ とおき、これらが $\hat{A}$ の共役勾配法の量であることを $k$ についての帰納法で示します。$\hat{\boldsymbol{r}}_0 = M^{-1/2}(\boldsymbol{b} - A\boldsymbol{x}_0) = \hat{\boldsymbol{b}} - \hat{A}\hat{\boldsymbol{x}}_0$、$\hat{\boldsymbol{p}}_0 = M^{1/2}M^{-1}\boldsymbol{r}_0 = \hat{\boldsymbol{r}}_0$ です。各段で
$$
\hat{\boldsymbol{r}}_k^{*}\hat{\boldsymbol{r}}_k = \boldsymbol{r}_k^{*}M^{-1}\boldsymbol{r}_k = \boldsymbol{r}_k^{*}\boldsymbol{z}_k,\qquad \hat{\boldsymbol{p}}_k^{*}\hat{A}\hat{\boldsymbol{p}}_k = \boldsymbol{p}_k^{*}A\boldsymbol{p}_k
$$
なので $\alpha_k$, $\beta_k$ は両方の算法で等しく、$\hat{\boldsymbol{x}}_{k+1} = \hat{\boldsymbol{x}}_k + \alpha_k\hat{\boldsymbol{p}}_k$、$\hat{\boldsymbol{r}}_{k+1} = M^{-1/2}(\boldsymbol{r}_k - \alpha_kA\boldsymbol{p}_k) = \hat{\boldsymbol{r}}_k - \alpha_k\hat{A}\hat{\boldsymbol{p}}_k$、$\hat{\boldsymbol{p}}_{k+1} = M^{1/2}(M^{-1}\boldsymbol{r}_{k+1} + \beta_k\boldsymbol{p}_k) = \hat{\boldsymbol{r}}_{k+1} + \beta_k\hat{\boldsymbol{p}}_k$ です。これは $\hat{A}$ に対する <Ref to="def-cg" /> そのものです。

**評価.** $\hat{\boldsymbol{x}}^{*} = M^{1/2}\boldsymbol{x}^{*}$ は $\hat{A}\hat{\boldsymbol{x}}^{*} = \hat{\boldsymbol{b}}$ の解で、$\hat{\boldsymbol{e}} = M^{1/2}\boldsymbol{e}$ について $\hat{\boldsymbol{e}}^{*}\hat{A}\hat{\boldsymbol{e}} = \boldsymbol{e}^{*}A\boldsymbol{e}$ です。$M^{-1}A = M^{-1/2}\hat{A}M^{1/2}$ は $\hat{A}$ と相似なので固有値が等しく、それは $\hat{A}$ の固有値（正の実数）です。<Ref to="thm-cg-chebyshev-bound" /> を $\hat{A}$ に使えば主張を得ます。
</Proof>

実装では $M^{1/2}$ を作る必要はなく、毎回 $M\boldsymbol{z} = \boldsymbol{r}$ を解くだけです。$M$ は「$A$ に近く、$M\boldsymbol{z} = \boldsymbol{r}$ が安く解ける」ように選びます。[定常反復法の章](/mathematics/numerical-linear-algebra/stationary-iterative-methods) の SSOR 法の分割行列や、不完全 Cholesky 分解がその代表で、次の Krylov 部分空間法 II の章で扱います。

<Remark id="rem-cg-finite-precision" title="有限精度での共役勾配法：収束の遅れと超 1 次収束">
<Ref to="cor-cg-finite-termination" /> の有限回での停止は厳密な算術での話で、実際の計算では残差の直交性が失われ、$n$ 回で止まるとは限りません。それでも共役勾配法は実用上よく収束します。Greenbaum は、有限精度の共役勾配法の誤差が、元の行列の各固有値の近くに固有値の小さな塊をもつ大きな行列に対する厳密な共役勾配法の誤差で説明できることを示しました（収束は遅れるが、失敗はしない）。また、Ritz 値が端の固有値に収束すると、その固有値は以後「取り除かれた」ように振る舞い、実効的な条件数が下がって収束が加速します（超 1 次収束。van der Sluis–van der Vorst）。どちらも本章の定理には使いません。Greenbaum『Iterative Methods for Solving Linear Systems』SIAM, 1997 — 第 4 章、Meurant–Strakoš（*Acta Numerica* 15, 2006）、A. van der Sluis, H. A. van der Vorst, "The rate of convergence of conjugate gradients", *Numerische Mathematik* 48 (1986) を参照してください。
</Remark>

## 7. 疎行列の固有値：シフト反転と ARPACK

研究でよく使う `scipy.sparse.linalg.eigsh` は、ARPACK の**陰的再開 Lanczos 法**を呼びます。Lanczos 法は端の固有値ほど速く収束する（<Ref to="thm-kaniel-paige" />）ので、中ほどや最小の固有値が欲しいときは**シフト反転**、すなわち $(A - \sigma I)^{-1}$ に Lanczos 法を使います。$\sigma$ に近い固有値 $\lambda$ は $(A - \sigma I)^{-1}$ では絶対値の大きい固有値 $1/(\lambda - \sigma)$ になり、しかも互いによく分離されます（[冪乗法の章](/mathematics/numerical-linear-algebra/power-and-inverse-iteration) の逆反復と同じ考え方）。$(A - \sigma I)^{-1}$ との積は、疎な LU 分解を 1 回作っておいて毎回代入で計算します。

```python
# Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1（x86-64, OpenBLAS）／乱数種 20261006
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import eigsh

m = 100                                   # 100×100 格子、h = 1/101
h = 1 / (m + 1)
T = sp.diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(m, m), dtype=float)
I = sp.identity(m, dtype=float)
A = (sp.kron(I, T) + sp.kron(T, I)).tocsc()

rng = np.random.default_rng(20261006)
v0 = rng.standard_normal(m * m)            # 開始ベクトルを固定
vals = eigsh(A, k=5, sigma=0, which="LM", v0=v0,
             return_eigenvectors=False)
vals = np.sort(vals)

t = 4 * np.sin(np.arange(1, 4) * np.pi * h / 2) ** 2  # T の固有値
exact = np.sort([t[j] + t[k] for j in range(3) for k in range(3)])[:5]
print(" 固有値（eigsh）     閉じた形")
for a, e in zip(vals, exact):
    print(f"{a:.10e}  {e:.10e}")
err = max(abs(vals - exact) / exact)
print("相対誤差の最大 < 1e-12:", err < 1e-12)
```

```text
 固有値（eigsh）     閉じた形
1.9348708320e-03  1.9348708320e-03
4.8362411488e-03  4.8362411488e-03
4.8362411488e-03  4.8362411488e-03
7.7376114656e-03  7.7376114656e-03
9.6687394780e-03  9.6687394780e-03
相対誤差の最大 < 1e-12: True
```

**検算.** (a) 2 次元の離散 Laplace 行列の固有値は $\tau_j + \tau_k$、$\tau_j = 4\sin^2(j\pi h/2)$ です（[定常反復法の章](/mathematics/numerical-linear-algebra/stationary-iterative-methods) のモデル問題の命題の証明と同じ計算。$2 - 2\cos$ ではなく $4\sin^2$ の形で計算して桁落ちを避けています）。小さい順に $(1,1)$、$(1,2)$ と $(2,1)$ の 2 重、$(2,2)$、$(1,3)$ と $(3,1)$ の片方、が並びます。(b) eigsh の結果は 10 桁一致し、相対誤差は $10^{-12}$ 未満です。2 重の固有値も 2 つとも求まっています。(c) $\sigma = 0$ のシフト反転では、最小の固有値 $1.9\times10^{-3}$ と次の固有値 $4.8\times10^{-3}$ が $A^{-1}$ では $517$ と $207$ になります。<Ref to="thm-kaniel-paige" /> の隔たりの比は、$A$ のまま最小の固有値を狙うと $(4.84 - 1.93)\times10^{-3}/(8.00 - 0.005) \approx 3.6\times10^{-4}$ ですが、$A^{-1}$ の最大の固有値では $(517 - 207)/207 \approx 1.5$ と 4 桁改善します。

<Remark id="rem-implicit-restart" title="陰的再開と大規模固有値問題の実務">
Lanczos 法をそのまま続けると、基底の保持と再直交化の手間が反復回数とともに増えます。ARPACK は、$k$ 本の Lanczos ベクトルをためるたびに、不要な Ritz 値をシフトとする陰的シフトつき QR 法を $T_k$ に施して基底を縮め、欲しい固有値の成分を残したまま再開します（Sorensen の陰的再開）。この仕組みは QR 法の章の陰的 Q の定理の応用です。D. C. Sorensen, "Implicit application of polynomial filters in a k-step Arnoldi method", *SIAM J. Matrix Anal. Appl.* 13 (1992)、R. B. Lehoucq, D. C. Sorensen, C. Yang『ARPACK Users' Guide』SIAM, 1998 を参照してください。
</Remark>

## 8. 演習

<Exercise id="exr-cg-iteration-count" difficulty="易">
正定値行列で $\kappa = 10^4$ とする。誤差の $A$ ノルムを $10^{-6}$ 倍にするのに十分な反復回数を、共役勾配法（<Ref to="thm-cg-chebyshev-bound" />）と最急降下法（<Ref to="prop-steepest-descent" />）の評価から求めよ。
<Solution>
共役勾配法は $2q^k \le 10^{-6}$、$q = (100 - 1)/(100 + 1)$ から $k \ge \ln(2\times10^{6})/\ln(101/99) \approx 14.509/0.02000 \approx 725.4$ で、$726$ 回で十分です。最急降下法は $\bigl((\kappa - 1)/(\kappa + 1)\bigr)^k \le 10^{-6}$ から $k \ge \ln 10^{6}/\ln(10001/9999) \approx 13.816/(2.000\times10^{-4}) \approx 6.9\times10^{4}$ 回です。比は約 $\sqrt{\kappa}$ 倍（ここでは約 95 倍）です。
</Solution>
</Exercise>

<Exercise id="exr-cg-low-rank" difficulty="標準">
$A = I + B$ で、$B$ は階数 $r$ の半正定値 Hermite 行列とする。共役勾配法は厳密な算術では高々 $r + 1$ 回で止まることを示せ。
<Solution>
$B$ の 0 でない固有値は高々 $r$ 個（重複を込めて）なので、$A$ の固有値は $1$ と、高々 $r$ 個の $1 + \mu_i$（$\mu_i > 0$）です。相異なる固有値の個数は高々 $r + 1$ なので、<Ref to="cor-cg-finite-termination" /> から共役勾配法は高々 $r + 1$ 回で止まります。多項式で言えば、$p(t) = (1 - t)\prod_i(1 - t/(1 + \mu_i))$ が $p(0) = 1$ を満たし $A$ のすべての固有値で $0$ になるので、<Ref to="thm-cg-chebyshev-bound" /> の真ん中の評価が $0$ になります。
</Solution>
</Exercise>

<Exercise id="exr-lanczos-two-steps" difficulty="標準">
$A = \operatorname{diag}(3, 2, 1)$、$\boldsymbol{b} = (1, 1, 1)^{\mathsf{T}}$ として Lanczos 法を 2 段行い、$T_2$ と Ritz 値 $\theta_1^{(2)}, \theta_2^{(2)}$ を求めよ。<Ref to="thm-ritz-optimality" /> の 2 の交互性と、<Ref to="thm-kaniel-paige" /> の評価が成り立っていることを確かめよ。
<Solution>
$\boldsymbol{q}_1 = (1,1,1)^{\mathsf{T}}/\sqrt{3}$、$\alpha_1 = \boldsymbol{q}_1^{\mathsf{T}}A\boldsymbol{q}_1 = 2$、$A\boldsymbol{q}_1 - 2\boldsymbol{q}_1 = (1, 0, -1)^{\mathsf{T}}/\sqrt{3}$ なので $\beta_1 = \sqrt{2/3}$、$\boldsymbol{q}_2 = (1, 0, -1)^{\mathsf{T}}/\sqrt{2}$、$\alpha_2 = (3 + 1)/2 = 2$ です。
$$
T_2 = \begin{pmatrix} 2 & \sqrt{2/3} \\ \sqrt{2/3} & 2 \end{pmatrix},\qquad \theta^{(2)} = 2 \pm \sqrt{2/3} \approx 2.8165,\ 1.1835 .
$$
交互性は $3 \ge 2.8165 \ge 2$、$2 \ge 1.1835 \ge 1$ で成り立ちます。Kaniel–Paige の評価では $\gamma = (3 - 2)/(2 - 1) = 1$、$T_1(3) = 3$、$\tan\angle(\boldsymbol{b}, \boldsymbol{e}_1) = \sqrt{2}$ なので、右辺は $(3 - 1)(\sqrt{2}/3)^2 = 4/9 \approx 0.444$ で、実際の誤差 $3 - 2.8165 \approx 0.1835$ はその内側にあります。
</Solution>
</Exercise>

## 参考文献

- G. H. Golub, C. F. Van Loan『Matrix Computations』4th ed., Johns Hopkins University Press, 2013 — §10.1（Lanczos 法）、§11.3（共役勾配法）。
- Y. Saad『Iterative Methods for Sparse Linear Systems』2nd ed., SIAM, 2003 — 第 6 章（Krylov 部分空間法）、第 9 章（前処理つき反復法）。
- A. Greenbaum『Iterative Methods for Solving Linear Systems』SIAM, 1997 — 第 2〜4 章。
- B. N. Parlett『The Symmetric Eigenvalue Problem』SIAM, 1998 — 第 12・13 章（Krylov 部分空間と Lanczos 法）。
- L. N. Trefethen, D. Bau III『Numerical Linear Algebra』SIAM, 1997 — Lectures 32–38。
- G. Meurant, Z. Strakoš, "The Lanczos and conjugate gradient algorithms in finite precision arithmetic", *Acta Numerica* 15 (2006), 471–542.