Mugichoko's blog

Mugichoko’s blog

プログラミングを中心としたメモ書き.

線形最小二乗 まとめ

基本事項

線形(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法」に詳しい.
  • 共役勾配法
  • その他の方法
    • [2]には「連立方程式の解の反復改良」として,丸め誤差の累積を解消し「完全な機械精度を回復するうまい手」が紹介されている.

単精度か倍精度か?

[2]によれば,変数の型を単精度(Float:32ビット浮動小数点)か倍精度(Double:64ビット浮動小数点)にするかを選ぶ基準は Nの大きさによるとのこと.

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

参考文献

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

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

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