基本事項

線形（1次）方程式の組（連立1次方程式）は以下のように表される．

$$ \begin{split} a_{11} x_1 + a_{12} x_2 + ... + a_{1N} x_N = b_1 \\ a_{21} x_1 + a_{22} x_2 + ... + a_{2N} x_N = b_2 \\ \vdots \qquad\qquad\qquad \\ a_{M1} x_1 + a_{M2} x_2 + ... + a_{MN} x_N = b_M. \\ \end{split} \tag{1} $$

ただし， $M$ は式の数， $N$ は未知数の数である．また， $a_{ij}$ と $b_i$ が既知の値 $(1 \leq i \leq M, 1 \leq j \leq N)$ ， $x_{j}$ を未知数とする．ここで考える問題は，どうやって $x_{j}$ を求めるのか，である．そのために，まず，式1の行列表現を考える．

$$ {\bf A x} = {\bf b} \tag{2} $$

ただし， ${\bf A}$ は $M \times N$ 行列の係数行列， ${\bf b}$ は $M$ 個の要素を持つ列ベクトルである．すなわち，式2は以下のように表せる．

$$ \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1N} \\ a_{21} & a_{22} & \cdots & a_{2N} \\ \vdots & & \ddots & \vdots \\ a_{M1} & a_{M2} & \cdots & a_{MN} \end{bmatrix} \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{N} \\ \end{bmatrix} = \begin{bmatrix} b_{1} \\ b_{2} \\ \vdots \\ b_{M} \\ \end{bmatrix} $$

この ${\bf x}$ （つまり， $x_{j}$ ）を求める解法には「直接解法」と「反復解法」の2つの方法が存在する[1]．

直接解法

「式を直接変形して解を求める[1]」のが直接解法である．具体的な解法には以下の方法がある．

Gauss-Jordan法 (Gauss-Jordan Elimination)
- 正方行列 $(M=N)$ に適用可能．「すべての右辺を同時に格納・操作しなければならない，逆行列がふようならこれより3倍速い方法がある[2]」という短所の一方で，直接解法の中では「最も安定な部類に属し，完全ピボット選択を用いれば他の方法より若干安定である[2]」という長所を持つ．尚，3倍速い方法とは下記のLU分解のことである．
- [2]の「2.1 Gauss-Jordan法」に詳しい．
Gauss消去法 (Gaussian Elimination)
- 正方行列 $(M=N)$ に適用可能．前進消去と後退代入を行い，解 ${\bf x}$ を求める．「ピボット選択なしのGauss法は絶対にしてはならない」旨が[2]のp.51に述べられている．
- [1]の「3.1 ガウス消去法」「3.2 部分ピボット選択付きガウス消去法」に詳しい．[2]の「2.2 Gaussの消去法と後退代入」には詳しいアルゴリズムは述べられていない．これは，この方法とLU分解法とがほとんど同じだからである．
LU分解 (LU Decomposition)
- 正方行列 $(M=N)$ に適用可能． ${\bf A}$ を下三角行列 (Lower triangular matrix) と上三角行列 (Upper triangular matrix) に分解することで解 ${\bf x}$ を求める方法，「このLU分解は，部分ピボット選択付きガウス消去法とほとんど同じです．唯一の違いは，LU分解が右辺ベクトル ${\bf b}$ に依存していないという点です．[1]」．つまり，「いったん ${\bf A}$ をLU分解してしまえば，新たな右辺について方程式を解く必要が生じても，LU分解をやり直す必要がない．[2]」という点で，LU分解は Gauss-Jordan法やGauss消去法より優れている．
- [1]の「3.3 LU分解」及び[2]の「2.3 LU分解」に詳しい．
特異値分解 (Singular Value Decomposition; SVD)
- $M×N$ 行列に適用可能． ${\bf A}$ を $M×N$ の列直行行列U， $N×N$ の対角行列 ${\bf W}$ （対角成分は非負）， $N×N$ の直行行列 ${\bf V}$ の転置 ${\bf V^{\rm T}}$ に分解して解 ${\bf x}$ を求める方法．基本的に， ${\bf x} = {\bf V}\lbrack \mathrm{diag}(1 / w_{j}) \rbrack ({\bf U}^{\rm T} {\bf b})$ で求まる．
- 「SVDは，ほとんどの線形最小2乗 (linear least squares) 問題をとくための，極めつけの方法でもある．[2]」とあるように，線形最小2乗の解法としても有名．
- [2]の「2.9 特異値分解」や[3]の「4.2.3 特異値分解と一般逆行列」に詳しい．
その他の方法
- 以下に列挙するのは特別な場合に用いられる方法である．「コレスキー分解法」 ${\bf A}$ が正定値行列の場合に用いる．[1]の「3.5 コレスキー分解法」に詳しい．3重対角な連立一次方程式の場合は[2]の「2.6 3重対角な連立1次方程式」に詳しい．Vandermonde行列とToeplitz行列の場合は[2]の「2.8 Vandermonde行列とToeplitz行列」に詳しい．疎な連立1次方程式の場合は[2]の「2.10 疎な連立1次方程式」に詳しい．

反復解法

「ある初期値を与えて，それをある規則に従って反復改良して解を求める[1]」のが反復解法である．具体的な解法には以下の方法がある．今のところ，あまり実例を見たことがない．

ヤコビ法 (Jacobi Method)
- [1]の「5.2 ヤコビ法」に詳しい．
ガウス・ザイデル法 (Gauss-Seidel Method)
- [1]の「5.3 ガウス・ザイデル法」に詳しい．
SOR法 (Successive Over-Relaxation Method)
- [1]の「5.4 SOR法」に詳しい．
共役勾配法
- [1]の「5.5 共役勾配法」に詳しい．
その他の方法
- [2]には「連立方程式の解の反復改良」として，丸め誤差の累積を解消し「完全な機械精度を回復するうまい手」が紹介されている．

単精度か倍精度か？

[2]によれば，変数の型を単精度（Float：32ビット浮動小数点）か倍精度（Double：64ビット浮動小数点）にするかを選ぶ基準は $N$ の大きさによるとのこと．

が20や
- 単精度
- 特異に近くないなら複雑な方法に頼らなくても機械的に解ける
- $N$ が5程度ならまず問題ない
- 特異な連立方程式の場合， $N$ が10でも複雑な方法に頼らなければならないことがある
- 特異値分解により，特異な問題を非特異な問題に変えられることがある
が数百
- 倍精度
- この辺になると，制度より時間が制約要因となる
が数千
- 倍精度
- 係数が疎（つまり，ほとんどの係数がゼロ）であれば，この条件を利用して解ける

参考文献

[1] 皆本晃弥著，C言語による数値計算入門〜解法・アルゴリズム・プログラム〜，サイエンス社，2007．

[2] William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery 著，丹慶勝市，奥村晴彦，佐藤俊郎，小林誠訳，NUMERICAL RECIPES in C [日本語版] ニューメリカルレシピ・イン・シー〜C言語による数値計算のレシピ〜，技術評論社，2012．

[3] 金谷健一著，これなら分かる最適化学〜基礎原理から計算手法まで〜，共立出版，2005．

Mugichoko’s blog

しがない研究者のプログラミングを中心としたメモ書き．

線形最小二乗まとめ

基本事項

直接解法

反復解法

単精度か倍精度か？

参考文献