4. 線形最小二乗法

4.1 概要

\(m\) 個の測定データ
\[
(x_i, y_i) = (x_i^* + e_{x_i}, y_i^* + e_{y_i}), i = 1, \dots, m
\] が与えられているものとします. ここでそれぞれのデータ \((x_i, y_i)\) は, 真値 \((x_i^*, y_i^*)\) に対して測定誤差 (\(e_{x_i} および e_{y_i}\)) を含んでいます.

各点の真値は, ある関数
\[
y_i^* = F(x_i^*)
\] を満たすものとします.

この関数 \(F(x)\) を \(n\) 個の基底関数 \(\phi_j\) の一次結合 \(f(x)\) により近似することにします.
\[
F(x) \simeq f(x) = c_1\phi_1(x) + c_2\phi_2(x) + \dots + c_n\phi_n(x)
\] この \(f(x)\) をモデル関数といいます. データ点が十分に多いとき (\(m \gg n\)), モデル関数がデータをできるだけうまくあてはめるようなパラメータ \(c_j (j = 1 \sim n)\) を求めることにします.

モデル関数の基底関数 \(\phi_j\) として最もよく使われるのは,
\[
\phi_j(x) = x^{j-1}
\] すなわち, 多項式近似です.
\[
f(x) = c_1 + c_2x + \dots + c_nx^{n-1}
\] 多項式近似以外でも一次結合で表すことができれば同じように扱うことができます. 例えば次のような場合です.

三角関数近似：
\[
f(x) = c_1sin(x) + c_2sin(2x) + \dots + c_nsin(nx)
\] 指数関数近似 (\(\lambda_j\) は定数とする)：
\[
f(x) = c_1exp(\lambda_1x) + c_2exp(\lambda_2x) + \dots + c_nexp(\lambda_nx)
\] さて, \(m\) 個の測定データをモデル関数に代入して縦に並べると次の連立一次方程式になります. これを観測方程式といいます.
\[
\boldsymbol{Ac} = \boldsymbol{y}
\] ただし, \(\boldsymbol{A}\) は \(m \times n\) 係数行列で, その要素は \(a_{ij} = \phi_j(x_i) (i = 1 \sim m, j = 1 \sim n)\) です. \(\boldsymbol{c}\) は \(n\) 解ベクトルで, パラメータ \(c_j (j = 1 \sim n)\) を要素とします. \(\boldsymbol{y}\) は \(m\) 右辺ベクトルで, 測定データ \(y_i (i = 1 \sim m)\) を要素とします.

この方程式は未知数の数よりも方程式の数の方が多い \((m > n)\) ので厳密に解くことはできない (不能な方程式) ので, 近似解を求めることにします. ひとつの方法として, ノルム \(||\boldsymbol{Ac} – \boldsymbol{y}||\) の二乗が最小になる解を求める方法が最小二乗法(*)です. このとき, 解 \(\boldsymbol{c}\) を最小二乗解といいます.

(*) すなわち, 残差 \(r_i = f(x_i) – y_i\) の二乗和 \(\sum_{i=1}^{m}r_i^2\) を最小にします.

4.2 線形最小二乗法の解法

最小二乗解は \(||\boldsymbol{Ac} – \boldsymbol{y}||^2\) の微分が 0 になる解として求めることができて, 次の方程式の解となります.
\[
\boldsymbol{A^TAc} = \boldsymbol{A^Ty}
\] これを正規方程式と呼び, \(未知数の数 = 方程式の数 = n\) となるので解くことができます. ただし, 正規方程式の係数行列の条件数は \(Cond(\boldsymbol{A^TA}) = (Cond(\boldsymbol{A}))^2\) と大きくなりやすく精度が悪くなることがあり, この解法は推奨されません.

代わって次の QR 分解が広く使われています.
\[
\boldsymbol{A} = \boldsymbol{QR}
\] ここで, \(\boldsymbol{Q}\) は \(m \times m\) 直交行列 (\(\boldsymbol{Q^TQ} = \boldsymbol{I}\)), \(\boldsymbol{R}\) は \(m \times n\) 上三角行列 (第 \(n+1 \sim m\) 行は 0) です.

QR 分解により元の連立一次方程式は次のように変形できます.
\[
\begin{align}
&\boldsymbol{Ac} = \boldsymbol{y} \\
&\boldsymbol{QRc} = \boldsymbol{y} \\
&\boldsymbol{Q^TQRc} = \boldsymbol{Q^Ty} \\
&\boldsymbol{Rc} = \boldsymbol{Q^Ty} \\
\end{align}
\] 最後の式において, 左辺の \(\boldsymbol{R}\) は上三角行列なので代入操作だけで容易に \(\boldsymbol{c}\) を求めることができます. こうして求めた \(\boldsymbol{c}\) は最小二乗解になっています.

