非線形最小二乗法

観測値の信頼性が等しくない場合、加重二乗和を最小化することができます。 $${\displaystyle S=\sum _{i=1}^{m}W_{ii}r_{i}^{2}.}$$

対角重み行列Wの各要素は、理想的には、測定の誤差分散の逆数に等しくなければならない。 ^{[ 1 ]} より一般的には、正規方程式は次のようになる。 $${\displaystyle \left(\mathbf {J} ^{\mathsf {T}}\mathbf {WJ} \right)\Delta {\boldsymbol {\beta }}=\mathbf {J} ^{\mathsf {T}}\mathbf {W} \ \Delta \mathbf {y} .}$$

条件不良や発散の問題は、最適値に近い初期パラメータ推定値を見つけることで修正できます。これを行う良い方法は、コンピュータシミュレーションです。観測データと計算データの両方が画面に表示されます。観測データと計算データの一致が十分に良好になるまで、モデルのパラメータを手動で調整します。これは主観的な判断になりますが、非線形改良の良い出発点を見つけるには十分です。初期パラメータ推定値は、変換または線形化を使用して作成できます。さらに良いのは、確率的漏斗アルゴリズムなどの進化アルゴリズムを使用すると、最適なパラメータ推定値を取り囲む凸状の吸引領域が得られます。 ランダム化とエリート主義を使用し、ニュートン法を併用するハイブリッドアルゴリズムは、有用で計算効率が高いことが示されています。

以下に説明する方法のいずれも、解決策を見つけるために適用できます。

収束の常識的な基準は、二乗和が反復ごとに増加しないことである。しかし、この基準は、さまざまな理由から、実際には実装が難しい場合が多い。有用な収束基準は、 $${\displaystyle \left|{\frac {S^{k}-S^{k+1}}{S^{k}}}\right|<0.0001.}$$ 0.0001という値はやや恣意的であり、変更する必要があるかもしれない。特に、実験誤差が大きい場合は、値を大きくする必要があるかもしれない。代替基準としては、 $${\displaystyle \left|{\frac {\Delta \beta _{j}}{\beta _{j}}}\right|<0.001,\qquad j=1,\dots ,n.}$$

繰り返しますが、数値はやや恣意的です。0.001は、各パラメータを0.1%の精度で精緻化することを指定するのと同等です。これは、パラメータの最大相対標準偏差よりも小さい場合に妥当です。

ヤコビ行列の要素の解析的表現を導出することが非常に困難、あるいは不可能なモデルも存在する。その場合、数値近似 $${\displaystyle {\frac {\partial f(x_{i},{\boldsymbol {\beta }})}{\partial \beta _{j}}}\approx {\frac {\delta f(x_{i},{\boldsymbol {\beta }})}{\delta \beta _{j}}}}$$ は、計算によって得られる。${\displaystyle f(x_{i},{\boldsymbol {\beta }})}$のために${\displaystyle \beta _{j}}$そして${\displaystyle \beta _{j}+\delta \beta _{j}}$増加分、${\displaystyle \delta \beta _{j}}$数値微分が大きすぎると近似誤差を生じ、小さすぎると丸め誤差を生じないように、サイズを選択する必要があります。

重み付き最小二乗法のページにある該当セクションに、関連情報が記載されています。

複数の極小値が発生する可能性のある状況は様々であり、その例としては以下のようなものがある。

