14. 非線形最小二乗法

14.1 概要

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

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

\(n\) ベクトル \(\boldsymbol{c}\) の要素は \(n\) 個のパラメータ \(c_1, c_2, \dots, c_n\) であるとします. 関数 \(F(x)\) をこの (\(n\) 個の) パラメータ \(\boldsymbol{c}\) を含む一般関数 \(f(x; \boldsymbol{c})\) (モデル関数という) で近似することにします.
\[
F(x) \simeq f(x; \boldsymbol{c}) = f(x; c_1, c_2, \dots, c_n)
\] このときモデル関数 \(f(x; \boldsymbol{c})\) がデータをできるだけうまくあてはめるパラメータ \(\boldsymbol{c}\) を求める方法のひとつが非線形最小二乗法で, 残差
\[
r_i(\boldsymbol{c}) = f(x_i; \boldsymbol{c}) – y_i
\] の二乗和 \(\sum_{i=1}^m r_i(\boldsymbol{x})^2\) が最小になるようなパラメータ \(\boldsymbol{c}\) を求めます.

14.2 非線形最小二乗法の解法

\(n\) ベクトル \(\boldsymbol{x}\) について, 次のような関数 \(R(\boldsymbol{x})\) を考えます.
\[
R(\boldsymbol{x}) = (1/2)\boldsymbol{r}(\boldsymbol{x})^T\boldsymbol{r}(\boldsymbol{x}) = \sum_{i=1}^m r_i(\boldsymbol{x})^2
\] ただし,
\[
\boldsymbol{r}(\boldsymbol{x}) =
\begin{pmatrix}
r_1(\boldsymbol{x}) \\
r_2(\boldsymbol{x}) \\
\vdots \\
r_m(\boldsymbol{x}) \\
\end{pmatrix}
\] すなわち, \(R(\boldsymbol{x})\) は \(r_i(\boldsymbol{x})\) の二乗和を表す関数とします.

このとき, \(R(\boldsymbol{x})\) の勾配 \(\nabla R(\boldsymbol{x})\) とヘッセ行列 \(\boldsymbol{H}(\boldsymbol{x})\) は次のように表すことができます.
\[
\nabla R(\boldsymbol{x}) = \boldsymbol{J}(\boldsymbol{x})^T\boldsymbol{r}(\boldsymbol{x}) \\
\boldsymbol{H}(\boldsymbol{x}) = \boldsymbol{J}(\boldsymbol{x})^T\boldsymbol{J}(\boldsymbol{x}) + \boldsymbol{S}(\boldsymbol{x}) \\
\] ただし, \(\boldsymbol{J}(\boldsymbol{x})\) はヤコビ行列, \(\boldsymbol{S}(\boldsymbol{x}) は \sum_{i=1}^m \partial^2r_i(\boldsymbol{x})/\partial x_j\partial x_k \space r_i(\boldsymbol{x}) \space (j, k = 1 \sim n)\) を表します.

すなわち, \(R(\boldsymbol{x})\) を目的関数として非線形最適化のニュートン法を適用すると, \(R(\boldsymbol{x})\) が \(r_i(\boldsymbol{x})\) の二乗和であるという特別な構造を持っていることにより, 反復式で必要な勾配とヘッセ行列はヤコビ行列を使って求めることができることになります.

ここで, 14.1 の説明における \(f(x; \boldsymbol{c})\) と \(y\) を使って改めて \(r_i\) を次のように表します.
\[
r_i(\boldsymbol{c}) = f(x_i; \boldsymbol{c}) – y_i \space (i = 1 \sim m)
\] すなわち, \(r_i\) はパラメータ \(\boldsymbol{c}\) の関数で, i 番目の測定データにおけるモデル関数値の残差を表します. このようにすると, 残差二乗和 \(R(\boldsymbol{c})\) を目的関数とした非線形最適化問題を解いて \(\boldsymbol{c}\) の非線形最小二乗解を得ることができます. その際には, 反復で必要な勾配とヘッセ行列はヤコビ行列を使って求めることができます.

