非線形回帰分析（カーブフィッティング）の手法：PythonによるLM法の実装

はじめに

今回は非線形な関数による回帰分析の手法について解説をします。考え方は線形と非線形の場合で違いはなく、最小二乗法の原理で残差平方和 ( 観測値と予測値の差の二乗和 ) を最小にするパラメータを求めます。異なる点は非線形関数の場合には最小にするパラメータが解析的に求められない場合が多いため、数値解析により計算する点です。本記事では数値解析による主な手法である、最急降下法、ニュートン法、ガウス・ニュートン法、レーベンバーグ・マルカート（LM）法について解説します。また、 $${\ 1/(ax+b)+c, \ ae^{bx}+c}$$ の２つの例を用いて実際に LM 法によりパラメータを求めるPythonのコードも掲載しております。

概要

$${n}$$ 個のデータ組 $${(x_{1},y_{1}) ,\ ...\ ,(x_{n},y_{n})}$$ が観測されたとして、$${y}$$ と $${x}$$ について、以下の近似を考えます。

$$
\begin{aligned}
y=f(x;\boldsymbol{a}) \hspace{25pt} (1)
\end{aligned}
$$

$${\boldsymbol{a}}$$ は近似に使用する $${m}$$ 個のパラメータであり、以下のように表します。

$$
\begin{aligned}
\boldsymbol{a} =\begin{pmatrix}
a_{1} \\ \vdots \\ a_{m}
\end{pmatrix} \hspace{25pt} (2)
\end{aligned}
$$

$${f}$$ は $${x}$$ とすべてのパラメータについて２次偏導関数が存在し、連続であるとします（$${C^{2}}$$ 級）。
残差平方和 $${\big(=}$$（実測値 － 予測値 ）の二乗和 $${\big)}$$を $${S(\boldsymbol{a})}$$ とすると、$${S(\boldsymbol{a})}$$ を最小にする $${\boldsymbol{a}(=a_{1},\ ...\ ,a_{m})}$$ の値を求めることが目標です。以下のように、残差平方和 $${S(\boldsymbol{a})}$$を定義します。

$$
\begin{aligned}
S(\boldsymbol{a})
&=\dfrac{1}{2}\sum_{i=1}^{n} \left (y_{i} -f(x_{i};\boldsymbol{a}) \right)^{2}
\hspace{25pt} (3)
\end{aligned}
$$

微分したときに生じる $${2}$$ が邪魔なため、 $${1/2}$$ をかけています。
$${S(\boldsymbol{a})}$$ が最小の場合、以下の式を満たします。