• パラメータを2乗以上に累乗します。たとえば、データをローレンツ曲線にフィッティングする場合などです。$${\displaystyle f(x_{i},{\boldsymbol {\beta }})={\frac {\alpha }{1+\left({\frac {\gamma -x_{i}}{\beta }}\right)^{2}}}}$$どこ${\displaystyle \alpha }$高さは、${\displaystyle \gamma }$位置と${\displaystyle \beta }$半分の高さにおける半分の幅は、半分の幅には 2 つの解があります。${\displaystyle {\hat {\beta }}}$そして${\displaystyle -{\hat {\beta }}}$これらは目的関数に対して同じ最適値を与える。
• 2 つのパラメータは、モデルの値を変更せずに交換できます。簡単な例として、モデルに 2 つのパラメータの積が含まれている場合が挙げられます。${\displaystyle \alpha \beta }$同じ値になります${\displaystyle \beta \alpha }$。
• パラメータは、例えば三角関数の中にあります。${\displaystyle \sin \beta }$同じ値を持つ${\displaystyle {\hat {\beta }}+2n\pi }$例として、レーベンバーグ・マルカート法を参照してください。

複数の最小値が存在する場合、目的関数の値が必ずしも等しいとは限りません。偽の最小値（局所的最小値とも呼ばれる）は、目的関数の値が、いわゆる大域的最小値における値よりも大きい場合に発生します。見つかった最小値が大域的最小値であることを確実にするには、パラメータの初期値を大きく変えて最適化を開始する必要があります。開始点に関わらず同じ最小値が見つかった場合、それは大域的最小値である可能性が高いと言えます。

複数の最小値が存在する場合、重要な結果が生じます。目的関数は、2 つの最小値の間のどこかに停留点 (例えば、最大値または鞍点) を持ちます。目的関数の停留点では、勾配がゼロになり、降下方向が一意に存在しないため、正規方程式行列は正定値ではありません。停留点に近い点 (パラメータ値の集合) からの改良は条件が悪く、開始点としては避けるべきです。たとえば、ローレンツ関数をフィッティングする場合、ローレンツ関数の半幅がゼロのとき、正規方程式行列は正定値ではありません。^{[ 2 ]}

非線形モデルは、場合によっては線形モデルに変換できます。このような近似は、例えば、最良推定量の近傍でよく適用され、ほとんどの反復最小化アルゴリズムにおける基本的な仮定の1つです。線形近似が有効な場合、モデルは一般化最小二乗法による推論に直接使用でき、線形テンプレート適合^{[ 3 ]}の式が適用されます。

線形近似のもう1つの例は、モデルが単純な指数関数である場合です。 $${\displaystyle f(x_{i},{\boldsymbol {\beta }})=\alpha e^{\beta x_{i}},}$$ これは対数を取ることで線形モデルに変換できる。 $${\displaystyle \log f(x_{i},{\boldsymbol {\beta }})=\log \alpha +\beta x_{i}}$$グラフ上では、これは片対数グラフ上 での作業に相当します。二乗和は次のようになります。 $${\displaystyle S=\sum _{i}(\log y_{i}-\log \alpha -\beta x_{i})^{2}.}$$ 誤差が乗法的で対数正規分布に従う場合を除き、この手順は避けるべきです。なぜなら、誤解を招く結果が生じる可能性があるからです。これは、 yの実験誤差がどのような値であっても、 log yの誤差は異なるためです。したがって、変換された二乗和を最小化すると、パラメータ値とその計算された標準偏差の両方で異なる結果が得られます。しかし、乗法的誤差が対数正規分布に従う場合、この手順によって偏りのない一貫性のあるパラメータ推定値が得られます。

別の例として、2つのパラメータを決定するために使用されるミカエリス・メンテン速度論が挙げられる。${\displaystyle V_{\max }}$そして${\displaystyle K_{m}}$: $${\displaystyle v={\frac {V_{\max }[S]}{K_{m}+[S]}}.}$$ ラインウィーバー・バークプロット$${\displaystyle {\frac {1}{v}}={\frac {1}{V_{\max }}}+{\frac {K_{m}}{V_{\max }[S]}}}$$ の${\textstyle {\frac {1}{v}}}$に対して${\textstyle {\frac {1}{[S]}}}$パラメータに関して線形である${\textstyle {\frac {1}{V_{\max }}}}$そして${\textstyle {\frac {K_{m}}{V_{\max }}}}$しかし、データエラーに非常に敏感で、独立変数の特定の範囲にデータを適合させる傾向が強い。${\displaystyle [S]}$。

