参考文献

1. Wendland, H. (2017). Numerical Linear Algebra: An Introduction (Cambridge Texts in Applied Mathematics). Cambridge: Cambridge University Press. doi:10.1017/9781316544938

線形連立方程式と最小化問題

線形連立方程式

を数値的に解く方法は，色々あるが，Aが対称正定値行列の場合は，共役勾配法という方法が使える．

共役勾配法とは，

の最小化問題を解く方法である．

どういうことかというと，

f(x) を微分してみると，

である．ただし，\nabla f = (\partial f/ \partial x_1, \dots, \partial f/ \partial x_N)^{\mathrm{T}}，Hf = (\partial^2 f / \partial x_i \partial x_j)_{ij} のことである．つまり勾配とヘッセ行列．

これを使って，f(x)をyの周りでテイラー展開してみると，

となって，y が Ax=b の解x^*の時にfが最小となることがわかる．
ここで，最後の不等式は，Aが正定値だからである．

つまり，最小化問題

の解x^*と，線形連立方程式

の解x^* は同じであるということである．

共役勾配法とは

ということで，タイトルに戻って，

共役勾配法とは，

を解く方法であって，求まった解xは，線形連立方程式

の解でもある．

具体的には，反復法を用いる．

つまり，

を満たす点列 x_1, x_2, \dots を構成する．

そして，f(x_{M+1}) = f(x_M) となるMがあったとき，x_Mが求めたかった x となる．

共役勾配法の詳細

具体的に点列 \{x_j\}_{j=1,2,\dots} を構成していく．

と表しておいて，\alpha_j と p_j を決める．

さて，f(x)の形を思いだすと，f(x)は，下に凸であるから，

と表しておくと，

を満たす \alpha を \alpha_jとすれば，x_jとp_jが与えられているときに，もっとも小さい f(x_{j+1}) を得ることができる．

したがって，

すなわち，

となる．

ここで，いくつか記号を導入しておく．

は，解x^* との誤差を表すベクトルである．

また，任意のx, y \in \mathbb{R}^N に対して，Aに関する内積を，

とくに，Aが単位行列Iのときの内積は，\langle x, y \rangle_I = \langle x,y\rangle_2 と表す．
これは，いつもの2乗して足していく内積である．
また，

が成り立つ．

これらの記号をつかうと，\alpha_jは，

と表せる．

探索方向 p_j

残るは，探索方向を決めるベクトル p_j を決定することである．

これは，A共役であるベクトルの組 p_0, \dots, p_{N-1} が選ばれる．
ここで，p_0, \dots, p_{N-1} が A共役であるとは，0 \leq i \not= j \leq N-1に対して，

であることを言う．

つまり，Aに関する内積による直交基底を選べ，と言うことである．

なぜ，A共役を選ぶかと言うと，上で選んだ \alpha_jに対して，解 x^* が，

と表すことができるからである．
すなわち，このように\alpha_j, p_j を選んでおくと，最悪でも，N 回で解に辿り着くことができる．

p_j の具体的な構成方法

では，p_j は具体的にはどのように選べば良いか．

A共役に選べば，最大で N 回の探索で十分であったが，できる限り探索回数は少なくしたい．
また，A 共役なベクトルの組 p_0, \dots, p_{N-1} は，最初にすべて計算しておくのではなく，反復するごとに，その都度計算した方が，計算効率がよい．

そこで，ここでは，p_0, \dots, p_i から p_{j+1} を構成する方法を述べる．
つまり，x_{j}が終了条件を満たさなかった時に，次の探索方向p_{j}を決める方法を述べる．
この場合，使用できるのは，

• これまでの探索方向，p_0,\dots, p_{j-1}，
• これまでの解の候補 x_0, \dots, x_j，
• 誤差 r_0, \dots, r_{j}

である．

少し天下り的であるが，以下のように表してみる．

p_0, \dots, p_{N-1} は A共役であるので，0 \leq l \leq j-1 に対して，

