連立一次方程式を数値的に解く
連立一次方程式を計算機で解く方法を解説する。 これは齊藤宣一『数値解析』の第3章の読書ノートである。
3.1 ガウスの消去法とLU分解
線形代数で勉強するように、連立一次方程式はガウスの消去法（掃き出し法）で解ける。ガウスの消去法は行基本変形によって正則行列を上三角行列に変形する。
ガウスの消去法が行の入れ替えなしで行える場合について考える。この場合、行列を上三角行列に変形するために必要な変形は、\(i < j\)に対して第\(i\)行の\(a\)倍を第\(j\)行に足す変形のみである。これを表す行列は下三角行列になっている。下三角行列の積は下三角行列なので、結局、ガウスの消去法は
\[L'A = U\]という操作をしていることになる。ただし、\(A\)は連立方程式の係数行列、\(L'\)は行基本変形の積を表す下三角行列、\(U\)はガウスの消去法によって得られる上三角行列である。\(L'\)が正則で、逆行列も下三角行列であることは容易にわかるので、
\[A = LU\]と表すことができる。ここで、\(L = L'^{-1}\)は\(L'\)の逆行列で、下三角行列である。正方行列\(A\)を下三角行列\(L\)と上三角行列\(U\)をもちいて\(A=LU\)と表すことを、LU分解という。
連立方程式の係数行列がLU分解できると、連立方程式\(Ax = b\)が小さい計算量で解ける。連立方程式\(Ax = b\)を2つの連立方程式\(Ly = b\)と\(Ux = y\)に分解する。連立方程式\(Ly = b\)は上の行から順に解き、\(Ux = y\)は下の行から順に解けばよい。
3.2 LU分解の十分条件
ここまでの議論では、ガウスの消去法が行の交換なしで実現できる場合について述べた。一般の場合には、ガウス消去法では行の交換が必要である。行の交換なしでガウスの消去法が行えるための十分条件について述べる。
定義（首座小行列）
\(C=(c_{ij}) \in \mathbb{R}^{n\times n},\, 1 \leq k \leq n-1\)に対して
\[C_k := \left(\begin{array}{ccc} c_{11} & \dots & c_{1k} \\ \vdots & \ddots & \vdots \\ c_{k1} & \dots & c_{kk} \end{array}\right) \in \mathbb{R}^{k\times k}\]を\(C\)の\(k\)次首座小行列という。また、その行列式を首座小行列式という。
定理3.4
正則行列\(A\)が（行交換なしで）LU分解ができるための必要十分条件は、任意の\(k\)に対して\(k\)次の首座小行列式が\(0\)でないことである。
（\(\because\) 第\(i\)行の\(a\)倍を第\(j\)行に足す変形で行列式は変化せず、上三角行列の行列式が\(0\)でないための必要十分条件はすべての対角成分が\(0\)でないことだから）
定理
正定値対称行列はLU分解できる。
証明
正定値対称行列の行列式は\(0\)でなく、正定値対称行列の任意の首座小行列は正定値対称行列だから。
定義3.7（優対角行列）
\(A=(a_{ij})\in\mathbb{R}^{n\times n}\)が優対角行列であるとは、すべての\(i\)に対して\(|a_{ii}| \geq \sum_{j\neq i} |a_{ij}|\) であることをいう。また、不等号が\(\geq\)でなく\(>\)で成り立つとき、狭義優対角行列であるという。
\(A=(a_{ij})\in\mathbb{R}^{n\times n}\)が一般化狭義優対角行列であるとは、対角成分が正の対角行列\(D\)が存在して\(AD\)が狭義優対角行列になることをいう。
命題3.9
一般化狭義優対角行列は正則である。
証明
\(\det(AD) = \det(A)\det(D)\)であるから、\(A\)自身が狭義優対角行列である場合についてのみ示せばよい。\(A\)を狭義優対角行列とする。
\(A\)が正則でないと仮定する。すると、\(0\)でないベクトル\(x\)が存在して、\(Ax=0\)となる。
いま、\(|x_i|\)が最大になる \(i\) をとると、\(Ax = 0\) の第 \(i\) 行について
\[\sum_{j} a_{ij}x_j = 0\]だから、
\[|a_{ii}x_i| = |\sum_{j\neq i} a_{ij}x_j| \leq \sum_{j\neq i} |a_{ij}||x_j|\leq \sum_{j\neq i}|a_{ij}||x_i|\]なので、
\[|a_{ii}| \leq \sum_{i\neq j} |a_{ij}|\]となってしまい、\(A\)が狭義優対角行列であるという仮定に反する。したがって、\(A\)は正則である。
命題3.10
一般化狭義優対角行列の首座小行列は一般化狭義優対角行列である。
証明
\(A\)が一般化狭義優対角行列で、\(AD\)が狭義優対角行列であるとする。すると、首座小行列の積\(A_kD_k\)が狭義優対角行列になる。よって、\(A_k\)が一般化狭義優対角である。
定理
一般化狭義優対角行列はLU分解できる。
（\(\because\) 命題3.9と命題3.10からしたがう）
3.3 ガウスの消去法とLU分解（一般の場合）
定理（PA=LU分解）
正則行列\(A\)は、一般には、置換行列（ガウスの消去法における行交換に対応）\(P\),下三角行列\(L\)、上三角行列\(U\)を用いて、
\[PA=LU\]と分解できる。
3.4 共役勾配法
共役勾配法では、正定値対称行列\(A\)に対する連立方程式
\[Ax = b\]を解くために、ベクトル列\(\{x_n\}\)を
\[x_0 \to x_1 \to \dots \to x_N\]と反復的に
\[Ax_N = b\]となるように構成する。
この記事では、初期値\(x_0 = 0\)の場合を解説する。
定義（内積）
\[(x, y) := \sum x_i y_i\] \[(x, y)_A := (Ax, y)\]注意
\[J(x) := \frac{1}{2}(Ax, x) - (b, x)\]とおくと、
\[Ax = b \iff J(x) = \min_{y} J(y)\]証明
\(A\)が正定値対称なので、直交行列で対角化すれば、
\[A = P\Lambda P^{-1}\]となる。ただし、\(\Lambda\) は対角成分が\(\lambda_1, \ldots, \lambda_n (\lambda_i > 0)\)の対角行列である。いま、\(x = Py\)とおくと、
\[\begin{eqnarray} J(x) & = &\frac{1}{2}(APy, Py) - (b, Py) \\ &=& \frac{1}{2}(\Lambda y, y) - (c, y) \\ &=& \frac{1}{2}(\lambda_1 y_1^2 + \dots + \lambda_n y_n^2) - (c_1y_1 + \dots + c_ny_n)\end{eqnarray}\]（ただし、\(c := P^{-1}b\)とおいた）となる。したがって、\(J(x)\)が最小となるのは、\(y_i = c_i/\lambda_i\)となるときで、これは
\[y = \Lambda^{-1}c = \Lambda^{-1}P^{-1}b\]となるときであり、すなわち
\[Ax = APy = AP\Lambda^{-1}P^{-1}b = A A^{-1}b = b\]となるときである。
上の事実をつかって、連立方程式\(Ax = b\)を解くために、\(J(x) = \frac{1}{2}(Ax, x) - (b, x)\)を最小化するのを目標とする。
共役勾配法のアイデア
これから、共役勾配法のアイデアを説明する。
\(J(x)\)の等高線は、一般には（回転した）楕円（のn次元版）になる。なぜならば、\(A\)が正定値対称なので、\((Ax, x)\)は直交行列による基底変換で\(\lambda_1y_1^2 + \dots \lambda_n y_n^2\)の形になるからである。
\(J(x)\)の等高線が円（のn次元版）になるときを考えてみる。このときは、\(J(x)\)は各軸の方向に逐次最小化していけば、全体として最小化できる。どういうことかというと、たとえば\(J(x) = (x_1-2)^2 + (x_2-2)^2\)となっているとき、
\[J(0 + (t, 0))\]を最小化する\(t = t_0\)をまず求め、次に
\[J(0 + (t_0, 0) + (0, s))\]を最小化する\(s = s_0\) を求めると、
\(J(t_0, s_0)\)が\(J(x)\)の最小値になる。
一方で、\(J(x)\)の等高線が（回転した）楕円になっているときは、このように各軸の方向に1回ずつ最小化しても、全体として最小にならない。共役勾配法のアイデアは、\(A\)の内積\((x, y)_A\)での直交基底をとり、その基底の方向に逐次最小化していくことで、等高線が円の場合と同じように最小化するというものである。
詳細は齊藤宣一『数値解析』にゆずる。
参考文献
- 齊藤宣一『数値解析』
- 齊藤宣一『数値解析入門』