$$
\left\{\begin{aligned}
\dfrac{\partial S(\boldsymbol{a})}{\partial a_{1}} &= S_{a_{1}}(\boldsymbol{a}) =0 \\
& \hspace{5pt} \vdots \hspace{105pt} (4)\\
\dfrac{\partial S(\boldsymbol{a})}{\partial a_{m}} &= S_{a_{m}}(\boldsymbol{a}) =0 \\
\end{aligned}\right.
$$

演算子 $${\nabla}$$ (ナブラ)を用いて以下のようにまとめて表現します。

$$
\begin{aligned}
\nabla S(\boldsymbol{a}) &= \boldsymbol{0} &\hspace{40pt} (5)\\
\nabla S(\boldsymbol{a}) = \dfrac{\partial S(\boldsymbol{a})}{\partial \boldsymbol{a}} &=
\begin{pmatrix}
\dfrac{\partial S(\boldsymbol{a})}{\partial a_{1}}\\
\vdots \\
\dfrac{\partial S(\boldsymbol{a})}{\partial a_{m}}
\end{pmatrix}
= \begin{pmatrix}
S_{a_{1}}(\boldsymbol{a}) \\
\vdots \\
S_{a_{m}}(\boldsymbol{a})
\end{pmatrix} &\hspace{40pt} (6)\\
\end{aligned}
$$

$${\nabla}$$ は多変数関数の１次偏導関数を表します。変数が $${m}$$ 個存在するため、ベクトルになります。勾配ベクトルと呼びます。

パラメータ $${\boldsymbol{a}}$$ について線形の最小二乗法では、上式を解くことができますが、非線形の場合、近似する関数 $${f}$$ によってはパラメータを解析的に求められません。そこで、$${\boldsymbol{a}}$$ について漸化式を立て、それを反復計算することにより、機械的に $${\nabla S(\boldsymbol{a}) = \boldsymbol{0}}$$ となる $${\boldsymbol{a}}$$ を求めます。手法については主として、最急降下法、ニュートン法、ガウス・ニュートン法、レーべンバーグ・マルカート法があり、これらについて解説します。

漸化式については上記のすべての方法について、以下の形で表されます。

$$
\begin{aligned}
　\boldsymbol{a^{(k+1)}}
&= \boldsymbol{a^{(k)}} - M^{-1}
\begin{pmatrix}
S_{a_{1}}(\boldsymbol{a^{(k)}}) \\
\vdots \\
S_{a_{m}}(\boldsymbol{a^{(k)}})
\end{pmatrix} \\
&=\boldsymbol{a^{(k)}} - M^{-1} \nabla S(\boldsymbol{a^{(k)}}) 　\hspace{25pt} (7)
\end{aligned}
$$

$${\boldsymbol{a^{(k+1)}}}$$ は更新後のパラメータ、$${\boldsymbol{a^{(k)}}}$$ は更新前のパラメータ、$${M}$$ は $${m\times m}$$ の正則行列です。

1. 最急降下法

最急降下法の漸化式は $${M=\dfrac{1}{\alpha}I}$$ としたもので、以下のようになります。

$$
\begin{aligned}
　\boldsymbol{a^{(k+1)}}
&=\boldsymbol{a^{(k)}} - \alpha\nabla S(\boldsymbol{a^{(k)}})　\hspace{25pt} (1.1)
\end{aligned}
$$

$${\alpha}$$ は自由に設定可能な係数です。極大値の場合は $${\alpha \lt 0}$$ 極小値の場合は $${\alpha \gt 0}$$ を用います。
$${\alpha}$$ は大きすぎると極値を超えて、発散する場合がありますが、小さすぎると、収束が遅くなります。
最急降下法は $${\alpha \gt 0}$$ として $${\alpha}$$ を十分に小さく設定すれば、更新後に $${S}$$ が増加することはありません。
計算量が少ないという利点がありますが、収束が非常に遅いという欠点があります。

2. ニュートン法

最急降下法の収束の遅さを改善したのがニュートン法です。
ニュートン法の漸化式は、$${M=H_{S}}$$ としたもので、以下のようになります。

$$
\begin{aligned}
　\boldsymbol{a^{(k+1)}}= \boldsymbol{a^{(k)}} - H_{S}^{-1}
\nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (2.1)
\end{aligned}
$$

$${H_{S}}$$ はヘッセ行列と呼び、多変数関数の二次偏導関数を表します。２階偏微分は $${m\times m}$$ 通り存在するため、$${m \times m}$$ の行列になります。偏微分の交換法則により、ヘッセ行列は対称行列になります。$${H_{S}}$$ は以下のように表せます。

$$
\begin{aligned}
H_{S}=\begin{pmatrix}
S_{a_{1}a_{1}}(\boldsymbol{a^{(k)}}) & \cdots & S_{a_{m}a_{1}}(\boldsymbol{a^{(k)}}) \\
\vdots & \ddots & \vdots \\
S_{a_{1}a_{m}}(\boldsymbol{a^{(k)}}) & \cdots & S_{a_{m}a_{m}}(\boldsymbol{a^{(k)}}) \\
\end{pmatrix} \hspace{25pt} (2.2)
\end{aligned}
$$

ニュートン法は $${S(\boldsymbol{a})}$$ を二次式で近似したときの極値を求めるという動作を反復します。
関数 $${S(\boldsymbol{a})}$$ を $${\boldsymbol{a}=\boldsymbol{a^{(k)}}}$$ の周りでテイラー展開し、２次式で近似すると、以下のようになります。

$$
\begin{aligned}
S(\boldsymbol{a})
&\approx S(\boldsymbol{a^{(k)}}) + \nabla S(\boldsymbol{a^{(k)}})^{T} (\boldsymbol{a}-\boldsymbol{a^{(k)}})
+\dfrac{1}{2} (\boldsymbol{a}-\boldsymbol{a^{(k)}})^T H_{S} (\boldsymbol{a}-\boldsymbol{a^{(k)}})
\hspace{25pt} (2.3)
\end{aligned}
$$

上式 $${(2.3)}$$ を $${\boldsymbol{a}}$$ で微分すると、

$$
\begin{aligned}
\nabla S(\boldsymbol{a})
&\approx \nabla S(\boldsymbol{a^{(k)}})
+ H_{S} (\boldsymbol{a}-\boldsymbol{a^{(k)}})=\boldsymbol{0} \hspace{25pt} (2.4)
\end{aligned}
$$

上式 $${(2.4)}$$ について $${\boldsymbol{a}=\boldsymbol{a^{(k+1)}}}$$ を代入し、$${\nabla S(\boldsymbol{a^{(k+1)}})=\boldsymbol{0}}$$ とすると、

$$
\begin{aligned}
\nabla S(\boldsymbol{a^{(k+1)}})
&\approx \nabla S(\boldsymbol{a^{(k)}})
+ H_{S} (\boldsymbol{a^{(k+1)}}-\boldsymbol{a^{(k)}})=\boldsymbol{0} \hspace{25pt} (2.5)
\end{aligned}
$$

よって、以下の式が導かれます。

$$
\begin{aligned}
H_{S} (\boldsymbol{a^{(k+1)}}-\boldsymbol{a^{(k)}}) =-\nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (2.6)
\end{aligned}
$$

以上より、項冒頭の漸化式 $${(2.1)}$$ が導出されます。
実際の計算では、逆行列を計算するよりも連立方程式を解く方が早いため、
$${(2.6)}$$ の連立方程式を解き、$${\boldsymbol{a^{(k+1)}} - \boldsymbol{a^{(k)}}}$$ を求めてから、$${\boldsymbol{a^{(k+1)}}}$$ を求めます。

また、ヘッセ行列 $${H_{S}}$$ の $${ij}$$ 成分 $${S_{a_{j}a_{i}}(\boldsymbol{a})}$$ については以下のように表せます。

$$
\begin{aligned}
S_{a_{j}a_{i}}(\boldsymbol{a})
&= \dfrac{\partial^{2} }{\partial a_{i} \partial a_{j}}S(\boldsymbol{a}) \\
&=\dfrac{\partial^{2} }{\partial a_{i} \partial a_{j}} \sum_{l=1}^{n} \dfrac{1}{2}
\left (y_{l} -f(x_{l};\boldsymbol{a}) \right)^{2} \\
&= \dfrac{\partial}{\partial a_{i}} \sum_{l=1}^{n}
\left (f(x_{l};\boldsymbol{a}) - y_{l} \right) f_{a_{j}}(x_{l};\boldsymbol{a}) \\
&= \sum_{l=1}^{n}
\left \{ f_{a_{i}}(x_{l};\boldsymbol{a})f_{a_{j}}(x_{l};\boldsymbol{a})
+ ( f(x_{l};\boldsymbol{a})
- y_{l})f_{a_{i}a_{j}}(x_{l};\boldsymbol{a}) \right\} \hspace{25pt} (2.7)\\
\end{aligned}
$$

ニュートン法は最急降下法と比較して、収束が早いのが利点ですが、欠点として $${\nabla S(\boldsymbol{a}) = \boldsymbol{0}}$$ となる点に向かっていくだけで、それが極大値、極小値、鞍点かどうかはわからず、$${S}$$ が増加する場合も減少する場合もあります。$${S}$$ が増加するか減少するかはヘッセ行列に依存します。
$${H_{S}}$$ が正定値行列の場合は、$${S(\boldsymbol{a^{(k+1)}})}$$ は減少方向(微小の変化なら減少)、
$${H_{S}}$$ が負定値行列の場合は、$${S(\boldsymbol{a^{(k+1)}})}$$ は増加方向(微小の変化なら増加)
$${H_{S}}$$ が上記以外の場合には、$${S(\boldsymbol{a^{(k+1)}})}$$ の増減方向は場合によって異なります。

以上より、ニュートン法により極小値を計算するには、ヘッセ行列が正定値行列である必要があります。さらに、２階偏微分の計算が必要で、計算量が多いという問題点もあります。

補足1：極大値、極小値とヘッセ行列

行列 $${P}$$ について 任意の実ベクトル $${\boldsymbol{x}\neq \boldsymbol{0}}$$ を用いて、正(負)定値行列の条件は以下のようになります。

1. $${P}$$ が正定値行列　 $${\leftrightarrow \boldsymbol{x}^{T}P\boldsymbol{x} \gt 0 \leftrightarrow}$$ 固有値がすべて正

2. $${P}$$ が半正定値行列 $${\leftrightarrow \boldsymbol{x}^{T}P\boldsymbol{x} \geq 0 \leftrightarrow}$$ 固有値がすべて $${0}$$ 以上

3. $${P}$$ が負定値行列　 $${\leftrightarrow \boldsymbol{x}^{T}P\boldsymbol{x} \lt 0 \leftrightarrow}$$ 固有値がすべて負

4. $${P}$$ が半負定値行列 $${\leftrightarrow \boldsymbol{x}^{T}P\boldsymbol{x} \leq 0 \leftrightarrow}$$ 固有値がすべて $${0}$$ 以下

多変数関数の停留点（勾配ベクトルが零ベクトルとなる点）について

１．極小値

• "$${\nabla S(\boldsymbol{a})=\boldsymbol{0}}$$ かつ $${H_{S}}$$ が半正定値行列" が必要条件

• 正定値行列であれば、必ず極小値

• すべてのパラメータ $${a_{i} (i=1,\ ...\ ,m)}$$ について極小値となる

２．極大値

• "$${\nabla S(\boldsymbol{a})=\boldsymbol{0}}$$ かつ $${H_{S}}$$ が半負定値行列" が必要条件

• 負定値行列であれば、必ず極大値

• すべてのパラメータ $${a_{i} (i=1,\ ...\ ,m)}$$ について極大値となる

３．鞍点

• $${\nabla S(\boldsymbol{a})=\boldsymbol{0}}$$ かつ $${H_{S}}$$ が半正定値行列でも半負定値行列でもなければ、必ず鞍点

• パラメータ $${a_{i} (i=1,\ ...\ ,m)}$$ について極大値となるパラメータと極小値となるパラメータが混在 (固有値がすべてゼロでない場合)

$${(2.3)}$$ の $${S(\boldsymbol{a})}$$ の２次近似式から、$${\boldsymbol{a^{(k)}}}$$ が極値であると仮定し、$${\nabla S(\boldsymbol{a^{(k)}})=\boldsymbol{0}}$$ を代入すると、以下のようになります。

$$
\begin{aligned}
S(\boldsymbol{a})
&\approx S(\boldsymbol{a^{(k)}}) +\dfrac{1}{2} (\boldsymbol{a}-\boldsymbol{a^{(k)}})^T H_{S} (\boldsymbol{a}-\boldsymbol{a^{(k)}}) \hspace{25pt} (2.8)
\end{aligned}
$$

式$${(2.8)}$$ から$${H_{S}}$$ が
正定値行列の場合 $${(\boldsymbol{a}-\boldsymbol{a^{(k)}})^T H_{S} (\boldsymbol{a}-\boldsymbol{a^{(k)}}) \gt 0}$$ となることから、極値 $${\boldsymbol{a^{(k)}}}$$ 付近の $${S(\boldsymbol{a})}$$ は $${S(\boldsymbol{a^{(k)}})}$$ より大きく、
負定値行列の場合 $${(\boldsymbol{a}-\boldsymbol{a^{(k)}})^T H_{S} (\boldsymbol{a}-\boldsymbol{a^{(k)}}) \lt 0}$$ となることから、極値 $${\boldsymbol{a^{(k)}}}$$ 付近の $${S(\boldsymbol{a})}$$ は $${S(\boldsymbol{a^{(k)}})}$$ より小さくなることが分かります。

補足2：正定値（負定値）行列とパラメータ更新後の S の増減

次に、$${M}$$ が正定値（負定値）行列の場合、更新後のパラメータ $${\boldsymbol{a^{(k+1)}}}$$ における残差平方和 $${S(\boldsymbol{a^{(k+1)}})}$$ が減少（増加）方向となることについて考えます。

まず、$${S(\boldsymbol{a})}$$ について、$${\boldsymbol{a^{(k)}}}$$ 周りで一次式で近似すると、

$$
\begin{aligned}
S(\boldsymbol{a})
&\approx S(\boldsymbol{a^{(k)}}) + \nabla S(\boldsymbol{a^{(k)}})^{T} (\boldsymbol{a}-\boldsymbol{a^{(k)}})
\hspace{25pt} (2.9)
\end{aligned}
$$

となります。

上式 $${(2.9)}$$ に $${\boldsymbol{a}= \boldsymbol{a^{(k+1)}}=\boldsymbol{a^{(k)}} - M^{-1} \nabla S(\boldsymbol{a^{(k)}})}$$ を代入すると、

$$
\begin{aligned}
S(\boldsymbol{a^{(k+1)}})
&=S(\boldsymbol{a^{(k)}}-M^{-1}\nabla S(\boldsymbol{a^{(k)}})) \\
&\approx S(\boldsymbol{a^{(k)}}) - \nabla S(\boldsymbol{a^{(k)}})^{T}M^{-1}\nabla S(\boldsymbol{a^{(k)}})
\hspace{25pt} (2.10) \\
\end{aligned}
$$

となります。$${M}$$ が正定値（負定値）の場合、$${M^{-1}}$$ も正定値（負定値）となることから、

$${M}$$ が正定値行列の場合は、$${\nabla S(\boldsymbol{a^{(k)}})^{T}M^{-1}\nabla S(\boldsymbol{a^{(k)}}) \gt 0}$$ となるため、$${S(\boldsymbol{a^{(k+1)}}) \lt S(\boldsymbol{a^{(k)}})}$$ を満たし、
$${M}$$ が負定値行列の場合は、$${\nabla S(\boldsymbol{a^{(k)}})^{T}M^{-1}\nabla S(\boldsymbol{a^{(k)}}) \lt 0}$$ となるため、$${S(\boldsymbol{a^{(k+1)}}) \gt S(\boldsymbol{a^{(k)}})}$$ を満たします。
以上より、$${M}$$ が正定値行列の場合は、減少方向、$${M}$$ が負定値行列の場合は、増加方向となります。

3. ガウス・ニュートン法

ニュートン法の問題点はヘッセ行列によって極値に収束しない場合があり、2階偏微分が必要で計算量が多いことでした。その２つの欠点を改善したのが、ガウス・ニュートン法です。
ガウス・ニュートン法の漸化式は、$${f}$$ のヤコビ行列 $${J_{f}}$$ を用いて $${M=J_{f}^{T}J_{f}}$$ としたもので、以下のようになります。

$$
\begin{aligned}
　\boldsymbol{a^{(k+1)}}= \boldsymbol{a^{(k)}} - (J_{f}^{T}J_{f})^{-1}
\nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (3.1)
\end{aligned}
$$

ここで、$${J_{f}}$$ は $${f}$$ のヤコビ行列であり、以下のように表せます。

$$
\begin{aligned}
J_{f}
&= \begin{pmatrix}
f_{a_{1}}(x_{1};\boldsymbol{a^{(k)}}) & \cdots & f_{a_{m}}(x_{1};\boldsymbol{a^{(k)}}) \\
\vdots & \ddots & \vdots \\
f_{a_{1}}(x_{n};\boldsymbol{a^{(k)}}) & \cdots & f_{a_{m}}(x_{n};\boldsymbol{a^{(k)}})
\end{pmatrix} \hspace{25pt} (3.2)
\end{aligned}
$$

$${J_{f}^{T}J_{f}}$$ は半正定値行列の条件を満たすため、更新後に $${S}$$ が減少方向となります。さらに2階偏微分の計算が不要で計算量も減少します。
ガウス・ニュートン法は、ニュートン法で計算するヘッセ行列の $${ij}$$ 成分 $${S_{a_{j}a_{i}}(\boldsymbol{a})}$$ について計算簡略化のために 2階偏微分の項を除去したものです。

$$
\begin{aligned}
S_{a_{j}a_{i}}(\boldsymbol{a})
&= \sum_{l=1}^{n}
\left \{ f_{a_{i}}(x_{l};\boldsymbol{a})f_{a_{j}}(x_{l};\boldsymbol{a})
+ f_{a_{i}a_{j}}(x_{l};\boldsymbol{a})( f(x_{l};\boldsymbol{a})
- y_{l}) \right\} \\
&\approx \sum_{l=1}^{n} f_{a_{i}}(x_{l};\boldsymbol{a})f_{a_{j}}(x_{l};\boldsymbol{a})
\hspace{25pt} (3.3a) \\
&= \begin{pmatrix}
f_{a_{i}}(x_{1};\boldsymbol{a}) & \cdots & f_{a_{i}}(x_{n};\boldsymbol{a})
\end{pmatrix} \begin{pmatrix}
f_{a_{j}}(x_{1};\boldsymbol{a}) \\ \vdots \\ f_{a_{j}}(x_{n};\boldsymbol{a})
\end{pmatrix} \hspace{25pt} (3.3b)
\end{aligned}
$$

上式からガウス・ニュートン法では、前述のヤコビ行列 $${J_{f}}$$ を用いて、ヘッセ行列 $${H_{S}}$$ を以下のように近似できます。

$$
\begin{aligned}
H_{S}
&\approx J_{f}^{T}J_{f} \hspace{25pt} (3.4)
\end{aligned}
$$

ガウス・ニュートン法においては、$${J_{f}^{T}J_{f}}$$ が半正定値行列の条件を満たすため、値を更新すると減少方向となります。

以下、半正定値行列である証明です。
任意の実ベクトル $${\boldsymbol{x}(\neq \boldsymbol{0})}$$ について、

$$
\begin{aligned}
\boldsymbol{x}^{T} J_{f}^{T}J_{f} \boldsymbol{x}
&= (J_{f}\boldsymbol{x})^{T}J_{f}\boldsymbol{x} \\
&= ||J_{f}\boldsymbol{x}||^{{2}} \geq 0 \hspace{25pt} (3.5)
\end{aligned}
$$

問題点として、逆行列が存在しない場合があります。その場合は漸化式により値を更新できなくなります。

また、 $${\nabla S(\boldsymbol{a^{(k)}})}$$ についてはヤコビ行列と残差のベクトル $${\boldsymbol{r}}$$ を用いて以下のように表す場合もあります。

$$
\begin{aligned}
\nabla S(\boldsymbol{a^{(k)}})
&= \begin{pmatrix}
S_{a_{1}}(\boldsymbol{a^{(k)}}) \\
\vdots \\
S_{a_{m}}(\boldsymbol{a^{(k)}})
\end{pmatrix} \\
&=\begin{pmatrix}
\sum_{l=1}^{n} f_{a_{1}}(x_{l};\boldsymbol{a^{(k)}})
\left (f(x_{l};\boldsymbol{a^{(k)}}) - y_{l} \right) \\
\vdots \\
\sum_{l=1}^{n} f_{a_{m}}(x_{l};\boldsymbol{a^{(k)}})
\left (f(x_{l};\boldsymbol{a^{(k)}}) - y_{l} \right) \\
\end{pmatrix} \\
&= \begin{pmatrix}
f_{a_{1}}(x_{1};\boldsymbol{a^{(k)}}) & \cdots & f_{a_{1}}(x_{n};\boldsymbol{a^{(k)}}) \\
\vdots & \ddots & \vdots \\
f_{a_{m}}(x_{1};\boldsymbol{a^{(k)}}) & \cdots & f_{a_{m}}(x_{n};\boldsymbol{a^{(k)}})
\end{pmatrix}\begin{pmatrix}
f(x_{1};\boldsymbol{a})-y_{1} \\
\vdots \\
f(x_{n};\boldsymbol{a})-y_{n}
\end{pmatrix} \\
&=J_{f}^{T} \boldsymbol{r} \hspace{25pt} (3.6)
\end{aligned}
$$

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

ガウス・ニュートン法の問題点である "逆行列が存在しない場合がある" 点を改善したのがレーベンバーグ・マルカート法です。
レーベンバーグ・マルカート法は $${\lambda \gt 0}$$ として、$${M=J_{f}^{T}J_{f}+\lambda I}$$ としたもので、以下のようになります。

$$
\begin{aligned}
　\boldsymbol{a^{(k+1)}}
= \boldsymbol{a^{(k)}} - &(J_{f}^{T}J_{f}+\lambda I)^{-1}
\nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (4.1)
\end{aligned}
$$

$${J_{f}^{T}J_{f}+\lambda I}$$ については、正定値行列の条件を満たすため（正定値行列は固有値に $${0}$$ がないため正則行列）、途中で値を更新できなくなることがありません。ただし、$${\lambda}$$ も反復計算の度に更新するため、 $${\lambda}$$ が非常に小さい値となると、逆行列が計算できない場合があります。その場合は、$${\lambda}$$ の下限値を設定した方がよいです。

以下、$${J_{f}^{T}J_{f}+\lambda I}$$ が正定値行列である証明です。
任意の実ベクトル $${\boldsymbol{x}(\neq \boldsymbol{0})}$$ について、$${\lambda \gt 0}$$ であることから、

$$
\begin{aligned}
\boldsymbol{x}^{T} (J_{f}^{T}J_{f}+\lambda I) \boldsymbol{x}
&= \boldsymbol{x}^{T} J_{f}^{T}J_{f}\boldsymbol{x} +\lambda \boldsymbol{x}^{T}\boldsymbol{x} \\
&= ||J_{f}\boldsymbol{x}||^{{2}} +\lambda ||\boldsymbol{x}||^{2}\gt 0 \hspace{25pt} (4.2)
\end{aligned}
$$

また、以下の式から、レーベンバーグ・マルカート法はガウス・ニュートン法と最急降下法を混合したものといえます。

$$
\begin{aligned}
& \lambda \approx 0 \\
&\hspace{15pt}\boldsymbol{a^{(k+1)}}
\approx \boldsymbol{a^{(k)}} - (J_{f}^{T}J_{f})^{-1}
\nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (4.3) \\
\\
& \lambda \gg 1 \\
&\hspace{15pt}\boldsymbol{a^{(k+1)}}
\approx \boldsymbol{a^{(k)}} - \dfrac{1}{\lambda}
\nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (4.4) \\
\end{aligned}
$$

以下、$${\lambda}$$ の更新手順の一例です。 $${\alpha \gt 1}$$ とします。
$${\lambda = \lambda_{k}}$$ として計算したパラメータ $${\boldsymbol{a^{(k)}}}$$ を $${\boldsymbol{a^{(k)}} | \lambda_{k}}$$ 、その残差平方和を $${S(\boldsymbol{a^{(k)}}|\lambda_{k})}$$ とします。

前提
$${k}$$ 回目の反復の結果、$${\lambda}$$ については $${\lambda_{k}}$$ 、パラメータについては $${\boldsymbol{a^{(k)}}|\lambda_{k}}$$ が得られたとします。

1. $${\lambda}$$ を $${\lambda_{k}/\alpha}$$ として漸化式から $${\boldsymbol{a^{(k+1)}}}$$ を計算

2. 更新前の残差平方和 $${S(\boldsymbol{a^{(k)}}|\lambda_{k})}$$ と更新後の残差平方和 $${S(\boldsymbol{a^{(k+1)}}|\lambda_{k+1}/\alpha)}$$ を比較

3. 更新後の方が小さい $${\left( S(\boldsymbol{a^{(k+1)}}|\lambda_{k}/\alpha) \boldsymbol{\lt} S(\boldsymbol{a^{(k)}}|\lambda_{k}) \right)}$$ とき、
   $${\to \lambda_{k+1}= \lambda_{k}/\alpha}$$ 更新後のパラメータは $${\boldsymbol{a^{(k+1)}}|\lambda_{k}/\alpha}$$

4. 更新後の方が大きい $${\left( S(\boldsymbol{a^{(k+1)}}|\lambda_{k}/\alpha) \boldsymbol{\gt} S(\boldsymbol{a^{(k)}}|\lambda_{k})\right)}$$ とき、
   $${\to}$$ 更新後の残差平方和の方が更新前よりも小さくなるまで以下の操作を繰り返す。
   　　　　1.　$${\lambda}$$ に $${\alpha}$$ をかける
   　　　　2.　$${\boldsymbol{a^{(k+1)}}}$$ を再計算
   　　　　3.　 更新前後の残差平方和を比較
   　　$${p}$$ を0以上の整数として $${\lambda_{k} \alpha^{p}}$$ の場合に初めて更新後の残差平方和の方が小さくなったとして
   　　$${\lambda_{k+1} = \lambda_{k} \alpha^{p}}$$ 、更新後のパラメータは $${\boldsymbol{a^{(k+1)}}|\lambda_{k}\alpha^{p}}$$

簡潔に言えば、$${p}$$ を $${-1}$$ 以上の整数として、$${S(\boldsymbol{a^{(k+1)}}|\lambda_{k}\alpha^{p}) \boldsymbol{\lt} S(\boldsymbol{a^{(k)}}|\lambda_{k})}$$ を満たす最小の $${p}$$ を求め、$${\lambda_{k+1} = \lambda_{k} \alpha^{p}}$$ 、更新後のパラメータを $${\boldsymbol{a^{(k+1)}}|\lambda_{k}\alpha^{p}}$$ とします。

まとめると、以下のようになります。

1. $${\lambda_{k}/\alpha}$$ として計算したパラメータの残差平方和が、更新前よりも小さくなる場合は $${\lambda}$$ を小さくする
   ⇒　ガウス・ニュートン法の比重を大きくする。

2. 1.を満たさず $${\lambda_{k}}$$ として計算したパラメータの残差平方和が更新前よりも小さくなる場合は $${\lambda}$$ はそのまま
   ⇒　比重はそのまま

3. $${\lambda_{k}}$$ として計算したパラメータの残差平方和が更新前よりも大きい場合は、$${S(\boldsymbol{a^{(k+1)}}) \lt S(\boldsymbol{a^{(k)}})}$$を満たすまで、$${\lambda_{k}}$$ に $${\alpha}$$ をかけ続け再計算 $${=\lambda}$$ を大きくする
   ⇒　最急降下法の比重を大きくする

例1. y = 1/(ax+b)+c の非線形回帰分析

例１として $${y=\dfrac{1}{ax+b}+c}$$ の場合について解説します。

$${\dfrac{a_{1}x+a_{2}}{a_{3}x+a_{4}}+a_{5}}$$ という形も $${\dfrac{1}{ax+b}+c}$$ の形に集約されます。

以下の手順で計算を行います。$${\boldsymbol{a}=(a,b,c)}$$ とします。

1. パラメータ $${a,b,c}$$ の初期値の計算

2. $${\nabla S(\boldsymbol{a^{(k)}})}$$ の計算

3. $${\lambda}$$ を $${\lambda_{k}/\alpha}$$ とする

4. 漸化式 $${(7)}$$ の行列 $${M}$$ の計算

5. $${\boldsymbol{a^{(k+1)}}}$$ を計算

6. $${\boldsymbol{a^{(k+1)}}}$$ の残差平方和 $${S(\boldsymbol{a^{(k+1)}})}$$ を計算

7. $${S(\boldsymbol{a^{(k+1)}}) \lt S(\boldsymbol{a^{(k)}})}$$ となるまで $${\lambda}$$ に $${\alpha}$$ をかけて 4-6 を繰り返す
   $${\to p}$$ を $${-1}$$ 以上の整数として、$${S(\boldsymbol{a^{(k+1)}}|\lambda_{k}\alpha^{p}) \boldsymbol{\lt} S(\boldsymbol{a^{(k)}})}$$ を満たす最小の $${p}$$ を求める

8. $${\lambda_{k+1} = \lambda_{k} \alpha^{p}}$$、更新後のパラメータを $${\boldsymbol{a^{(k+1)}}|\lambda_{k}\alpha^{p}}$$ とする

9. $${\nabla S(\boldsymbol{a^{(k+1)}})}$$ を計算し、許容値以下なら終了、許容値より大きい場合は $${k+1}$$ を $${k}$$ として、3 へ戻る

$${\boldsymbol{a^{(k+1)}} | \lambda_{k} \alpha^{p}}$$：$${\lambda = \lambda_{k}\alpha^{p}}$$ として計算したパラメータ $${\boldsymbol{a^{(k+1)}}}$$

1. パラメータ a,b,c の初期値の計算の補足

良好な計算結果を得るには、初期値を適切に設定する必要があります。初期値をそれぞれ $${a_{0},b_{0},c_{0}}$$ とします。
本例における非線形関数では、$${y=\dfrac{1}{ax+b}+c}$$ より $${z=\dfrac{1}{y-c} (=ax+b)}$$ として $${z}$$ と $${x}$$ が直線になることを利用します。詳細は以下の記事をご覧ください

手順は以下の通りです。

1. 分数関数の形状は"上に凸 / 下に凸"、"増加 / 減少"で $${2 \times 2 = 4}$$ 通り存在　
   → 減少、下に凸 （$${y=1/x}$$ の第１象限の形状） となるよう $${x,y}$$ の符号を調整

2. $${x}$$ と $${z=\dfrac{1}{y-c}}$$ の相関係数 $${\rho_{xz}}$$ が最大となる $${c_{0}}$$ をニュートン法により求める

3. $${z=\dfrac{1}{y-c}}$$ と $${x}$$ の直線近似により線形の最小二乗法から $${a_{0},b_{0}}$$ を求める

相関係数を最大(実際には極大)にする $${c_{0}}$$ については、相関係数 $${\rho_{xz}}$$ の微分 $${\dfrac{d\rho_{xz}}{dc}}$$ から $${c}$$ に無関係なものを除外したものを $${g_{2}(c)}$$ とし、 $${g_{2}(c)=0}$$ となる点 $${c_{0}}$$ をニュートン法により計算します。漸化式は以下の通りです。

$$
\begin{aligned}
c^{(k+1)}= c^{(k)} + \dfrac{g_{2}(c)}{|g_2'(c)|} \hspace{25pt} (\mathrm{ex}1.1)
\end{aligned}
$$

$${g_{2}(c)}$$ が許容値以下になるまで上記の漸化式を反復計算します。極大値を求めるので、$${g_{2}(c)}$$ が正（負）のとき、$${c}$$ が増加（減少）するように $${g'_{2}(c)}$$ に絶対値をとります。$${g_2(c),g'_{2}(c)}$$ は以下の通りです。

$$
\begin{aligned}
g_{2}(c)
&=\left\{ \sum_{i=1}^{n} z_{i}^{2} (x_{i} -\bar{x}) \right\}
\sum_{i=1}^{n}( z_{i} -\bar{z} )^2
-\sum_{i=1}^{n} (x_{i} -\bar{x})z_{i}
\left\{ \sum_{i=1}^{n} z_{i}^{2} ( z_{i} -\bar{z} ) \right\} \hspace{25pt} (\mathrm{ex}1.2) \\
g_{2}'(c)
&=\left\{ 2\sum_{i=1}^{n} z_{i}^3 (x_{i} -\bar{x}) \right\}
\left\{\sum_{i=1}^{n}( z_{i} -\bar{z} )^2 \right\}
+\left\{ \sum_{i=1}^{n} z_{i}^{2} (x_{i} -\bar{x}) \right\}
\left\{ \sum_{i=1}^{n} z_{i}^{2} ( z_{i} -\bar{z} ) \right\} \\
&\hspace{15pt}
-\left\{ \sum_{i=1}^{n} (x_{i} -\bar{x})z_{i} \right\}
\sum_{i=1}^{n} z_{i}^{2} \left\{ 2z_{i}( z_{i} -\bar{z} )
+z_{i}^{2} -\overline{z^{2}} \right\} \hspace{25pt} (\mathrm{ex}1.3)\\
\end{aligned}
$$

$${a,b}$$ の初期値 $${a_{0},b_{0}}$$ については $${z=\dfrac{1}{y-c_{0}},x}$$ の最小二乗法より求めます。

$$
\begin{aligned}
a_{0}
&= \dfrac{\sum_{i=1}^{n} (x_{i}-\bar{x})z_{i}}{\sum_{i=1}^{n} (x_{i}-\bar{x})^{2}} &\hspace{25pt} (\mathrm{ex}1.4)\\
b_{0}
&=\bar{z} -a_{0}\bar{x} &\hspace{25pt} (\mathrm{ex}1.5)\\
\end{aligned}
$$

2. ∇S(a^(k)) の計算の補足

残差平方和 $${S(a,b,c)}$$ および、$${S(a,b,c)}$$ の偏微分はそれぞれ、以下のようになります。

$$
\left\{\begin{aligned}
S(a,b,c) &=\dfrac{1}{2}\sum_{i=1}^{n} \left(y_{i}-\dfrac{1}{a x_{i} + b} -c \right)^{2}
& \hspace{25pt} (\mathrm{ex}1.6) \\
S_{a}(a,b,c) &=\sum_{i=1}^{n} \dfrac{x_{i}}{(a x_{i} + b)^{2}}
\left(y_{i}-\dfrac{1}{a x_{i} + b} -c \right) &\hspace{25pt} (\mathrm{ex}1.6a)\\
S_{b}(a,b,c) &=\sum_{i=1}^{n} \dfrac{1}{(a x_{i} + b)^{2}}
\left(y_{i}-\dfrac{1}{a x_{i} + b} -c \right) &\hspace{25pt} (\mathrm{ex}1.6b) \\
S_{c}(a,b,c) &=\sum_{i=1}^{n} \left(\dfrac{1}{a x_{i} + b} +c -y_{i} \right) &\hspace{25pt} (\mathrm{ex}1.6c)
\end{aligned}\right.
$$

3. 漸化式の行列 M の計算の補足

$${f(x;a,b,c)}$$ および、$${f(x;a,b,c)}$$ の偏微分はそれぞれ、以下のようになります。

$$
\left\{\begin{aligned}
f(x;a,b,c) &=\dfrac{1}{a x + b} +c &\hspace{25pt} (\mathrm{ex}1.7) \\
f_{a}(x;a,b,c) &= -\dfrac{x}{(a x + b)^{2}} &\hspace{25pt} (\mathrm{ex}1.7a) \\
f_{b}(x;a,b,c) &= -\dfrac{1}{(a x + b)^{2}} &\hspace{25pt} (\mathrm{ex}1.7b) \\
f_{c}(x;a,b,c) &= 1 &\hspace{25pt} (\mathrm{ex}1.7c)
\end{aligned}\right.
$$

$${J_{f}^{T}J_{f}}$$ については、以下の式より計算します。

$$
\begin{aligned}
(J_{f}^{T}J_{f})_{ij} =
\sum_{l=1}^{n} f_{a_{i}}(x_{l};\boldsymbol{a^{(k)}})f_{a_{j}}(x_{l};\boldsymbol{a^{(k)}})
&\hspace{25pt} (\mathrm{ex}1.8) \\
\end{aligned}
$$

$$
\begin{aligned}
J_{f}^{T}J_{f}
&= \begin{pmatrix}
\sum_{l=1}^{n} f_{a}^{2} & \sum_{l=1}^{n} f_{a}f_{b} & \sum_{l=1}^{n} f_{a}f_{c} \\
\sum_{l=1}^{n} f_{a}f_{b} & \sum_{l=1}^{n} f_{b}^{2} & \sum_{l=1}^{n} f_{b}f_{c} \\
\sum_{l=1}^{n} f_{a}f_{c} & \sum_{l=1}^{n} f_{b}f_{c} & \sum_{l=1}^{n} f_{c}^{2} \\
\end{pmatrix} \\
&=\begin{pmatrix}
\sum_{l=1}^{n} \dfrac{x_{l}^{2}}{(ax_{l}+b)^{4}} & \sum_{l=1}^{n}\dfrac{x_{l}}{(ax_{l}+b)^{4}} & -\sum_{l=1}^{n} \dfrac{x}{(a x + b)^{2}} \\
\sum_{l=1}^{n} \dfrac{x_{l}}{(ax_{l}+b)^{4}} & \sum_{l=1}^{n}\dfrac{1}{(ax_{l}+b)^{4}} & -\sum_{l=1}^{n} \dfrac{1}{(a x + b)^{2}} \\
-\sum_{l=1}^{n} \dfrac{x}{(a x + b)^{2}} & -\sum_{l=1}^{n} \dfrac{1}{(a x + b)^{2}} & n \\
\end{pmatrix} &\hspace{25pt} (\mathrm{ex}1.9)\\
\end{aligned}
$$

以上より、

$$
\begin{aligned}
M
&=J_{f}^{T}J_{f} + \lambda I
\\
&=\begin{pmatrix}
\sum_{l=1}^{n} \dfrac{x_{l}^{2}}{(ax_{l}+b)^{4}}+\lambda & \sum_{l=1}^{n}\dfrac{x_{l}}{(ax_{l}+b)^{4}} & -\sum_{l=1}^{n} \dfrac{x}{(a x + b)^{2}} \\
\sum_{l=1}^{n} \dfrac{x_{l}}{(ax_{l}+b)^{4}} & \sum_{l=1}^{n}\dfrac{1}{(ax_{l}+b)^{4}}+\lambda & -\sum_{l=1}^{n} \dfrac{1}{(a x + b)^{2}} \\
-\sum_{l=1}^{n} \dfrac{x}{(a x + b)^{2}} & -\sum_{l=1}^{n} \dfrac{1}{(a x + b)^{2}} & n+\lambda \\
\end{pmatrix} &\hspace{25pt} (\mathrm{ex}1.10)\\
\end{aligned}
$$

となります。

4. パラメータ更新の際の計算の補足

$${M}$$ の逆行列の計算は計算量が多いので、以下の連立方程式を解いてから、$${\boldsymbol{a^{(k+1)}}}$$ を求める形にします。

$$
\begin{aligned}
M(\boldsymbol{a^{(k+1)}}-\boldsymbol{a^{(k)}})=-\nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (\mathrm{ex}1.11)
\end{aligned}
$$

連立方程式の解の計算には numpy の linalg.slove 関数を使用します。

Python ソースコード

csv(カンマ区切)、txt(タブ区切)のファイルを読み込み、$${y=\dfrac{1}{ax+b}+c}$$ の非線形回帰のパラメータをレーベンバーグ・マルカート法により計算するプログラムです。近似曲線のグラフと計算結果のcsvファイルを出力します。

バージョン情報
python : 3.12.0
numpy : 2.1.0
pandas : 2.2.2
matplotlib: 3.8.4

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import time
import warnings
import os
import glob

#関数定義
# f, fa , fb , fc , S , Sa , Sb , Sc 

def f(a,x):
    y = 1/(a[0]*x + a[1]) + a[2]
    return y

def f_a(a,x):
    y = - x/(a[0]*x + a[1]) **2 
    return y

def f_b(a,x):
    y = -1/(a[0]*x + a[1])**2
    return y    

def f_c(a,x):
    y = np.ones(len(x))
    return y    

def S(a,x,y):
    r = y - 1/(a[0]*x + a[1]) -a[2]
    return 0.5*np.sum(r**2)

def S_a(a,x,y):
    r = x/(a[0]*x + a[1]) **2 * (y - 1/(a[0]*x+a[1]) -a[2])
    return np.sum(r)

def S_b(a,x,y):
    r = 1/(a[0]*x + a[1]) **2 * (y - 1/(a[0]*x+a[1]) -a[2] )
    return np.sum(r)

def S_c(a,x,y):
    r = 1/(a[0]*x+a[1]) + a[2] -y 
    return np.sum(r)


#数字入力用関数　数字以外は除外
def ENTER_NUM(phrase):
    
    while(1): #数字が入力されるまでループ

        #数字を入力
        num = input(phrase)
        
        if num.isdecimal() : #数値判定
            num = int(num)
            break

        print('数値を入力してください')
    return num

#入力された文字を含むファイル名を出力する関数
# 検索結果が複数の場合は数字を入力
def SEARCH(phrase,exts):
    
    while(1): #ファイル選択完了までループ
        
        word = input(phrase) #検索するファイル名の入力
        word = word.strip() #strip()⇒両端の空白文字除去
        if word == '':
            continue

        #glob でファイルの検索結果（ファイル名）を格納するlist       
        filelist = []

        # 入力となる拡張子すべてに対し検索を実施
        for ext in exts:
            filelist.extend(glob.glob(f'*{word}*{ext}'))

        #ファイルが複数存在する場合は、”番号：ファイル名”を出力し、番号を選択する
        if len(filelist) >= 2:
            print('指定のファイルが複数存在します。')
            for i in range(len(filelist)):
                print(str(i+1)+ ':' +filelist[i])
            print('0:再選択')

            while(1):  #適切なファイル番号が入力されるまでループ

                file_no = ENTER_NUM('番号を選択してください：')
                
                if file_no == 0:
                    break #0が選択されたら、検索の文字入力まで戻る
                
                elif 1 <= file_no <= len(filelist):
                    return filelist[file_no - 1]

                else:
                    continue
            
            continue # 0が入力され、文字検索をやり直す用

        #ファイルが見つからない場合は再度文字の入力から
        elif len(filelist) == 0:
            print('ファイルが見つかりません')
            continue
        
        break
    
    return filelist[0]


# メインの関数
def main():

    # ファイル名の入力
    exts = ['csv','txt'] #検索するファイルの拡張子

    #入力された文字を含むファイルを検索
    infilename = SEARCH('点列データのリンクを入力してください(txt(タブ区切) / csv(カンマ区切))',exts)

    # 選択されたファイルの表示
    print(f'ファイル:"{infilename}"を選択しました。')

    # ファイル名から拡張子の抽出
    ext = infilename.split('.')[-1]

    #拡張子によって区切り文字を選択
    if ext == 'txt':
        delim = '\t'
    elif ext == 'csv':
        delim = ','
    else:
        print(f'拡張子:{ext}は使用できません。txt(タブ区切) / csv(カンマ区切)のファイルを選択してください。')
        return
    

    # 行 / 列番号の選択
    col_x = ENTER_NUM('x:説明変数の列番号を入力：')
    col_y = ENTER_NUM('y:目的変数の列番号を入力：')
    start_row = ENTER_NUM('データが始まる行の番号を入力：')


    #ファイルからx,yのデータを読み込み
    #エラー発生時はmain関数終了
    try:
        x = np.loadtxt(infilename,delimiter=delim,skiprows=start_row-1, usecols=col_x-1)
        y = np.loadtxt(infilename,delimiter=delim,skiprows=start_row-1, usecols=col_y-1)

    except Exception as e:
        print(f'ファイルのデータに数字以外が含まれているか、データの数が一致しません。\n'\
              f'ファイルの中身を確認ください。最初に戻ります。\n' \
              f'ERROR:{e}')
        return
    
    # x,y の平均値を定義
    x_bar = np.average(x)
    y_bar = np.average(y)

    # x,yを 二次式で近似し、二次の項の係数で下に凸か下に凸かを判定
    # yi=β2xi^2 + β0xi +β0 をまとめて 行列 y=Xβ として X の作成
    X = []
    for i in range(len(x)):
        x_i = [1,x[i],x[i]**2] 
        X.append(x_i)

    X =np.array(X)

    # 正規方程式より近似係数を計算 二次項の係数 = beta[2]
    beta = np.linalg.solve(X.T @ X, X.T @ y)    

    # 共分散により増加関数か減少関数か判定
    cov_xy = np.sum((x-x_bar)*y)

    # 減少関数、下に凸となるようにx,yを変換
    sgn_x = -np.sign(cov_xy)*np.sign(beta[2])
    sgn_y = np.sign(beta[2])

    x1 = sgn_x * x
    y1 = sgn_y * y

    # c, g2(c) , g2'(c) の初期値の設定
    c = min(y1)-0.001  # c の初期値は最小値より若干小さい値とする
    g2 =1 #g2(c)の初期値 
    g2_prime = -1 #g2'(c) の初期値
    g2_tol = 1e-3 #g2の許容値　許容値以下になるまで反復繰返し

    rep =0 #反復回数
    rep_limit = 200 #反復回数の上限

    # g2'(c)が負, g2(c) の値が許容値以下になるまで、ニュートン法による反復計算を繰り返す。
    while( (abs(g2) >= g2_tol) or (g2_prime>0) ):

        # y1 -> z に変換
        z = 1/(y1-c)

        #平均,二乗平均等を計算
        x1_bar = np.average(x1)
        z_bar = np.average(z)
        z_sqbar = np.average(z**2)

        #g2(c)の計算
        g2 = np.sum(z**2*(x1-x1_bar)) * np.sum((z-z_bar)**2) \
             - np.sum(z*(x1-x1_bar)) * np.sum(z**2*(z-z_bar))

        #g2'(c)の計算 
        g2_prime = 2*np.sum(z**3*(x1-x1_bar)) * np.sum((z-z_bar)**2) \
               +np.sum(z**2*(x1-x1_bar)) * np.sum(z**2*(z-z_bar)) \
             - np.sum(z*(x1-x1_bar)) * np.sum(z**2*(2*z*(z-z_bar)+z**2-z_sqbar))

        #ニュートン法により値を更新
        c = c + g2/abs(g2_prime)

        # c が y1 の最小値以上の場合、無限大か負値となるため、下限を設定
        if c >= min(y1):
            c = min(y1)-0.001

        #反復回数の更新        
        rep +=1

        # 反復回数が上限に達した場合は計算終了
        if rep == rep_limit:
            print('反復回数の上限に達しました。計算を終了します。')
            break

    # cの計算が終了したら、反復回数と、c,g2(c)の値を表示
    #print(f'{rep}回目：g2={g2}\nc={c}') 

    # z=1/(y-c) と x1 を直線近似し、計算結果をa,bの初期値とする
    c0 = c
    z = 1/(y1-c)

    # x1,z の平均値計算
    x1_bar = np.average(x1)
    z_bar = np.average(z)

    # a,bの初期値a0,b0の計算

    # 最小二乗法
    a0 = (np.dot(x1-x1_bar,z) / np.sum((x1-x1_bar) ** 2))
    b0 = z_bar - x1_bar * a0

    # 符号の修正
    a0 = sgn_x*sgn_y * a0
    b0 = sgn_y * b0
    c0 = sgn_y * c0

    # 3係数を配列にまとめる
    a = np.array([a0,b0,c0])
    #a = np.array([-3,1,1]) #初期値を手動設定する場合

    # 残差平方和の計算   
    sq_err = S(a,x,y)

    # Sの勾配 ∇S(a)の計算
    nabla_S = np.array([S_a(a,x,y),S_b(a,x,y),S_c(a,x,y)])

    #初期値計算結果の表示
    print(f'\n--------初期値の計算結果--------\n'+\
          f'[a,b,c]={a}\n'+\
          f'∇S(a,b,c)={nabla_S}\n'+ \
          f'二乗誤差:{sq_err}\n'+\
          f'--------------------------------')


    # レーベンバーグ・マルカート法による係数の計算、∇S(a)=(Sa,Sb,Sc)が許容値以下になるまで反復計算
    rep = 1 #反復回数
    rep_limit = 1000 #反復回数の上限
    
    error_limit = 1e-10 #許容値
    
    lamb = 1 #λの初期値
    alpha = 2 #αの値
    lamb_lowlim = 1e-7 #λの下限値の設定

    # ∇S(a)の全要素が許容値以下になるまで反復計算
    while ( (abs(nabla_S[0]) >= error_limit) or (abs(nabla_S[1]) >= error_limit) or (abs(nabla_S[2]) >= error_limit) ):

        # ヘッセ行列 =Jf^{T}Jf の各要素の計算 
        Haa = np.sum(f_a(a,x) **2)
        Hab = np.sum(f_a(a,x) * f_b(a,x))
        Hac = np.sum(f_a(a,x) * f_c(a,x))
        Hbb = np.sum(f_b(a,x) **2)
        Hbc = np.sum(f_b(a,x) * f_c(a,x))
        Hcc = np.sum(f_c(a,x) **2)

        # 行列化 = 二次元配列化
        Hs = np.array([[Haa,Hab,Hac],[Hab,Hbb,Hbc],[Hac,Hbc,Hcc]])

        # λをλ/αとする
        lamb = lamb / alpha

        #  λI を加算
        M = Hs + lamb * np.eye(3)

        # 漸化式より係数 a の値を更新
        a1 = a + np.linalg.solve(M,-nabla_S)

        # 更新後、残差平方和が小さくなる場合、λと更新したパラメータはそのままλ=λ/α として、次の更新作業に
        if sq_err > S(a1,x,y):
            pass

        # 更新後、残差平方和が大きくなる場合、
        #更新後の残差平方和が更新前よりも小さくなるまでλにαをかけ続ける
        
        else:
            while(sq_err < S(a1,x,y)): #更新後の残差平方和が更新前よりも小さくなるまでλにαをかけ続ける
                lamb = lamb * alpha
                M = Hs + lamb * np.eye(3)
                a1 = a + np.linalg.solve(M,-nabla_S)      

        a = a1

        #Sの勾配ベクトル ∇S(a,b,c) の値の計算
        nabla_S = np.array([S_a(a,x,y),S_b(a,x,y),S_c(a,x,y)])

        # 残差平方和 S(a,b,c) の計算
        sq_err = S(a,x,y)

        # λが下限値以下になったら、下限値に再設定
        if lamb <= lamb_lowlim:
            lamb = lamb_lowlim
        
        # 反復回数の更新
        rep +=1
        
        # 反復回数が上限に達した場合は計算終了        
        if rep >= rep_limit:
            print('反復回数の上限に達しました。計算を終了します。')
            break

    
    #出力結果格納用配列
    dataset = [a[0],a[1],a[2],nabla_S,sq_err]
    print(f'\n------------計算結果------------\n'+\
          f'反復回数：{rep}\n'+\
          f'[a,b,c]={a}\n'+\
          f'∇S(a,b,c)={nabla_S}\n'+ \
          f'二乗誤差:{sq_err}\n'+\
          f'--------------------------------')


    #'''グラフ出力（黒実線：近似曲線、●：入力データ） グラフ出力が不要の場合は行頭の#を削除
    dx = (max(x)-min(x))/1000
    x1=np.arange(min(x),max(x),dx)
    y_NLLSM = f(a,x1)
    plt.scatter(x,y)
    plt.plot(x1,y_NLLSM,color='black')
    plt.show()
    #'''
                    

    # ファイル出力
    ## 表題作成
    index1 = ['a', 'b','c','∇S','残差平方和']
 
    pd_out = pd.DataFrame(dataset, index=index1)

    # 名前の設定
    outname0 = f'x={int(min(x))}-{int(max(x))}_Non_Linear_LSM'
    outname = outname0 + '.csv'

    #同名ファイルが存在する場合の名前変更
    filenum = 1
    while( os.path.exists(outname) ):
        outname = outname0 + f'_{filenum}.csv'
        filenum += 1

    # csvデータ出力　#shift-jis でないとエクセルで文字化け  
    pd.DataFrame(pd_out).to_csv(outname,header=None, encoding='shift-jis')
    print(f'ファイルを保存しました。ファイル名:{outname}')

       
if __name__ == "__main__":

    while(1):
        main()
        #input('ENTERキーで最初に戻ります。')

例2. y = a exp(bx) + c の回帰分析

例２として、$${y=ae^{bx}+c}$$ の場合について解説します。

$${a_{1}e^{a_{2}(x-a_{3})}+a_{4}}$$ という形も $${ae^{bx}+c}$$ の形に集約されます。

以下の手順で計算を行います。$${\boldsymbol{a}=(a,b,c)}$$ とします。例１とほぼ同じです。初期値の計算方法のみ異なります。

1. パラメータ $${a,b,c}$$ の初期値の計算

2. $${\nabla S(\boldsymbol{a^{(k)}})}$$ の計算

3. $${\lambda}$$ を $${\lambda_{k}/\alpha}$$ とする

4. 漸化式 $${(7)}$$ の行列 $${M}$$ の計算

5. $${\boldsymbol{a^{(k+1)}}}$$ を計算

6. $${\boldsymbol{a^{(k+1)}}}$$ の残差平方和 $${S(\boldsymbol{a^{(k+1)}})}$$ を計算

7. $${S(\boldsymbol{a^{(k+1)}}) \lt S(\boldsymbol{a^{(k)}})}$$ となるまで $${\lambda}$$ に $${\alpha}$$ をかけて 4-6 を繰り返す
   $${\to p}$$ を $${-1}$$ 以上の整数として、$${S(\boldsymbol{a^{(k+1)}}|\lambda_{k}\alpha^{p}) \boldsymbol{\lt} S(\boldsymbol{a^{(k)}})}$$ を満たす最小の $${p}$$ を求める

8. $${\lambda_{k+1} = \lambda_{k} \alpha^{p}}$$、更新後のパラメータを $${\boldsymbol{a^{(k+1)}}|\lambda_{k}\alpha^{p}}$$ とする

9. $${\nabla S(\boldsymbol{a^{(k+1)}})}$$ を計算し、許容値以下なら終了、許容値より大きい場合は $${k+1}$$ を $${k}$$ として、3 へ戻る

$${\boldsymbol{a^{(k+1)}} | \lambda_{k} \alpha^{p}}$$：$${\lambda = \lambda_{k}\alpha^{p}}$$ として計算したパラメータ $${\boldsymbol{a^{(k+1)}}}$$

1. パラメータ a,b,c の初期値の計算の補足

$${a,b,c}$$ の初期値をそれぞれ $${a_{0},b_{0},c_{0}}$$ とします。
例２では、$${y=a e^{bx}+c}$$ より $${z=\ln (y-c) (=bx+\ln a)}$$ として $${z}$$ と $${x}$$ が直線になることを利用します。詳細以下の記事をご覧ください

手順は以下の通りです。

1. 指数関数の形状は "上に凸 / 下に凸"、"増加 / 減少"で $${2 \times 2 = 4}$$ 通り存在　
   → 増加、下に凸 $${(a\gt 0,b \gt 0)}$$ となるよう $${x,y}$$ の符号を調整

2. $${x}$$ と $${z=\ln (y-c)}$$ の相関係数 $${\rho_{xz}}$$ が最大となる $${c_{0}}$$ をニュートン法により求める

3. $${z=\ln (y-c)}$$ と $${x}$$ の直線近似により線形の最小二乗法から $${a_{0},b_{0}}$$ を求める

相関係数を最大(実際には極大)にする $${c_{0}}$$ については、例１と同様に相関係数 $${\rho_{xz}}$$ の微分 $${\dfrac{d\rho_{xz}}{dc}}$$ から $${c}$$ に無関係なものを除外したものを $${g_{2}(c)}$$ として、 $${g_{2}(c)=0}$$ となる点 $${c_{0}}$$ をニュートン法により計算します。漸化式は以下のような式になります。

$$
\begin{aligned}
c^{(k+1)}= c^{(k)} + \dfrac{g_{2}(c)}{|g_2'(c)|} \hspace{25pt} (\mathrm{ex}2.1)
\end{aligned}
$$

$${g_{2}(c)}$$ が許容値以下になるまで上記の漸化式を反復計算します。極大値を求めるので、$${g_{2}(c)}$$ が正（負）のとき、$${c}$$ が増加（減少）するように $${g'_{2}(c)}$$ に絶対値をとります。$${g_2(c),g'_{2}(c)}$$ は以下の通りです。

$$
\begin{aligned}
g_{2}(c)
&=-\left\{ \sum_{i=1}^{n} e^{-z_{i}} (x_{i} -\bar{x}) \right\}
\sum_{i=1}^{n}( z_{i} -\bar{z} )^2
+\sum_{i=1}^{n} (x_{i} -\bar{x})z_{i}
\left\{ \sum_{i=1}^{n} e^{-z_{i}} ( z_{i} -\bar{z} ) \right\}
&\hspace{25pt} (\mathrm{ex}2.2) \\
g_{2}'(c)
&=-\left\{ \sum_{i=1}^{n} e^{-2z_{i}} (x_{i} -\bar{x}) \right\}
\left\{\sum_{i=1}^{n}( z_{i} -\bar{z} )^2 \right\}
+\left\{ \sum_{i=1}^{n} e^{-z_{i}} (x_{i} -\bar{x}) \right\}
\left\{ \sum_{i=1}^{n} e^{-z_{i}} ( z_{i} -\bar{z} ) \right\} \\
&\hspace{15pt}
+\left\{ \sum_{i=1}^{n} (x_{i} -\bar{x})z_{i} \right\}
\sum_{i=1}^{n} e^{-2z_{i}} \left\{ z_{i} -\bar{z}
-1 + e^{z_{i}} \overline{e^{-z}} \right\} \hspace{25pt} &(\mathrm{ex}2.3)\\
\end{aligned}
$$

$${a,b}$$ の初期値 $${a_{0},b_{0}}$$ については $${z=\ln (y-c_{0}),x}$$ の最小二乗法より求めます。

$$
\begin{aligned}
b_{0}
&= \dfrac{\sum_{i=1}^{n} (x_{i}-\bar{x})z_{i}}{\sum_{i=1}^{n} (x_{i}-\bar{x})^{2}}
& \hspace{25pt} (\mathrm{ex}2.4)\\
a_{0}
&=\exp(\bar{z} -b_{0}\bar{x} ) & \hspace{25pt} (\mathrm{ex}2.5)\\
\end{aligned}
$$

2. ∇S(a^(k)) の計算の補足

残差平方和 $${S(a,b,c)}$$ および、$${S(a,b,c)}$$ の偏微分はそれぞれ、以下のようになります。

$$
\left\{\begin{aligned}
S(a,b,c)
&=\sum_{i=1}^{n} \dfrac{1}{2}\left(y_{i}-ae^{bx_{i}} -c \right)^{2} &\hspace{25pt} (\mathrm{ex}2.6) \\
S_{a}(a,b,c)
&=\sum_{i=1}^{n} e^{bx_{i}} ( ae^{bx_{i}} +c -y_{i} ) &\hspace{25pt} (\mathrm{ex}2.6a) \\
S_{b}(a,b,c)
&=\sum_{i=1}^{n} ax_{i}e^{bx_{i}}( ae^{bx_{i}} +c - y_{i} ) &\hspace{25pt} (\mathrm{ex}2.6b) \\
S_{c}(a,b,c)
&=\sum_{i=1}^{n} (ae^{bx_{i}} +c -y_{i} ) &\hspace{25pt} (\mathrm{ex}2.6c)
\end{aligned}\right.
$$

3. 漸化式の行列 M の計算の補足

$${f(x;a,b,c)}$$ および、$${f(x;a,b,c)}$$ の偏微分はそれぞれ、以下のようになります。

$$
\left\{\begin{aligned}
f(x;a,b,c) &=ae^{bx} +c &\hspace{25pt} (\mathrm{ex}2.7) \\
f_{a}(x;a,b,c) &=e^{bx} &\hspace{25pt} (\mathrm{ex}2.7a) \\
f_{b}(x;a,b,c) &= axe^{bx} &\hspace{25pt} (\mathrm{ex}2.7b) \\
f_{c}(x;a,b,c) &= 1 &\hspace{25pt} (\mathrm{ex}2.7c)
\end{aligned}\right.
$$

$${J_{f}^{T}J_{f}}$$ については、$${J_{f}}$$ は作成せず、以下の式より計算します。

$$
\begin{aligned}
J_{f}^{T}J_{f}
&= \begin{pmatrix}
\sum_{l=1}^{n} f_{a}^{2} & \sum_{l=1}^{n} f_{a}f_{b} & \sum_{l=1}^{n} f_{a}f_{c} \\
\sum_{l=1}^{n} f_{a}f_{b} & \sum_{l=1}^{n} f_{b}^{2} & \sum_{l=1}^{n} f_{b}f_{c} \\
\sum_{l=1}^{n} f_{a}f_{c} & \sum_{l=1}^{n} f_{b}f_{c} & \sum_{l=1}^{n} f_{c}^{2} \\
\end{pmatrix} \\
&=\begin{pmatrix}
\sum_{l=1}^{n} e^{2bx} & \sum_{l=1}^{n} axe^{2bx} & \sum_{l=1}^{n} e^{bx} \\
\sum_{l=1}^{n} axe^{2bx} & \sum_{l=1}^{n} a^{2}x^{2}e^{2bx} & \sum_{l=1}^{n} axe^{bx} \\
\sum_{l=1}^{n} e^{bx} & \sum_{l=1}^{n} axe^{bx} & n \\
\end{pmatrix} &\hspace{25pt} (\mathrm{ex}2.8) \\
\end{aligned}
$$

以上より、

$$
\begin{aligned}
M
&=J_{f}^{T}J_{f} + \lambda I \\
&=\begin{pmatrix}
\sum_{l=1}^{n} e^{2bx} +\lambda & \sum_{l=1}^{n} axe^{2bx} & \sum_{l=1}^{n} e^{bx} \\
\sum_{l=1}^{n} axe^{2bx} & \sum_{l=1}^{n} a^{2}x^{2}e^{2bx} +\lambda & \sum_{l=1}^{n} axe^{bx} \\
\sum_{l=1}^{n} e^{bx} & \sum_{l=1}^{n} axe^{bx} & n+\lambda \\
\end{pmatrix} &\hspace{25pt} (\mathrm{ex}2.9)
\end{aligned}
$$

となります。

4. パラメータ更新の際の計算の補足

例１と同様に、$${M}$$ の逆行列の計算は計算量が多いので、以下の連立方程式を解いてから、$${\boldsymbol{a^{(k+1)}}}$$ を求める形にします。

$$
\begin{aligned}
M(\boldsymbol{a^{(k+1)}}-\boldsymbol{a^{(k)}})=-\nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (\mathrm{ex}2.10)
\end{aligned}
$$

連立方程式の解の計算には numpy の linalg.slove 関数を使用します。

Python ソースコード

csv(カンマ区切)、txt(タブ区切)のファイルを読み込み、$${y=ae^{bx}+c}$$ の非線形回帰のパラメータをレーベンバーグ・マルカート法により計算するプログラムです。 近似曲線のグラフと計算結果のcsvファイルを出力します。

バージョン情報
python : 3.12.0
numpy : 2.1.0
pandas : 2.2.2
matplotlib: 3.8.4

import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd
import time
import warnings
import os
import glob

#関数定義
# f, fa , fb , fc , S , Sa , Sb , Sc

def f(a,x):
    y = a[0]*np.exp(a[1]*x) + a[2]
    return y

def f_a(a,x):
    y = np.exp(a[1]*x) 
    return y

def f_b(a,x):
    y = a[0]*x*np.exp(a[1]*x) 
    return y    

def f_c(a,x):
    y = np.ones(len(x))
    return y    

def S(a,x,y):
    r = y - a[0]*np.exp(a[1]*x) - a[2]
    return 0.5*np.sum(r**2)

def S_a(a,x,y):
    r = np.exp(a[1]*x) * (a[0]*np.exp(a[1]*x) +a[2] -y)
    return np.sum(r)

def S_b(a,x,y):
    r = a[0]*x*np.exp(a[1]*x) * (a[0]*np.exp(a[1]*x) +a[2] - y)
    return np.sum(r)

def S_c(a,x,y):
    r = a[0]*np.exp(a[1]*x) +a[2] -y
    return np.sum(r)


#数字入力用関数　数字以外は除外
def ENTER_NUM(phrase):
    
    while(1): #数字が入力されるまでループ

        #数字を入力
        num = input(phrase)
        
        if num.isdecimal() : #数値判定
            num = int(num)
            break

        print('数値を入力してください')
    return num


#入力された文字を含むファイル名を出力する関数
# 検索結果が複数の場合は数字を入力
def SEARCH(phrase,exts):
    
    while(1): #ファイル選択完了までループ
        
        word = input(phrase) #検索するファイル名の入力
        word = word.strip() #strip()⇒両端の空白文字除去
        if word == '':
            continue

        #glob でファイルの検索結果（ファイル名）を格納するlist       
        filelist = []

        # 入力となる拡張子すべてに対し検索を実施
        for ext in exts:
            filelist.extend(glob.glob(f'*{word}*{ext}'))

        #ファイルが複数存在する場合は、”番号：ファイル名”を出力し、番号を選択する
        if len(filelist) >= 2:
            print('指定のファイルが複数存在します。')
            for i in range(len(filelist)):
                print(str(i+1)+ ':' +filelist[i])
            print('0:再選択')

            while(1):  #適切なファイル番号が入力されるまでループ

                file_no = ENTER_NUM('番号を選択してください：')
                
                if file_no == 0:
                    break #0が選択されたら、検索の文字入力まで戻る
                
                elif 1 <= file_no <= len(filelist):
                    return filelist[file_no - 1]

                else:
                    continue
            
            continue # 0が入力され、文字検索をやり直す用

        #ファイルが見つからない場合は再度文字の入力から
        elif len(filelist) == 0:
            print('ファイルが見つかりません')
            continue
        
        break
    
    return filelist[0]

# メインの関数
def main():

    # ファイル名の入力
    exts = ['csv','txt'] #検索するファイルの拡張子

    #入力された文字を含むファイルを検索
    infilename = SEARCH('点列データのリンクを入力してください(txt(タブ区切) / csv(カンマ区切))',exts)

    # 選択されたファイルの表示
    print(f'ファイル:"{infilename}"を選択しました。')

    # ファイル名から拡張子の抽出
    ext = infilename.split('.')[-1]

    #拡張子によって区切り文字を選択
    if ext == 'txt':
        delim = '\t'
    elif ext == 'csv':
        delim = ','
    else:
        print(f'拡張子:{ext}は使用できません。txt(タブ区切) / csv(カンマ区切)のファイルを選択してください。')
        return

    
    # 行 / 列番号の選択
    col_x = ENTER_NUM('x:説明変数の列番号を入力：')
    col_y = ENTER_NUM('y:目的変数の列番号を入力：')
    start_row = ENTER_NUM('データが始まる行の番号を入力：')


    #ファイルからx,yのデータを読み込み
    #エラー発生時はmain関数終了
    try:
        x = np.loadtxt(infilename,delimiter=delim,skiprows=start_row-1, usecols=col_x-1)
        y = np.loadtxt(infilename,delimiter=delim,skiprows=start_row-1, usecols=col_y-1)

    except Exception as e:
        print(f'ファイルのデータに数字以外が含まれているか、データの数が一致しません。\n'\
              f'ファイルの中身を確認ください。最初に戻ります。\n' \
              f'ERROR:{e}')
        return

    # x,y の平均値を定義
    x_bar = np.average(x)
    y_bar = np.average(y)

    # x,yを 二次式で近似し、二次の項の係数で下に凸か下に凸かを判定
    # yi=β2xi^2 + β0xi +β0 をまとめて 行列 y=Xβ として X の作成
    X = []
    for i in range(len(x)):
        x_i = [1,x[i],x[i]**2] 
        X.append(x_i)

    X =np.array(X)

    # 正規方程式より近似係数を計算 二次項の係数 = beta[2]
    beta = np.linalg.solve(X.T @ X, X.T @ y)    

    # 共分散により増加関数か減少関数か判定
    cov_xy = np.sum((x-x_bar)*y)

    # 増加関数、下に凸 (a>0,b>0)となるようにx,yを変換
    sgn_x = np.sign(cov_xy)*np.sign(beta[2])
    sgn_y = np.sign(beta[2])

    x1 = sgn_x * x
    y1 = sgn_y * y


    # ニュートン法により相関係数が最も大きくなるcの計算
    
    # c, g2(c) , g2'(c) の初期値の設定
    c = min(y1)-0.001  # c の初期値は最小値より若干小さい値とする
    g2 =1 #g2(c)の初期値 
    g2_prime = -1 #g2'(c) の初期値
    g2_tol = 1e-3 #g2の許容値　許容値以下になるまで反復繰返

    rep =0 #反復回数
    rep_limit = 200 #反復回数の上限

    # g2'(c)が負, g2(c) の値が許容値以下になるまで、ニュートン法による反復計算を繰り返す。
    while( (abs(g2) >= g2_tol) or (g2_prime>0) ):

        # y1 -> z に変換
        z = np.log((y1-c))

        #平均等を計算        
        x1_bar = np.average(x1)
        z_bar = np.average(z)
        exp_minusz = np.exp(-z)
        exp_minusz_bar = np.average(np.exp(-z))

        #g2(c)の計算
        g2 = -np.sum(np.exp(-z)*(x1-x1_bar)) * np.sum((z-z_bar)**2) \
             + np.sum(z*(x1-x1_bar)) * np.sum(np.exp(-z)*(z-z_bar))

        #g2'(c)の計算 
        g2_prime = -np.sum(np.exp(-2*z)*(x1-x1_bar)) * np.sum((z-z_bar)**2) \
               +np.sum(np.exp(-z)*(x1-x1_bar)) * np.sum(np.exp(-z)*(z-z_bar)) \
             + np.sum(z*(x1-x1_bar)) * np.sum(np.exp(-2*z)*(z-z_bar -1+np.exp(z)*exp_minusz_bar))

        #ニュートン法により値を更新
        c = c + g2/abs(g2_prime)


        # c が y1 の最小値以上の場合、無限大か負値となるため、下限を設定
        if c >= min(y1):
            c = min(y1)-0.001

        #反復回数の更新        
        rep +=1

        # 反復回数が上限に達した場合は計算終了
        if rep == rep_limit:
            print('反復回数の上限に達しました。計算を終了します。')
            break

        
    # cの計算が終了したら、反復回数と、c,g2(c)の値を表示
    #print(f'{rep}回目：g2={g2}\nc={c}') 

    # z=ln(y-c) と x1 を直線近似し、計算結果をa,bの初期値とする
    c0 = c
    z = np.log(y1-c)


    # x1,z の平均値計算
    x1_bar = np.average(x1)
    z_bar = np.average(z)

    # a,bの初期値a0,b0の計算

    # 最小二乗法による係数の計算
    slp = (np.dot(x1-x1_bar,z) / np.sum((x1-x1_bar) ** 2))
    intcep = z_bar - x1_bar * slp

    a0 = np.exp(intcep)
    b0 = slp
    
    # 符号の修正
    a0 = sgn_y * a0
    b0 = sgn_x * b0
    c0 = sgn_y * c0

    # 3係数を配列にまとめる
    a = np.array([a0,b0,c0])
    #a = np.array([-1,-1,-16]) #初期値を手動設定する場合

    # 残差平方和の計算   
    sq_err = S(a,x,y)

    # Sの勾配 ∇S(a)の計算
    nabla_S = np.array([S_a(a,x,y),S_b(a,x,y),S_c(a,x,y)])

    #初期値計算結果の表示
    print(f'\n--------初期値の計算結果--------\n'+\
          f'反復回数：{rep}\n'+\
          f'[a,b,c]={a}\n'+\
          f'∇S(a,b,c)={nabla_S}\n'+ \
          f'二乗誤差:{sq_err}\n'+\
          f'--------------------------------')

    # レーベンバーグ・マルカート法による係数の計算、∇S(a)=(Sa,Sb,Sc)が許容値以下になるまで反復計算
    rep = 1 #反復回数
    rep_limit = 1000 #反復回数の上限
    
    error_limit = 1e-8 #許容値
    
    lamb = 1 #λの初期値
    alpha = 2 #αの値
    lamb_lowlim = 1e-6 #λの下限値の設定

    # ∇S(a)の全要素が許容値以下になるまで反復計算
    while ( (abs(nabla_S[0]) >= error_limit) or (abs(nabla_S[1]) >= error_limit) or (abs(nabla_S[2]) >= error_limit) ):

        # ヘッセ行列 =Jf^{T}Jf の各要素の計算 
        Haa = np.sum(f_a(a,x) **2)
        Hab = np.sum(f_a(a,x) * f_b(a,x))
        Hac = np.sum(f_a(a,x) * f_c(a,x))
        Hbb = np.sum(f_b(a,x) **2)
        Hbc = np.sum(f_b(a,x) * f_c(a,x))
        Hcc = np.sum(f_c(a,x) **2)

        # 行列化 = 二次元配列化
        Hs = np.array([[Haa,Hab,Hac],[Hab,Hbb,Hbc],[Hac,Hbc,Hcc]])

        # λをλ/αとする
        lamb = lamb / alpha

        #  λI を加算
        M = Hs + lamb * np.eye(3)

        # 漸化式より係数 a の値を更新
        a1 = a + np.linalg.solve(M,-nabla_S)

        # 更新後、残差平方和が小さくなる場合、λと更新したパラメータはそのままλ=λ/α として、次の更新作業に
        if sq_err > S(a1,x,y):
            pass

        # 更新後、残差平方和が大きくなる場合、
        #更新後の残差平方和が更新前よりも小さくなるまでλにαをかけ続ける
        
        else:
            while(sq_err < S(a1,x,y)): #更新後の残差平方和が更新前よりも小さくなるまでλにαをかけ続ける
                lamb = lamb * alpha
                M = Hs + lamb * np.eye(3)
                a1 = a + np.linalg.solve(M,-nabla_S)      

        a = a1

        #Sの勾配ベクトル ∇S(a,b,c) の値の計算
        nabla_S = np.array([S_a(a,x,y),S_b(a,x,y),S_c(a,x,y)])

        # 残差平方和 S(a,b,c) の計算
        sq_err = S(a,x,y)

        # λが下限値以下になったら、下限値に再設定
        if lamb <= lamb_lowlim:
            lamb = lamb_lowlim
        
        # 反復回数の更新
        rep +=1
        
        # 反復回数が上限に達した場合は計算終了        
        if rep >= rep_limit:
            print('反復回数の上限に達しました。計算を終了します。')
            break

    
    #出力結果格納用配列
    dataset = [a[0],a[1],a[2],nabla_S,sq_err]
    print(f'\n------------計算結果------------\n'+\
          f'反復回数：{rep}\n'+\
          f'[a,b,c]={a}\n'+\
          f'∇S(a,b,c)={nabla_S}\n'+ \
          f'二乗誤差:{sq_err}\n'+\
          f'--------------------------------')


    #'''グラフ出力（黒実線：近似曲線、●：入力データ） グラフ出力が不要の場合は行頭の#を削除
    dx = (max(x)-min(x))/1000
    x1=np.arange(min(x),max(x),dx)
    y_NLLSM = f(a,x1)
    plt.scatter(x,y)
    plt.plot(x1,y_NLLSM,color='black')
    plt.show()
    #'''
             
    # ファイル出力
    ## 表題作成
    index1 = ['a', 'b','c','∇S','残差平方和']
 
    pd_out = pd.DataFrame(dataset, index=index1)

    # 名前の設定
    outname0 = f'x={int(min(x))}-{int(max(x))}_Non_Linear_LSM'
    outname = outname0 + '.csv'

    # 同名ファイルが存在する場合の名前変更
    filenum = 1
    while( os.path.exists(outname) ):
        outname = outname0 + f'_{filenum}.csv'
        filenum += 1
    
    # csvデータ出力　#shift-jis でないとエクセルで文字化け
    pd.DataFrame(pd_out).to_csv(outname,header=None, encoding='shift-jis')
    print(f'ファイルを保存しました。ファイル名:{outname}')


if __name__ == "__main__":

    while(1):
        main()
        #input('ENTERキーで最初に戻ります。')