パラメータ数 \(n\) はむやみに大きくすればよいというものではありません. 例えば, 明らかに直線にのっているデータを 2 次以上の多項式で近似しても意味がありません. \(n\) が大きすぎると係数行列 \(\boldsymbol{A}\) の列が一次従属になってしまい, QR 分解が失敗します. これをランク落ちといいます (ランク数は一次独立な列の数です). このような場合にはモデル関数の見直しが必要です.

4.3 XLPack を使った線形最小二乗法の解き方

VBA サブルーチン Dgels あるいはワークシート関数 WDgels を使うことにより精度よく最小二乗解を求めることができます.

これらは, 線形計算ライブラリ LAPACK のサブルーチン DGELS を使用しており, QR 分解により最小二乗解を求めます. ただし, ランク落ちがあってはいけません.

例題: 多項式近似

次のデータは アメリカ国立標準技術研究所 (NIST) のホームページ (http://www.itl.nist.gov/div898/strd/lls/data/Norris.shtml) に掲載されているテストデータです.

 
     y          x
    0.1        0.2
  338.8      337.4
  118.1      118.2
  888.0      884.6
    9.2       10.1
  228.1      226.5
  668.5      666.3
  998.5      996.3
  449.1      448.6
  778.9      777.0
  559.2      558.2
    0.3        0.4
    0.1        0.6
  778.1      775.5
  668.8      666.9
  339.3      338.0
  448.9      447.5
   10.8       11.6
  557.7      556.0
  228.3      228.1
  998.0      995.8
  888.8      887.6
  119.6      120.2
    0.3        0.3
    0.6        0.3
  557.6      556.8
  339.3      339.1
  888.0      887.2
  998.5      999.0
  778.9      779.0
   10.2       11.1
  117.6      118.3
  228.9      229.2
  668.4      669.1
  449.2      448.9
    0.2        0.5

 
これを 1 次多項式(直線)
\[
f(x) = c_1 + c_2x
\] により近似します.

4.3.1 ワークシート関数を使用した解き方

(1) ワークシートの適当な場所にデータを入力します(オレンジ色のセル).

(2) 係数行列 \(\boldsymbol{A} \space (a_{ij} = \phi_j(x_i) (i = 1 \sim m, j = 1 \sim n))\) を作ります. この場合は, \(a_{i1} = 1, a_{i2} = x_i (i = 1 \sim m)\) となります.

(3) 求めるパラメータ \(\boldsymbol{c}\) およびその他の結果を入れる \(n + 2\) 個のセルを選択してワークシート関数 WDgels を設定します (緑色のセル).

WDgels の必要なパラメータは, Trans, Cov, M, N, A, B, Nrhs です.
Trans: “N”: \(\boldsymbol{Ax} = \boldsymbol{b}\) を解く. “T”: \(\boldsymbol{A^Tx} = \boldsymbol{b}\) を解く.
Cov: “N”: 分散共分散行列を計算しない. “D”: 分散 (分散共分散行列の対角要素) を計算する. “C”: 分散共分散行列を計算する.
M: データ数.
N: パラメータ数 (この例では 2).
A: 行列 \(\boldsymbol{A}\) のセル範囲.
B: 右辺ベクトル (= ベクトル \(\boldsymbol{y}\)) のセル範囲.
Nrhs: 右辺ベクトルの本数 (\(\boldsymbol{y}\) だけ変えて複数回計算する場合に入力. 省略時は 1 とみなす).

入力が終了したら Ctrl + Shift + Enter を押すと結果が得られます. ここでは, データを丸で, 計算値を直線で表すグラフも作成してみました.

2つのパラメータが求められました (\(c_1 = -0.26232, c_2 = 1.002117\)). これらは, NIST の値とぴったり一致しています.

4.3.2 VBAプログラムを使用した解き方

上と同じ例題を VBA プログラムにより解きます. Dgels を使ったプログラム例を示します.

Sub Start()
    Const MMax = 100, NMax = 5
    Dim M As Long, N As Long
    Dim A(MMax, NMax) As Double, B(MMax) As Double
    Dim Info As Long, I As Long, J As Long
    '--- Input data
    M = 36: N = 2
    For I = 0 To M - 1
        For J = 0 To N - 1
            A(I, J) = Cells(5 + I, 3 + J)
        Next
        B(I) = Cells(5 + I, 2)
    Next
    '--- Compute least squares solution
    Call Dgels("N", M, N, A(), B(), Info)
    '--- Output result
    For J = 0 To N - 1
        Cells(5 + J, 5) = B(J)
    Next
    Cells(7, 5) = Info
End Sub

所定の位置にデータを入力し, マクロ Start を実行すると次のように上と同じ結果が得られます. ワークシート関数を使用したときと異なり, 係数行列 \(\boldsymbol{A}\) と右辺ベクトル \(\boldsymbol{b}\) の値を入力しただけでは結果が得られず, VBA プログラムを実行してやる必要があります.