以下では, 再び \(\boldsymbol{c}\) を \(\boldsymbol{x}\) で表して \(R(\boldsymbol{x})\) と表記することにします.

14.2.1 ガウス・ニュートン法

ニュートン法を使って \(R(\boldsymbol{x})\) を目的関数とした非線形最適化問題を解きます. ただし, ヘッセ行列を求めるときに \(\boldsymbol{S}(\boldsymbol{x})\) の項を省略して (\(= 0\) として) 近似します.

そうすると, 方向ベクトル \(\boldsymbol{d_k}\) は次の連立一次方程式を解いて求めることができます.
\[
\boldsymbol{J}(\boldsymbol{x_k})^T\boldsymbol{J}(\boldsymbol{x_k})\boldsymbol{d_k} = -\boldsymbol{J}(\boldsymbol{x_k})^T\boldsymbol{r}(\boldsymbol{x_k})
\] この式は線形最小二乗法の正規方程式の形をしているので, このまま計算せずに線形最小二乗法と同じように QR 分解などを用いて求めます. (4 章を参照.)

14.2.2 レーベンバーグ・マルカート法

ガウス・ニュートン法では反復の途中でヤコビ行列がランク落ちしてしまうことがあり, 連立一次方程式を解くことができなくなるという問題があります.

それを解決するために, パラメータ \(\lambda_k\) を導入し \(\lambda_k\boldsymbol{I}\) の項を追加して, 係数行列が正定値になるようにした方法がレーベンバーグ・マルカート法です.
\[
(\boldsymbol{J}(\boldsymbol{x_k})^T \boldsymbol{J}(\boldsymbol{x_k}) + \lambda_k \boldsymbol{I})\boldsymbol{d_k} = -\boldsymbol{J}(\boldsymbol{x_k})^T\boldsymbol{r}(\boldsymbol{x_k})
\] パラメータ \(\lambda_k = 0\) であればガウス・ニュートン法そのものになります. また, \(\lambda_k\) の値が大きくなると最急降下法に近づいていきます.

パラメータの値により収束性が異なり, その決め方にはいくつかの方法が提案されています.

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

レーベンバーグ・マルカート法の VBA サブルーチン Lmdif1 を使うことができます. これは微分 (ヤコビ行列) を差分近似で計算し, ユーザーによる微分計算を不要にしたサブルーチンです. Lmdif1 は XLPack ソルバーから使うこともできます.

非線形最適化の場合と同じように, 目的の解にできるだけ近い初期値を与えるようにすべきであることは変わりません.

例題