正規方程式 $${\displaystyle \left(\mathbf {J} ^{\mathsf {T}}\mathbf {WJ} \right)\Delta {\boldsymbol {\beta }}=\left(\mathbf {J} ^{\mathsf {T}}\mathbf {W} \right)\Delta \mathbf {y} }$$ 解決できるかもしれない${\displaystyle \Delta {\boldsymbol {\beta }}}$線形最小二乗法で説明されているように、コレスキー分解によってパラメータが反復的に更新されます。 $${\displaystyle {\boldsymbol {\beta }}^{k+1}={\boldsymbol {\beta }}^{k}+\Delta {\boldsymbol {\beta }}}$$ ここで、kは反復回数である。この方法は単純なモデルには十分かもしれないが、発散が発生すると失敗する。したがって、発散を防ぐ対策が不可欠である。

二乗和の最小値は、正規方程式を立てる必要のない方法で求めることができます。線形化モデルを用いた残差は次のように表すことができます。 $${\displaystyle \mathbf {r} =\Delta \mathbf {y} -\mathbf {J} \,\Delta {\boldsymbol {\beta }}.}$$ ヤコビ行列は直交分解の対象となり、QR分解はその過程を説明するために用いられる。 $${\displaystyle \mathbf {J} =\mathbf {QR} }$$ ここでQは直交行列である${\displaystyle m\times m}$行列であり、 Rは${\displaystyle m\times n}$分割された行列${\displaystyle n\times n}$ブロック、${\displaystyle \mathbf {R} _{n}}$、そして${\displaystyle (m-n)\times n}$ゼロブロック。${\displaystyle \mathbf {R} _{n}}$上三角形です。

$${\displaystyle \mathbf {R} ={\begin{bmatrix}\mathbf {R} _{n}\\\mathbf {0} \end{bmatrix}}}$$

残差ベクトルは左から乗算される${\displaystyle \mathbf {Q} ^{\mathsf {T}}}$。

$${\displaystyle \mathbf {Q} ^{\mathsf {T}}\mathbf {r} =\mathbf {Q} ^{\mathsf {T}}\ \Delta \mathbf {y} -\mathbf {R} \ \Delta {\boldsymbol {\beta }}={\begin{bmatrix}\left(\mathbf {Q} ^{\mathsf {T}}\ \Delta \mathbf {y} -\mathbf {R} \ \Delta {\boldsymbol {\beta }}\right)_{n}\\\left(\mathbf {Q} ^{\mathsf {T}}\ \Delta \mathbf {y} \right)_{m-n}\end{bmatrix}}}$$

これは二乗和には影響しない。${\displaystyle S=\mathbf {r} ^{\mathsf {T}}\mathbf {Q} \mathbf {Q} ^{\mathsf {T}}\mathbf {r} =\mathbf {r} ^{\mathsf {T}}\mathbf {r} }$Qは直交しているため、 Sの最小値は上側のブロックがゼロのときに得られます。したがって、シフトベクトルは、次の式を解くことによって求められます。 $${\displaystyle \mathbf {R} _{n}\ \Delta {\boldsymbol {\beta }}=\left(\mathbf {Q} ^{\mathsf {T}}\ \Delta \mathbf {y} \right)_{n}.}$$

Rは上三角関数であるため、これらの式は容易に解くことができる。

直交分解法の変形として、特異値分解があり、これはRをさらに直交変換によって対角化するものである。

