線型方程式 $Ax=b$ への解法としての共役勾配法のメモ．
問題設定と最適化問題への変換
線型方程式
\begin{equation}
Ax=b
\tag{1}
\end{equation}
を考える．ここで，$A\in\mathbb{R}^{n\times n}$ は正定値対称行列とする．
次の最小化問題の最適解は線型方程式 (1) の解となっている．
\begin{equation}
\mathop{\text{minimize}}_{x\in\mathbb{R}^n} \quad f(x) := \frac12x^\top Ax - b^\top x
\tag{2}
\end{equation}
これは，$A$ が正定値対称行列であるため，$f(x)$ が狭義凸関数で $\nabla f(x) = Ax - b$ となることからわかる．
以降，最適化問題 (2) について考える．
最適化問題に対する反復法
$f(x)$ を反復法により最小化することを考える．つまり，点列 ${x^0, x^1, \ldots}$ を
x^{k+1} = x^k + t_k d^k
により生成することを考える．ここで，$t_k \in \mathbb R$ はステップ幅，$d^k \in \mathbb R^n$ は探索方向である．
$t_k$ を正確な直線探索により計算する，つまり
t_k = \mathop{\text{arg min}}_{t \in \mathbb R} f\left(x^k + t d^k\right)
で計算するとすると，$\frac{\mathrm{d}}{\mathrm{d}t} f\left(x^k + t d^k\right) = 0$ より，
t_k = -\frac{\nabla f\left(x^k\right)^\top d^k}{\left(d^k\right)^\top A d^k}
となる．
Cholesky 分解を用いた問題の変形と探索方向の考察
問題の Cholesky 分解を用いた変形
正定値対称行列 $A$ が $A = LL^\top$ と下三角行列 $L\in\mathbb R^{n \times n}$ を用いて Cholesky 分解されるとする．新たに $y := L^\top x$ とすると，問題 (2) は
\mathop{\text{minimize}}_{x \in \mathbb R^n} \quad \frac12 x^\top LL^\top x - b^\top L^{-\top} L^\top x = \frac12 y^\top y -\left(L^{-1} b\right)^\top y
と書き換えられる．$b' = L^{-1}b$ とすると，$g(y) := \frac12 y^\top y - \left(L^{-1} b\right)^\top y = \sum_{i=1}^n \left(\frac12y_i^2 - b'_iy_i\right)$ は成分ごとに最小化することができる．
言い換えれば，各単位ベクトル $e^1, e^2, \ldots, e^n$ を探索方向とすることで，$y^n$ を最適解とする数列 ${y^0, y^1, \ldots, y^n}$ を
t_k = \mathop{\text{arg min}}_{t \in \mathbb R} g\left(y^k + t e^{k+1}\right),\quad y^{k+1} = y^k + t_k e^{k+1}
により生成することができる．
問題の直交変換を用いた変形
さらに，直交行列 $P$ による直交変換 $z = Py$ (つまり，$y = P^\top z$) を考える．
\mathop{\text{minimize}}_{x \in \mathbb R^n} \quad \frac12 y^\top y -\left(L^{-1} b\right)^\top y = \frac12 z^\top z - \left(PL^{-1}b\right)^\top z
より，先ほどと同様に $z$ の各成分について最小化できることがわかる．
つまり，$s^k = P^\top e^{k+1}$ とすると，直交するベクトル列 $s^0,s^1,\ldots,s^{n-1}$ を探索方向とすることで，$y^n$ を最適解とする数列 ${y^0, y^1, \ldots, y^n}$ を
t_k = \mathop{\text{arg min}}_{t \in \mathbb R} g\left(y^k + t s^k\right),\quad y^{k+1} = y^k + t_k s^k
により生成することができる．
元の問題への帰着
ここで，$y = L^\top x$ に注意して元の問題に適用することを考える．
$y^{k+1} = y^k + t_k d^k$ の両辺の左から $L^{-\top}$ をかけて整理すると，
x^{k+1} = x^k + t_k L^{-\top} s^k
となる．新たに $d^k = L^{-\top} s^k$ とすると，$s^k$ の直交条件は，
0 = \left(s^k\right)^\top s^k = \left(d^k\right)^\top LL^\top d^k = \left(d^k\right)^\top A d^k
とかける．内積 $\langle a, b\rangle_A := a^\top A b$ について直交することを $A$-共役という．
$A$-共役なベクトル列 ${d^0, d^1, \ldots, d^{n-1}}$ を探索方向とすることで，$x^n$ を最適解とする数列 ${x^0, x^1, \ldots, x^n}$ を
t_k = \mathop{\text{arg min}}_{t \in \mathbb R} f\left(x^k + t d^k\right),\quad x^{k+1} = x^k + t_k d^k
により生成することができる．
A-共役な探索方向の計算
最急降下法における探索方向 $-\nabla f\left(x^0\right), -\nabla f\left(x^1\right), \ldots, -\nabla f\left(x^{n-1}\right)$ を Gram-Schmidt の直交化法により A-共役にすることを考える．つまり，
\begin{align}
d^0 &= -\nabla f\left(x^0\right)\\
d^1 &= -\nabla f\left(x^1\right) + \frac{\left\langle \nabla f\left(x^1\right), d^0 \right\rangle_A}{\left\langle d^0, d^0 \right\rangle_A} d^0\\
&\vdots\\
d^{n-1} &= -\nabla f\left(x^{n-1}\right) + \sum_{i=0}^{n-2} \frac{\left\langle \nabla f\left(x^{n-1}\right), d^i \right\rangle_A}{\left\langle d^i, d^i \right\rangle_A} d^i
\end{align}
を考える．$k = 0, 1, \ldots, n-1$ について $\nabla f\left(x^k\right) \neq 0$ を仮定すると，$d^k \neq 0$ であり，
$$\nabla f\left(x^k\right)^\top \nabla f\left(x^\ell\right) = \nabla f\left(x^k\right)^\top d^\ell = 0\quad (\ell = 0, 1, \ldots, k-1)$$
となることを示すことができる．これは，正確な直線探索より $\nabla f\left(x^{k+1}\right)^\top d^{k} = \left.\frac{\mathrm{d}}{\mathrm{d} t}f\left(x^{k} + td^{k}\right)\right|_{t=t_k} = 0$ となることを用いて，数学的帰納法より証明することができる．
また，
\begin{align}
\left\langle \nabla f\left(x^k\right), d^\ell \right\rangle_A &= \frac1{t_\ell}\left\langle \nabla f\left(x^k\right), x^{\ell+1} - x^{\ell} \right\rangle_A\\
&= \frac1{t_\ell}\left\langle \nabla f\left(x^k\right), Ax^{\ell+1} - Ax^{\ell} \right\rangle\\
&= \frac1{t_\ell}\left\langle \nabla f\left(x^k\right), Ax^{\ell+1} - b - \left(Ax^{\ell} - b\right)\right\rangle\\
&= \frac1{t_\ell}\left(\nabla f\left(x^k\right)^\top \nabla f\left(x^{\ell+1}\right) - \nabla f\left(x^k\right)^\top \nabla f\left(x^\ell\right)\right)\\
&= \begin{cases}
\frac1{t_{k-1}} \nabla f\left(x^k\right)^\top \nabla f\left(x^k\right) & (\ell = k-1)\\
0 & (\ell \leq k-2)
\end{cases}\\
&= \begin{cases}
-\frac{\left\langle d^{k-1}, d^{k-1}\right\rangle_A}{\nabla f\left(x^{k-1}\right)^\top d^{k-1}}\nabla f\left(x^k\right)^\top \nabla f\left(x^k\right) & (\ell = k-1)\\
0 & (\ell \leq k-2)
\end{cases}
\end{align}
より，
$$d^{k+1} = -\nabla f\left(x^{k+1}\right) + \beta_k d^k \quad (k\geq0)$$
が成り立つ．ここで，
$$\beta_k = - \frac{\nabla f\left(x^{k+1}\right)^\top \nabla f\left(x^{k+1}\right)}{\nabla f\left(x^{k}\right)^\top d^{k}}$$
とした．また，$\beta_k$ の分母について
$$\nabla f\left(x^{k}\right)^\top d^{k} = \nabla f\left(x^{k}\right)^\top \left(-\nabla f\left(x^{k}\right) + \beta_{k-1} d^{k-1}\right) = -\nabla f\left(x^{k}\right)^\top \nabla f\left(x^{k}\right)$$
となるので，
$$\beta_k = \frac{\nabla f\left(x^{k+1}\right)^\top \nabla f\left(x^{k+1}\right)}{\nabla f\left(x^{k}\right)^\top \nabla f\left(x^{k}\right)}$$
が得られる．
まとめ
線型方程式 $Ax = b$ の解は以下の手順で求められる．
- $x^0$ を初期化する．
- $d^0 = r^0 = -\nabla f\left(x^0\right) = b - Ax^0$ で探索方向と残差を初期化する．
- for $k = 0, 1, \ldots, n-1$ do
- ステップ幅 $t_k = \frac{\left(r^k\right)^\top d^k}{\left(d^k\right)^\top A d^k}$ を計算する．
- $x^{k+1} = x^k + t_k d^k$ で更新する．
- $r^{k+1} = r^k - t_kAd^k$ で残差を更新する．
- $\beta_k = \frac{\left(r^{k+1}\right)^\top r^{k+1}}{\left(r^{k}\right)^\top r^{k}}$ を計算する．
- $d^{k+1} = r^{k+1} + \beta_k d^k$ で探索方向を更新する．
- end for
また，このアルゴリズムを一般の非線形関数 $f$ に拡張することができる．
- $x^0$ を初期化する．
- $d^0 = -\nabla f\left(x^0\right)$ で探索方向を初期化する．
- for $k = 0, 1, \ldots, n-1$ do
- ステップ幅 $t_k$ を計算する．
- $x^{k+1} = x^k + t_k d^k$ で更新する．
- $\beta_k = \frac{\nabla f\left(x^{k+1}\right)^\top \nabla f\left(x^{k+1}\right)}{\nabla f\left(x^{k}\right)^\top \nabla f\left(x^{k}\right)}$ を計算する．
- $d^{k+1} = -\nabla f\left(x^{k+1}\right) + \beta_k d^k$ で探索方向を更新する．
- end for
おまけ
係数 $\beta_k$ は様々なバリエーションが提案されている．ここでは，主要な 6 種類を紹介する．
$y^k = \nabla f\left(x^{k+1}\right) - \nabla f\left(x^{k}\right)$ とする．
- $\beta_k^\text{FR} = \frac{\nabla f\left(x^{k+1}\right)^\top \nabla f\left(x^{k+1}\right)}{\nabla f\left(x^{k}\right)^\top \nabla f\left(x^{k}\right)}$
- $\beta_k^\text{HS} = \frac{\nabla f\left(x^{k+1}\right)^\top y^k}{\left(d^k\right)^\top y^k}$
- $\beta_k^\text{PR} = \frac{\nabla f\left(x^{k+1}\right)^\top y^k}{\nabla f\left(x^{k}\right)^\top \nabla f\left(x^{k}\right)}$
- $\beta_k^\text{CD} = \frac{\nabla f\left(x^{k+1}\right)^\top \nabla f\left(x^{k+1}\right)}{-\nabla f\left(x^{k+1}\right)^\top d^k}$
- $\beta_k^\text{LS} = \frac{\nabla f\left(x^{k+1}\right)^\top y^k}{-\nabla f\left(x^{k+1}\right)^\top d^k}$
- $\beta_k^\text{DY} = \frac{\nabla f\left(x^{k+1}\right)^\top \nabla f\left(x^{k+1}\right)}{\left(d^k\right)^\top y^k}$
これらは，$f(x) := \frac12x^\top Ax - b^\top x$ の場合には等価である．