であるから，

このように，p_j を構成するメリットは，実は，\beta_j^{(0)} = \cdots = \beta_j^{(j-2)} = 0 となることである． つまり，\beta_j = \beta_j^{(j-1)} とすると，

となって，これはメモリー使用量的にもとても良い形である．

さらに，この形は，他の探索方向の候補よりも効率がよい．次はこれをより正確に述べよう．

Krylov空間

より正確に述べるために，いくつか記号を定義する．

Krylov空間とは，ある行列A \in M_n(\mathbb{R})とベクトルp \in \mathbb{R}^Nに対して，i 個のベクトルの組 p, Ap, \dots, A^{i-1}pが生成する部分空間

のことである．

上記で構成した，r_0, \dots, r_{N-1} と p_0, \dots, p_{N-1} に関して，

が成り立ち，\dim \mathcal{K}_i(A, r_0) = i である．

また，x-y \in \mathcal{K}_i(A, p) のことを，x \in y + \mathcal{K}_i(A, p)と書く．

これらの記号を使うと，j番目の解候補 x_j は，x_0 + \mathcal{K}_j(A,r_0)の元であると言える．
この中の一つが，上で構成した共役勾配法による x_jということになる．

解との距離

2つのベクトル x, y の距離を先ほど導入したAによる内積によって作られるノルム|| \cdot ||_Aによって定める．つまり，

として，x と y の距離を，||x-y||_A，と定める．

共役勾配法によって得られる，j回目の反復のときの解の候補 x_j は，他の解の候補 x \in x_0 + \mathcal{K}_j(A, r_0) よりも，真の解 x^* に近い．つまり，

である．

最大の反復回数

共役勾配法の最大の反復回数は，行列 A の相異なる固有値の個数dである．

若干の変形

数値計算しやすいように，表示を少し変えておく．

ここで，数値計算しやすい形とは，

• 行列とベクトルの積を含む内積，\langle \cdot \rangle_A よりも \langle \cdot \rangle_2 が現れるような形
• 内積よりもノルムで計算できるような形

という感じである．

今，

であって，これらを変形しておく．

最後の等号は\alpha_jを決めるときの等式(\phi^\prime(\alpha) = 0)による．

また，

より，

となる．

共役勾配法のアルゴリズム

ここまでの話をアルゴリズムにすると以下のようになる．

Choose x[0] and epsilon
Set p[0] = r[0] = b - Ax[0]
for j=0 to j=N-1 do
    t[j] = Ap[j]
    alpha[j] = ||r[j]||^2/<t[j],p[j]>
    x[j+1] = x[j] + alpha[j] * p[j]
    r[j+1] = r[j] - alpha[j] * t[j]
    beta[j+1] = ||r[j+1]||^2/||r[j]||^2
    p[j+1] = r[j+1] + beta[j+1] * p[j]
    if ||r[j+1|| < epsilon then
        break
    j++

見にくいな．良い書き方はないか．

この後の展開

この後の展開としては，以下がある．

一般の行列に対する反復解法

行列 A は対称正定値行列に限られていたので，これを一般の行列Aにまで拡張させようというもの

CGNR法

Ax = b の両辺に転置A^\mathrm{T}をかけて

をとすれば，CG法が使えるぞっという方法．

であるから，結局 r_j を最小化させていることとなる．
これより，Conjugate Gradient on th Normal equation to minimise the Residual と呼ばれる．

GMRES and MINRES

上の発想から，はじめから

を考えようというものが，GMRESとMINRES
まだよくわかっていない．スクラップを改めて書いていく予定．

制約付きの最小化問題

制約がついている場合にも対応可能．

リーマン多様体上の共役勾配法

制約を\mathbb{R}^Nに埋め込まれた多様体として表現しておいて，この多様体上で共役勾配法をすれば，無制約の最小化問題に帰着できるというもの．

これがやりたくてこのスクラップを書いている．