共役勾配法(CG法)

記事の内容


今回はCG法をとりあげます. 収束性の議論や前処理は扱いません. また, 細かい証明も省きます.

CG法の概要

経緯と概要

CG法(Conjugate Gradient Method; 共役勾配法)は, 連立一次方程式の数値解法で, 1952年にM.R.HestenesとE.Stiefelによって提案されました. 多くの望ましい性質を有することから, 当初から相当注目されたそうです. 元々, NBSという大規模な国家的な研究プロジェクトの中で提案されたことから, 期待も大きく, 世紀の大発見として持て囃されたとか. 戸川隼人先生[1]によれば, "数値計算法の一解法が, これほど大々的に報道されたことは, 前にも後にも例がない"とのこと. しかし不幸にも, 当初期待されていたほどの性能が発揮できなかったことから, 次第に忘れられていきます. CG法の提案から数年経った1960年ごろ頃, 非線形最適化問題の数値解法として再登場. さらに1970年頃には, 係数行列が疎の場合の収束性が注目され, その重要性が再認識されます. 1970年代後半, 前処理技術も発達し, 応用に向けて研究も活発に行われました. そして1980年頃にはKrylov部分空間法として再定式化され, 変種のアルゴリズムの開発が加速します. 今では応用でも利用される重要なアルゴリズムの1つであり, 20世紀の最も重要なアルゴリズムの1つとして知られます.

それでは, 一体なぜ, CG法はそれほどまでに注目されたのでしょうか. 一般に認識されているCG法のメリットをいくつか挙げます.

【CG法のメリット】

  • 原理的には有限回の反復で厳密回に到達する.
  • 計算量は直接法よりも少ない.
  • 行列の全成分を知らなくても, 行列-ベクトル積が分かれば計算可能. (例えばFFTなどと相性が良い)
  • fill-inを回避できる. (疎性を保ったまま計算できる)

CG方が提案された際, 特に注目されたのは1つ目の点です. 係数行列が\(n\times n\)行列の場合, 理論上は最大でも\(n\)反復以内に厳密回に到達します. しかし実際には丸め誤差の影響から, より多くの反復を必要とします. 収束しない場合もあります. 連立一次方程式\(A\boldsymbol{x}=\boldsymbol{b}\)に対するCG法の具体的なアルゴリズムを以下に示します. ただし, 係数行列は対称行列とします.

【CG法のアルゴリズム】

Set the threshold \(\varepsilon > 0\).
Initialize \(\boldsymbol{x}_0\) and set \(\boldsymbol{r}_0=\boldsymbol{b}-A\boldsymbol{x}_0\), \(\boldsymbol{p}_0=\boldsymbol{r}_0\), \(k=0\).
while \(\|\boldsymbol{r}_k\| < \varepsilon\)
\(\displaystyle \alpha_k = \frac{\langle \boldsymbol{r}_k, \boldsymbol{p}_k\rangle}{\langle \boldsymbol{p}_k, A\boldsymbol{p}_k\rangle}\)
\(\boldsymbol{x}_{k+1} = \boldsymbol{x}_k+\alpha_k \boldsymbol{p}_k\)
\(\boldsymbol{r}_{k+1} = \boldsymbol{r}_k-\alpha_k A\boldsymbol{p}_k\)
\(\displaystyle \beta_k = -\frac{\langle \boldsymbol{r}_{k+1}, A\boldsymbol{p}_k\rangle}{\langle \boldsymbol{p}_k, A\boldsymbol{p}_k\rangle}\)
\(\boldsymbol{p}_{k+1} = \boldsymbol{r}_k + \beta_k \boldsymbol{p}_k\)
\(k = k+1\)
end

行列-ベクトル積\(A\boldsymbol{p}_k\)が何度か登場しますが, これは各反復で一度計算すれば十分です. 追加の計算コストを削減できます. CG法で用いる行列Aに関する情報は, この行列-ベクトル積のみです. したがって, 行列-ベクトル積の計算方法さえ分かれば, CG法は適用できます. また, Gaussの消去法などの直接法と異なり, 計算過程で元々0の要素に非ゼロ要素が混入するfill-in現象は起こりません. したがって, JuliaやPythonであればスパース行列を用いて非ゼロ要素のみを保持しておけば計算できます. これは大規模でスパースな問題では特に有利です.

2つのストーリー