NIST (アメリカ国立標準技術研究所) のホームページ (https://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml) に掲載されているテストデータを使用します.

次のようなモデル関数を使用します.
\[
f(x) = c_1(1 – exp(-c_2x))
\] データは次のとおりです.
\[
\begin{array}{cccc}
\hline
y & x \\ \hline
10.07 & 77.6 \\
14.73 & 114.9 \\
17.94 & 141.1 \\
23.93 & 190.8 \\
29.61 & 239.9 \\
35.18 & 289.0 \\
40.02 & 332.8 \\
44.82 & 378.4 \\
50.76 & 434.8 \\
55.05 & 477.3 \\
61.01 & 536.8 \\
66.40 & 593.1 \\
75.47 & 689.1 \\
81.78 & 760.0 \\
\hline
\end{array}
\] このデータをプロットすると次のようになります.

テストデータの初期値は 2 種類指定されており, (500, 0.0001) と (250, 0.0005) です. 示されている解は \(c_1 = 2.3894212918 \times 10^2, c_2 = 5.5015643181 \times 10^{-4}\) です.

14.3.1 VBA プログラムを使用した解き方 (1)

Lmdif1 を使ったプログラム例を示します.

Const Ndata = 14, Nparam = 2
Dim Xdata(Ndata - 1) As Double, Ydata(Ndata - 1) As Double

Sub F(M As Long, N As Long, X() As Double, Fvec() As Double, IFlag As Long)
    Dim I As Long
    For I = 0 To M - 1
        Fvec(I) = Ydata(I) - X(0) * (1 - Exp(-Xdata(I) * X(1)))
    Next
End Sub

Sub Start()
    Dim X(Nparam - 1) As Double, Fvec(Ndata - 1) As Double, Tol As Double
    Dim Info As Long, I As Long
    '--- Input data
    For I = 0 To Ndata - 1
        Xdata(I) = Cells(9 + I, 1)
        Ydata(I) = Cells(9 + I, 2)
    Next
    '--- Start calculation
    Tol = 0.0000000001  '** 1e-10
    For I = 0 To 1
        '--- Initial values
        X(0) = Cells(6 + I, 1)
        X(1) = Cells(6 + I, 2)
        '--- Nonlinear least squares
        Call Lmdif1(AddressOf F, Ndata, Nparam, X(), Fvec(), Tol, Info)
        '--- Output result
        Cells(6 + I, 3) = X(0)
        Cells(6 + I, 4) = X(1)
        Cells(6 + I, 5) = Info
    Next
End Sub

目的関数を外部サブルーチンとして用意し, 初期値と要求精度を指定して Lmdif1 を呼び出します. Lmdif1 は, 初期値として与えられた近似解を出発点に反復計算を行い, その近傍の最小二乗解を求めます.

このプログラムを実行すると次の結果が得られました. なお, 得られたパラメータを使った計算値(オレンジの線)とデータ(青い点)をグラフにしました.

この例ではどちらの初期値の場合も正しい解が求められました.

14.3.2 VBA プログラムを使用した解き方 (2)

リバースコミュニケーション版 (RCI) の Lmdif1_r を使ったプログラム例を示します.

Sub Start_r()
    Const M As Long = 14, N As Long = 2
    Dim X(Nparam - 1) As Double, Fvec(Ndata - 1) As Double, Tol As Double
    Dim Info As Long, I As Long, J As Long
    Dim XX(Nparam - 1) As Double, YY(Ndata - 1) As Double, IRev As Long, IFlag As Long
    Dim Xdata(Ndata - 1) As Double, Ydata(Ndata - 1) As Double
    '--- Input data
    For I = 0 To Ndata - 1
        Xdata(I) = Cells(9 + I, 1)
        Ydata(I) = Cells(9 + I, 2)
    Next
    '--- Start calculation
    Tol = 0.0000000001  '** 1e-10
    For I = 0 To 1
        '--- Initial values
        X(0) = Cells(6 + I, 1)
        X(1) = Cells(6 + I, 2)
        '--- Nonlinear least squares
        IRev = 0
        Do
            Call Lmdif1_r(Ndata, Nparam, X(), Fvec(), Tol, Info, XX(), YY(), IRev)
            If IRev = 1 Or IRev = 2 Then
                For J = 0 To Ndata - 1
                    YY(J) = Ydata(J) - XX(0) * (1 - Exp(-Xdata(J) * XX(1)))
                Next
            End If
        Loop While IRev <> 0
        '--- Output result
        Cells(6 + I, 3) = X(0)
        Cells(6 + I, 4) = X(1)
        Cells(6 + I, 5) = Info
    Next
End Sub

RCI サブルーチンでは, 目的関数を外部関数として与えるのではなく, IRev = 1 または 2 のときに XX() の値を使って関数値を計算し YY() に入れて再度 Lmdif1_r を呼び出します. RCIの詳細については こちら を参照してください.

このプログラムを実行すると上と同じ結果が得られます.

14.3.3 ソルバーを使用した解き方

XLPackソルバーアドインの「非線形最小二乗法」を使って解くこともできます. C 列で B8, C8 のパラメータ値を使ってモデル関数値を計算します. 例えば, 式 =B$8*(1-EXP(-C$8*A11)) がC11セルに入力されています. D 列では残差を求めています. 例えば, 式 =B11-C11 がD11セルに入力されています. ソルバーはこの残差から最小二乗解を求めます.

ソルバーについては こちら も参照ください.