$${\displaystyle \mathbf {J} =\mathbf {U} {\boldsymbol {\Sigma }}\mathbf {V} ^{\mathsf {T}}}$$ どこ${\displaystyle \mathbf {U} }$直交しており、${\displaystyle {\boldsymbol {\Sigma }}}$は特異値の対角行列であり、${\displaystyle \mathbf {V} }$は、の固有ベクトルの直交行列である。${\displaystyle \mathbf {J} ^{\mathsf {T}}\mathbf {J} }$または同等に右特異ベクトル${\displaystyle \mathbf {J} }$この場合、シフトベクトルは次のように与えられる。 $${\displaystyle \Delta {\boldsymbol {\beta }}=\mathbf {V} {\boldsymbol {\Sigma }}^{-1}\left(\mathbf {U} ^{\mathsf {T}}\ \Delta \mathbf {y} \right)_{n}.}$$

この式の比較的単純さは、非線形最小二乗法の理論的解析において非常に有用である。特異値分解の応用については、LawsonとHanson ^{[ 6 ]で詳しく議論されている。}

There are many examples in the scientific literature where different methods have been used for non-linear data-fitting problems.

• Inclusion of second derivatives in The Taylor series expansion of the model function. This is Newton's method in optimization. $${\displaystyle f(x_{i},{\boldsymbol {\beta }})=f^{k}(x_{i},{\boldsymbol {\beta }})+\sum _{j}J_{ij}\,\Delta \beta _{j}+{\frac {1}{2}}\sum _{j}\sum _{k}\Delta \beta _{j}\,\Delta \beta _{k}\,H_{jk_{(i)}},\ H_{jk_{(i)}}={\frac {\partial ^{2}f(x_{i},{\boldsymbol {\beta }})}{\partial \beta _{j}\,\partial \beta _{k}}}.}$$ The matrix H is known as the Hessian matrix. Although this model has better convergence properties near to the minimum, it is much worse when the parameters are far from their optimal values. Calculation of the Hessian adds to the complexity of the algorithm. This method is not in general use.
• Davidon–Fletcher–Powell method. This method, a form of pseudo-Newton method, is similar to the one above but calculates the Hessian by successive approximation, to avoid having to use analytical expressions for the second derivatives.
• Steepest descent. Although a reduction in the sum of squares is guaranteed when the shift vector points in the direction of steepest descent, this method often performs poorly. When the parameter values are far from optimal the direction of the steepest descent vector, which is normal (perpendicular) to the contours of the objective function, is very different from the direction of the Gauss–Newton vector. This makes divergence much more likely, especially as the minimum along the direction of steepest descent may correspond to a small fraction of the length of the steepest descent vector. When the contours of the objective function are very eccentric, due to there being high correlation between parameters, the steepest descent iterations, with shift-cutting, follow a slow, zig-zag trajectory towards the minimum.
• Conjugate gradient search. This is an improved steepest descent based method with good theoretical convergence properties, although it can fail on finite-precision digital computers even when used on quadratic problems.^[8]

Direct search methods depend on evaluations of the objective function at a variety of parameter values and do not use derivatives at all. They offer alternatives to the use of numerical derivatives in the Gauss–Newton method and gradient methods.

• Alternating variable search.^[4] Each parameter is varied in turn by adding a fixed or variable increment to it and retaining the value that brings about a reduction in the sum of squares. The method is simple and effective when the parameters are not highly correlated. It has very poor convergence properties, but may be useful for finding initial parameter estimates.
• ネルダー・ミード（単体）探索。この文脈における単体とは、 n次元空間におけるn + 1 個の頂点を持つ多面体、すなわち平面上の三角形、3 次元空間における正四面体などを指します。各頂点は、特定のパラメータセットに対する目的関数の値に対応します。単体の形状とサイズは、最も高い頂点における目的関数の値が常に減少するようにパラメータを変化させることによって調整されます。二乗和は最初は急速に減少する可能性がありますが、MJD Powell の例によれば、準凸問題では非定常点に収束する可能性があります。  

これらの方法やその他の方法に関するより詳細な説明は、『数値計算レシピ』に掲載されており、様々なプログラミング言語によるコンピュータコードも併せて提供されています。