次節から, 上記アルゴリズムがどのように導出できるかを説明します. 理解の仕方は大まかに2通りあります. 1つ目は最適化問題から出発して, 適当な仮定を追加していくことでアルゴリズムを構成していく見方. 2つ目は, Krylov部分空間法という枠組みから, CG法がどのように位置付けられるかを俯瞰する見方. 次節ではまず前者のストーリーを追っていきます.

CG法の解説は, 数値解析の教科書ならば大抵載っています. [1][2][3][4][5]あたりは変種も含め詳しい記述があります. 特に, [3]は非対称行列に対するアルゴリズムなど, 話題が豊富です. 最適化問題の解法としての記述は[6]が詳しく, [7]は関数解析をベースにHilbert空間内で議論を展開しています.

CG法を構成する

この節では, CG法を具体的に構成していきます. よくある説明ですが, 勾配降下法から出発して, 途中で適当な仮定や操作を加え, アルゴリズムを組み立てていく方針です. 細かい証明は省きます. ここでは, 係数行列が正定値対称の場合を考えます.

解の特徴づけ

まずここで用いる記号を確認します. \(n\)を2以上の自然数とし, \(A\in\mathbb{R}^{n\times n}\), \(\boldsymbol{b}\in\mathbb{R}^{n}\)とします. 本節"CG法を構成する"では, \(A\)を正定値対称行列とします. また, \(\mathbb{R}^{n}\)の標準内積を\(\langle \cdot, \cdot\rangle\)で表し, 次のような内積も定めます. \(A\)の対称性から, 内積であることが分かります.

\begin{equation} \langle \boldsymbol{x}, \boldsymbol{y}\rangle_A = \langle \boldsymbol{x}, A\boldsymbol{y}\rangle \end{equation}

次の連立一次方程式を考えます. \(\boldsymbol{x}\in\mathbb{R}^n\)は, この方程式の解とします. \(A\)は正定値なので正則であり, 解が一意に存在します.

\begin{equation} A\boldsymbol{x} = \boldsymbol{b} \end{equation}

以下, この方程式の解\(\boldsymbol{x}\)の近似解列\(\{\boldsymbol{x}_k\}_{k=0}^{\infty}\)を構成する方法を考えます. まず, 解の特徴を掴むため, 次の関数を考えます.

\begin{equation} \Phi(\boldsymbol{\xi}) = \frac{1}{2}(\boldsymbol{\xi}-\boldsymbol{x})^{\mathrm{T}}A(\boldsymbol{\xi}-\boldsymbol{x}) \end{equation}

\(A\)の正定値性から, 関数\(\Phi\)は\(\boldsymbol{x}\)で唯一の最小値を持ちます. したがって, \(\Phi\)の最小値を探すことで, 連立一次方程式の解にたどり着きます. そこで, 初期解\(\boldsymbol{x}_0\)から始めて, 次のような直線探索の更新式で近似解列を構成することにします.

\begin{equation} \boldsymbol{x}_{k+1} = \boldsymbol{x}_{k}+\alpha_k \boldsymbol{p}_k \end{equation}

ここで, ステップ幅\(\alpha_k\)と探索方向ベクトル\(\boldsymbol{p}_k\)は現時点では未知です. \(\alpha_k\)は, 各\(\boldsymbol{x}_k\)において\(\boldsymbol{\Phi}\)が最小となるよう定めるのが自然です. \(\Phi\)を\(\alpha_k\)に関する関数とみて微分し, 微分係数を0とおいて整理すると, 次式が得られます.

\begin{equation} \alpha_k = \frac{\langle \boldsymbol{p}_k, \boldsymbol{r}_k\rangle}{\langle \boldsymbol{p}_k, A\boldsymbol{p}_k\rangle} \end{equation}

ただし, \(\boldsymbol{r}_k\)を次式で定め, 残差ベクトルと呼びます.

\begin{equation} \boldsymbol{r}_{k} = \boldsymbol{b}-A\boldsymbol{x}_{k} \end{equation}

あとは, 探索方向\(\boldsymbol{p}_k\)を定めるだけです. 最急降下法の場合, 負の勾配 \begin{equation} -\nabla\boldsymbol{\Phi}(\boldsymbol{x}_k)=\boldsymbol{r}_k \end{equation} とします. つまり, 残差ベクトル方向に進めば, 関数値が最も減少します. ただし, 最急降下法は精度があまり良くないので, 違う方向に進むことにします. 残差ベクトル方向は良い候補ですが, ここでは修正を加えます.

直交化

探索方向ベクトルは, 残差ベクトル方向に少し修正を加えて構成します. 探索はなるべく少ない回数で完了した方が良さそうです. そこで, 探索方向ベクトルは互いに1次独立であるとします. これは, 一度探索した方向は二度探索しないための要請です. さらに, 直交性も課しておけば, 探索の効率性も期待できます. なお, 直交性を仮定しておけば1次独立性は保証されるので, ここでは直交性のみを要請します. この要請を満たすために, 探索方向\(\{\boldsymbol{p}_k\}_{k=0}^{n-1}\)を次のように構成します. これが探索方向ベクトルの更新式になります.

\begin{align} \boldsymbol{p}_0 &= \boldsymbol{r}_0\\ \boldsymbol{p}_k &= \boldsymbol{r}_k - \sum_{i=0}^{k-1}\frac{\langle\boldsymbol{r}_k, \boldsymbol{p}_i \rangle_A}{\langle \boldsymbol{p}_i, \boldsymbol{p}_i\rangle_A}\boldsymbol{p}_i,\quad (k=1,\cdots n-1) \end{align}

これは, 内積\(\langle \cdot, \cdot\rangle_A\)に関するGram-Shmidtの直交化です. すなわち, 残差ベクトル方向にGram-Shmidtの直交化法を適用することで, 探索方向ベクトルを更新できます. しかし, 更新式の右辺は少々複雑です. 次の補題1を用いると簡単な式になります. 証明は, [2]や[3]にあります.

補題1

自然数\(k=2,3,\cdots, n-1\), 自然数\(i=0,1,\cdots, k-2\)に対して, 次式が成り立つ.

\begin{equation} \langle \boldsymbol{r}_k, \boldsymbol{p}_i \rangle_A = 0 \end{equation}

以上より, 探索方向は次式のように表されます.

\begin{align} \boldsymbol{p}_0 &= \boldsymbol{r}_0\\ \boldsymbol{p}_k &= \boldsymbol{r}_k - \frac{\langle\boldsymbol{r}_k, \boldsymbol{p}_{k-1} \rangle_A}{\langle \boldsymbol{p}_{k-1}, \boldsymbol{p}_{k-1}\rangle_A}\boldsymbol{p}_{k-1},\quad (k=1,\cdots n-1) \end{align}

以上の議論により, CG法のアルゴリズムが導出されます. もう一度CG法のアルゴリズムを見返してみると, ステップ幅の更新, 直線探索による近似解の更新, 探索ベクトル方向の直交化という流れになっています. なお, 以上のように解を更新すると, 残差ベクトルは標準内積の意味で直交します[2].

補題2

自然数\(i,j(i\neq j)\)に対して, 次式が成り立つ.

\begin{equation} \langle \boldsymbol{r}_i, \boldsymbol{r}_j \rangle= 0 \end{equation}

したがって, 最初のn反復で残差ベクトルはn本のベクトルからなる直交系をなします. n+1反復目の残差ベクトルはさらにその直交系と直交します. 一方, \(\mathbb{R}^n\)はn次元なので, n+1反復目の残差ベクトルは自動的に零ベクトルになります. したがって, n反復で残差は0に達し, 真の解にたどり着きます.

ちなみに関数\(\Phi\)はエネルギーの形をしており, CG法はエネルギーを効率よく最小化するアルゴリズムとも言えます.

CG法を俯瞰する

この節では, Krylov部分空間法という枠組み方CG法を俯瞰します. 前節のようにアルゴリズムを組み立てていくというよりは, より大きな視点から出発して, いくつかの仮定をおき, CG法がどのような位置付けであるかを追っていく方針です. 細かい証明は省きます. ここでは, 係数行列が対称の場合を考えます.

Krylov部分空間

本節で用いる記号は前節と同じです. ただし, 本節"CG法を俯瞰する"では, 行列Aは対称行列とし, 必ずしも正定値とは限りません. 冒頭からKrylov部分空間という名前だけは出していますが, 未定義でした. まずは以下に定義を述べます.

定義【Krylov部分空間】

非ゼロベクトル\(\boldsymbol{v}\in\mathbb{R}^{n}\), 行列\(A\in\mathbb{R}^{n\times n}\), 自然数\(m(\leq n)\)に対して, \(\boldsymbol{v}, A\boldsymbol{v}, A^2\boldsymbol{v}, \cdots,A^{m-1}\boldsymbol{v}\)が1次独立のとき, \begin{equation} K_m(A;\boldsymbol{v}) = \mathrm{span} \{\boldsymbol{v}, A\boldsymbol{v}, A^2\boldsymbol{v}, \cdots,A^{m-1}\boldsymbol{v}\} \end{equation} を行列\(A\)とベクトル\(\boldsymbol{v}\)に関する\(m\)次のKrylov部分空間という.

Krylov部分空間を用いた枠組みは, Krylov部分空間法と呼ばれます.

定義【Krylov部分空間法】

Krylov部分空間の何らかの意味での直交系を構成し, 次式によって連立一次方程式\(A\boldsymbol{x}=\boldsymbol{b}\)の近似解列\(\{\boldsymbol{x}_k\}_{k=0}^{\infty}\)を構成する方法を一般にKrylov部分空間法という.

\begin{equation} \boldsymbol{x}_k = \boldsymbol{x}_{0} + \boldsymbol{z}_k,\quad \boldsymbol{z}_k\in K_k(A;\boldsymbol{r}_0) \end{equation}

行列の性質にもよりますが, Krylov部分空間は次数が大きくなると集合として大きくなります. 上の定義で, \(\boldsymbol{z}_k\)は初期解からの正味の移動を表します. ステップを踏むごとに正味の移動を表すベクトルが探索できる範囲が拡大するので, 効率的な探索が期待できます.

直交化

CG法において, 直交化の操作は本質的です. Krylov部分空間法の枠組みでも, 直交化が解法を特徴付けます. Krylov部分空間では, 次のような直交化法を用います. 1つのベクトルをインプットとし, AをかけるごとにGram-Schmidtと同じ要領で直交化します. インプットである\(\boldsymbol{v}\in\mathbb{R}^n\)は, 任意のベクトルとします. アウトプットであるベクトル集合\(\{\boldsymbol{v}_1,\cdots,\boldsymbol{v}_n\}\)は標準内積の意味でKrylov部分空間の正規直交系です.

【Lanczos過程】

Initialize \(\displaystyle \boldsymbol{v}_1 = \frac{\boldsymbol{v}}{\|\boldsymbol{v}\|}\).
for \(k=1,2,\cdots,n-1\)
\(\displaystyle \boldsymbol{v}_{k+1}^{\prime} = A\boldsymbol{v}_{k} - \frac{\langle A\boldsymbol{v}_{k}, \boldsymbol{v}_{k}\rangle}{\langle \boldsymbol{v}_{k}, \boldsymbol{v}_{k}\rangle}\boldsymbol{v}_{k}-\frac{\langle A\boldsymbol{v}_{k}, \boldsymbol{v}_{k-1}\rangle}{\langle \boldsymbol{v}_{k-1}, \boldsymbol{v}_{k-1}\rangle}\boldsymbol{v}_{k-1}\)
\(\displaystyle \boldsymbol{v}_{k+1} = \frac{\boldsymbol{v}_{k+1}^{\prime}}{\|\boldsymbol{v}_{k+1}^{\prime}\|}\)
end

Lanczos過程は, Arnoldi過程と呼ばれるものの特別な場合です. Arnoldi過程はGram-Schmidtと似た形をしています. Arnoldi過程を対称行列に適用すると, Lanczos過程になります. これらの方法は, 固有値の数値計算でも用いられます.

CG法の位置づけ

Krylov部分空間法の枠組みでCG法を考えます. まず, 初期残差ベクトル\(\boldsymbol{r}_0\)から始まるLanczos過程を用いて, 残差ベクトル集合\(\{\boldsymbol{r}_0, \cdots, \boldsymbol{r}_{n-1}\}\)がKrylov部分空間の標準内積の意味で直交基底をなすようにします. これは前節での補題2に対応します. さらに, 次を仮定します. これはRitz-Galerkin条件と呼ばれます.

仮定【直交条件】

各\(k=1,\cdots,n-1\)に対して, 次式を仮定する.

\begin{equation} \boldsymbol{r}_k \bot K_k(A;\boldsymbol{r}_0) \end{equation}

この条件は, 前節での補題1に対応します. これらを用いて計算することで, CG法のアルゴリズムが導出できます[4]. 計算は少々混み入っています. 逆に, どのような意味でKrylov部分空間の直交基底を構成するか, 直交条件の代わりにどのような条件を仮定するかで, 導出されるアルゴリズムは異なります. それら変種のアルゴリズムは, [3][4]あたりに詳しく書いてあります.

CG法を使ってみる

この節では, CG法を実装して計算してみます.

実装

実装の前に, 多少の式変形によって, 次のように書けます.

\begin{equation} \alpha_k = \frac{\langle \boldsymbol{p}_k, \boldsymbol{r}_k\rangle}{\langle \boldsymbol{p}_k, A\boldsymbol{p}_k\rangle} = \frac{\langle \boldsymbol{r}_k, \boldsymbol{r}_k\rangle}{\langle \boldsymbol{p}_k, A\boldsymbol{p}_k\rangle} \end{equation} \begin{equation} \beta_k = -\frac{\langle \boldsymbol{p}_k, A\boldsymbol{r}_{k+1}\rangle}{\langle \boldsymbol{p}_k, A\boldsymbol{p}_k\rangle} = \frac{\langle \boldsymbol{r}_{k+1}, \boldsymbol{r}_{k+1}\rangle}{\langle \boldsymbol{r}_k, \boldsymbol{r}_k\rangle} \end{equation}

今回は, こちらを用います. また, 今回書いたプログラムは素朴で, あまり良いコードではありません(もちろんいつものことですが...). 次の連立一次方程式を考えます.

\begin{equation} \begin{bmatrix} 2 & 3 & 1\\ 3 & -1 & 4\\ 1 & 4 & 1\end{bmatrix} \boldsymbol{x} = \begin{bmatrix} 1 \\ 3 \\ 2 \end{bmatrix} \end{equation}

最大の反復回数は1000回, 収束判定は残差ノルムが\(10^{-6}\)を下回った時点で終了としました. 各反復の残差ノルムの値を以下に示します. ちなみに解は, [-0.6666666666666674, 0.3333333333333345, 1.333333333333334]となりました.

【コード3の実行結果】

コード

【Juliaコード1; インポート】
using Plots
using LinearAlgebra
pyplot()
【Juliaコード2; インポート】
n = 3
A = [
    2. 3. 1.;
    3. -1. 4.;
    1. 4. 1.
]
b = [1., 3., 2.]
x_true = A\b

#行列-ベクトル積
function multAvec(A::Array{Float64,2}, vec::Array{Float64,1})
    return A*vec
end
multA(vec::Array{Float64,1}) = multAvec(A,vec)

#CG法
function CG(b::Array{Float64, 1}, x₀::Array{Float64, 1}, maxiter::Int64, ϵ::Float64)
    #探索方向ベクトル, 残差ベクトル, 行列-ベクトル積, 近似解
    r = b- multA(x₀)
    p = r
    Ap = multA(p)
    x = x₀

    #残差ベクトルのノルム
    normr_vec = zeros(maxiter+1)
    normr = norm(r)
    normr_sq = normr^2
    normr_vec[1] = normr
    if normr < ϵ
        return x, normr_vec[1:k+1]
    end

    #メインループ
    for k in 1:maxiter
        #更新
        α = normr_sq/dot(p,Ap)
        x = x + α*p
        r = r - α*Ap
        normr = norm(r)
        normr_vec[k+1] = normr

        #収束のチェック
        if normr < ϵ
            return x, normr_vec[1:k+1]
        end

        β = normr^2/normr_sq
        p = r + β*p
        Ap = multA(p)
        normr_sq = normr^2
    end
    return x, normr_vec
end
【Juliaコード3; インポート】
x₀ = zeros(n)
maxiter = 1000
ϵ = 1e-6
x, normr_vec = CG(b, x₀, maxiter, ϵ)
println(x)

fig = plot( title="Residual norm", xlabel="step", ylabel="norm of residual", leg=false)
plot!(0:length(normr_vec)-1, normr_vec, yscale=:log10, marker=:auto)
savefig(fig, "figs-CG/fig1.png")
参考文献

      [1] 戸川隼人, 共役勾配法, 教育出版, 1977
      [2] 森正武, 数値解析, 共立出版, 2002
      [3] 杉原正顕, 室田一雄, 線形計算の数理, 岩波書店, 2016
      [4] 藤野清次, 張紹良, 反復法の数理, 1996
      [5] R.Barrett, T.F.Chan , J.Donato, M.Berry, J.Demmel, (長谷川里美, 長谷川秀彦, 藤野清次 訳), 反復法Templates, 朝倉書店, 1996
      [6]D.Ruenberger, (増渕正美, 嘉納秀明 訳), 関数解析による最適理論, コロナ社, 1973
      [7]J.Nocedal, S.Wright, Numerical Oprimization , Springer, 2009