【非線形回帰分析】信頼領域法について

はじめに

今回は非線形な関数による回帰分析（モデルとなる関数のパラメータの計算）の手法の中で信頼領域法について解説をします。信頼領域法は ”対象となる関数（以下、目的関数）を二次式で近似し、その近似関数を最小とするパラメータを求める” という動作を反復します。従来の方法と異なる点は更新の際に加算するベクトルの大きさに制限を設ける点です。制限する大きさについては、目的関数とその近似関数のパラメータ更新による減少量の比によって決定します。最急降下法やガウス・ニュートン法など従来の計算方法に関する知識が必要ですので、こちらの記事を予めご一読頂ければ幸いです。

1. 前置き

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

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

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

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

$f$ は $x$ とパラメータ $\boldsymbol{a}$ について二次偏導関数が存在し、連続であるとします（$C^{2}$ 級）。
また、残差平方和 $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} (1.3) \end{aligned}$$

$S(\boldsymbol{a})$ を最小とする $\boldsymbol{a}$ を信頼領域法により求めるのが今回の目標です。$S(\boldsymbol{a})$ が最小の場合、以下の式を満たします。

$$\begin{aligned} \nabla S(\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} = \boldsymbol{0} \hspace{40pt} (1.4) \end{aligned}$$

2. 概要

信頼領域法は従来の方法と同様に漸化式の反復計算により求めます。
漸化式については、$k$ 回目の計算によって得られたパラメータを $\boldsymbol{a^{(k)}}$ 、更新の際に加算するベクトル（以下、探索ベクトル）を $\boldsymbol{d^{(k)}}$ として、以下のように表せます。

$$\begin{aligned} \boldsymbol{a^{(k+1)}} = \boldsymbol{a^{(k)}} + \boldsymbol{d^{(k)}} \hspace{40pt} (2.1) \end{aligned}$$

信頼領域法では、探索ベクトル $\boldsymbol{d^{(k)}}$ の大きさが $r_{k}$ 以下となるように制約を設けた上で、$S(\boldsymbol{a})$ の二次近似式を最小とする $\boldsymbol{d^{(k)}}$ を計算します。$r_{k}$ を信頼半径、探索ベクトルの大きさが信頼半径以下である領域を信頼領域と呼びます。信頼半径は以下に示す評価関数 $\rho_{k}$ の値によって決定されます。 $\rho_{k}$は "目的関数（=残差平方和 $S(\boldsymbol{a})$）の減少量" $÷$ "二次近似関数の減少量" で計算される値です。$S(\boldsymbol{a})$ の二次近似関数を $S_{(2)}(\boldsymbol{a})$ とします。

$$\begin{aligned} \rho_{k} &= \dfrac{S(\boldsymbol{a^{(k)}+d^{(k)}})-S(\boldsymbol{a^{(k)}})} {S_{(2)}(\boldsymbol{a^{(k)}+d^{(k)}})-S_{(2)}(\boldsymbol{a^{(k)}})} \hspace{40pt} (2.2)\\ \end{aligned}$$

信頼領域法のおおまかな手順は、以下のようになります。

1. パラメータの初期値 $\boldsymbol{a^{(0)}}$ 、信頼半径の初期値 $r_{0}$ 、閾値 $\eta$ を設定する
2. $S(\boldsymbol{a})$ を二次式で近似した関数 $S_{(2)}(\boldsymbol{a})$ を最小とする探索ベクトル $\boldsymbol{d^{(k)}}$ を $\boldsymbol{d^{(k)}} \leq r_{k}$ の範囲で計算する
   $\to$ ドッグレッグ法により計算
3. 目的関数とその近似関数のズレを評価する関数 $\rho_{k}$ を計算
4. 信頼半径を更新
   [1] $\rho_{k} \lt 0.25$ の場合、$r_{k+1}=0.25 r_{k}$
   [2] $\rho_{k} \gt 0.75$ かつ $||\boldsymbol{d^{(k)}}|| = r_{k}$ の場合、$r_{k+1}=2r_{k}$
   [3] 1,2を満たさない場合、 $r_{k+1}=r_{k}$
5. パラメータの値を更新
   [1] $\rho_{k} \geq \eta$ の場合、 $\boldsymbol{a^{(k+1)}} = \boldsymbol{a^{(k)}} + \boldsymbol{d^{(k)}}$ よりパラメータの値を更新
   [2] $\rho_{k} \lt \eta$ の場合、 $\boldsymbol{a^{(k+1)}} = \boldsymbol{a^{(k)}}$ ⇒パラメータはそのまま
6. $\nabla S(\boldsymbol{a^{(k+1)}})$ を計算し、許容値以下なら終了、そうでない場合は $k+1 \to k$ として 2 へ戻る

3. ドッグレッグ法

信頼領域法は上記の手順２で示したように、"二次近似関数が最小値となる探索ベクトルを信頼領域の範囲内で求める" といった制限付きの最適化問題が存在します。これを解くにあたってはドッグレッグ法と呼ばれる手法を用います。

まず、$S(\boldsymbol{a})$ を二次式で近似します。 $S$ のヘッセ行列 $H_{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)}}) +\dfrac{1}{2} (\boldsymbol{a}-\boldsymbol{a^{(k)}})^T H_{S} (\boldsymbol{a}-\boldsymbol{a^{(k)}}) \hspace{25pt} (3.1) \end{aligned}$$

ヘッセ行列については二階偏微分の計算の複雑さからガウス・ニュートン法同様に、関数 $f$ に対するヤコビ行列 $J_{f}$ を用いて以下のように近似します。

$$\begin{aligned} H_{S} = J_{f}^{T}J_{f} \hspace{25pt} (3.2) \end{aligned}$$

ヤコビ行列 $J_{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.3) \end{aligned}$$

以上をまとめると、二次近似関数 $S_{(2)}(\boldsymbol{a})$ は以下のようになります。

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

$||\boldsymbol{a}-\boldsymbol{a^{(k)}}||=||\boldsymbol{d^{(k)}}|| \leq r_{k}$ の範囲で $S$ の二次近似式 $S_{(2)}$ を最小とする探索ベクトル $\boldsymbol{d^{(k)}}$ をドッグレッグ法により計算します。

ドッグレッグ法の簡易手順を以下に示します。

1. ガウス・ニュートンステップ $\boldsymbol{d_{gn}}$ を計算
2. 1 で求めたガウス・ニュートンステップ $\boldsymbol{d_{gn}}$ と信頼半径 $r_{k}$ の大小を比較し、
   [1] $||\boldsymbol{d_{gn}}|| \leq r_{k}$ の場合 　⇒　$\boldsymbol{d^{(k)}}=\boldsymbol{d_{gn}}$
   [2] $||\boldsymbol{d_{gn}}|| \gt r_{k}$ の場合 　⇒　3 へ
3. コーシー点 $\boldsymbol{d_{c}}$ を計算
4. 3 で求めたコーシー点 $\boldsymbol{d_{c}}$ の大きさと信頼半径 $r_{k}$ の大小を比較し、
   [1] $||\boldsymbol{d_{c}}|| \lt r_{k}$ の場合　⇒　$\boldsymbol{d^{(k)}}=\boldsymbol{d_{dl}}$ ⇒ドッグレッグステップ
   [2] $||\boldsymbol{d_{c}}|| \geq r_{k}$ の場合　⇒　$\boldsymbol{d^{(k)}}=\dfrac{r_{k}}{||\boldsymbol{d_{c}}||}\boldsymbol{d_{c}}$ (方向はコーシー点と同じで大きさが $r_{k}$)

ガウス・ニュートンステップ $\boldsymbol{d_{gn}}$ 、コーシー点 $\boldsymbol{d_{c}}$ 、ドッグレッグステップ $\boldsymbol{d_{dl}}$ の式は以下の通りです。

$$\left\{\begin{aligned} \boldsymbol{d_{gn}} &= -(J_{f}^{T}J_{f})^{-1}\nabla S(\boldsymbol{a^{(k)}}) &\hspace{2pt} (3.5a) \\ \boldsymbol{d_{c}} &= -\alpha_{0} \nabla S(\boldsymbol{a^{(k)}}) *_{1} &\hspace{2pt} (3.5b) \\ \boldsymbol{d_{dl}} &= \beta \boldsymbol{d_{gn}}+(1-\beta)\boldsymbol{d_{c}} *_{2} &\hspace{2pt} (3.5c) \end{aligned}\right.$$

$*_{1}$ $\alpha_{0} = \dfrac{||\nabla S(\boldsymbol{a^{(k)}})||^{2}}{||J_{f}\nabla S(\boldsymbol{a^{(k)}})||^{2}}$ 　　$\alpha_{0}$ は最急降下法の係数で、$S_{(2)}(\boldsymbol{a} -\alpha \nabla S(\boldsymbol{a^{(k)}}))$ を最小にする $\alpha$
$*_{2}$ $\beta$ は $||\boldsymbol{d_{dl}}||=r_{k}$ を満たすように設定

3.1 ガウス・ニュートンステップ

まず最初に、ガウス・ニュートンステップを計算します。ガウス・ニュートンステップは二次近似関数 $S_{(2)}(\boldsymbol{a})$ の勾配ベクトル $\nabla S_{(2)}(\boldsymbol{a})=\boldsymbol{0}$ となる点のことです。ガウス・ニュートン法と同じものです。

$S_{(2)}(\boldsymbol{a})$ の勾配ベクトルを求め、$\boldsymbol{0}$ とすると、

$$\begin{aligned} \nabla S_{(2)}(\boldsymbol{a}) &= \nabla S(\boldsymbol{a^{(k)}}) + J_{f}^{T}J_{f} (\boldsymbol{a}-\boldsymbol{a^{(k)}}) = \boldsymbol{0} \hspace{25pt} (3.1.1) \end{aligned}$$

よって、勾配ベクトル $=\boldsymbol{0}$ となる $\boldsymbol{a}$ は以下のようになります。ガウス・ニュートン法の漸化式と全く同じです。

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

以上より、ガウス・ニュートンステップ $\boldsymbol{d_{gn}}$ は以下のようになります。

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

3.2 コーシー点

ガウス・ニュートンステップが信頼領域外に存在する場合は、コーシー点を計算します。コーシー点は最急降下法と同様に探索ベクトルを $-\alpha_{0} \nabla S(\boldsymbol{a^{(k)}})$ とします。勾配ベクトルにかけ合わせる係数 $\alpha_{0}$ については、$S_{(2)}(\boldsymbol{a}-\alpha \nabla S(\boldsymbol{a^{(k)}}))$ を $\alpha$ の関数として考えた場合に $S_{(2)}$ を最小にする $\alpha$ です。

$S_{(2)}(\boldsymbol{a}-\alpha \nabla S(\boldsymbol{a^{(k)}}))$ を整理すると以下のようになります。

$$\begin{aligned} S_{(2)}(\boldsymbol{a}-\alpha \nabla S(\boldsymbol{a^{(k)}})) &= S(\boldsymbol{a^{(k)}}) - \nabla S(\boldsymbol{a^{(k)}})^{T} \alpha \nabla S(\boldsymbol{a^{(k)}}) +\dfrac{1}{2} (\alpha \nabla S(\boldsymbol{a^{(k)}}))^T (J_{f}^{T}J_{f}) (\alpha \nabla S(\boldsymbol{a^{(k)}})) \\ &=S(\boldsymbol{a^{(k)}})-\alpha ||\nabla S(\boldsymbol{a^{(k)}})||^{2} +\dfrac{\alpha^{2}}{2}||J_{f}\nabla S(\boldsymbol{a^{(k)}})||^{2} \\ &=S(\boldsymbol{a^{(k)}})+\dfrac{||J_{f}\nabla S(\boldsymbol{a^{(k)}})||^{2}}{2} \left(\alpha - \dfrac{||\nabla S(\boldsymbol{a^{(k)}})||^{2}}{||J_{f}\nabla S(\boldsymbol{a^{(k)}})||^{2}} \right)^{2} +o(\alpha) \hspace{25pt} (3.2.1) \end{aligned}$$

上式を最小とする $\alpha_{0}$ は

$$\begin{aligned} \alpha_{0} = \dfrac{||\nabla S(\boldsymbol{a^{(k)}})||^{2}}{||J_{f}\nabla S(\boldsymbol{a^{(k)}})||^{2}} \hspace{25pt} (3.2.2) \end{aligned}$$

となります。以上より、コーシー点は以下のように表せます。

$$\begin{aligned} \boldsymbol{d_{c}} = - \dfrac{||\nabla S(\boldsymbol{a^{(k)}})||^{2}} {||J_{f}\nabla S(\boldsymbol{a^{(k)}})||^{2}} \nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (3.2.3) \end{aligned}$$

3.3 ドッグレッグステップ

コーシー点が信頼領域内にある場合、 ガウス・ニュートンステップ、コーシー点を結ぶ直線と信頼領域の境界との交点を使用します。この交点をドッグレッグステップと呼びます。数式で表すと以下のようになります。$\beta$ は $||\boldsymbol{d_{dl}}|| = r_{k}$ となるように設定します。

$$\begin{aligned} \boldsymbol{d_{dl}} &= \boldsymbol{d_{c}} + \beta( \boldsymbol{d_{gn}}- \boldsymbol{d_{c}})\\ &= \beta\boldsymbol{d_{gn}} + (1-\beta) \boldsymbol{d_{c}}\\ &= - \left\{ \beta (J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta) I \right\} \nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (3.3.1) \end{aligned}$$

上式 $(3.3.1)$ よりドッグレッグステップはレーベンバーグ・マルカート（LM）法と類似した形となります。

$\beta$ については $||\boldsymbol{d_{dl}}||^{2}$ を二次方程式の形に変形し、解の公式から求めます。($0 \lt \beta \lt 1$)

まず、$||\boldsymbol{d_{dl}}||^{2}$ を変形すると

$$\begin{aligned} ||\boldsymbol{d_{dl}}||^{2} &=\boldsymbol{d_{dl}}^{T}\boldsymbol{d_{dl}}\\ &= \{ \boldsymbol{d_{c}}^{T} + \beta( \boldsymbol{d_{gn}}- \boldsymbol{d_{c}})^{T}\} \{ \boldsymbol{d_{c}} + \beta( \boldsymbol{d_{gn}}- \boldsymbol{d_{c}})\} \\ &=||\boldsymbol{d_{gn}}-\boldsymbol{d_{c}}||^{2}\beta^{2} +2\beta \boldsymbol{d_{c}}^{T}( \boldsymbol{d_{gn}}- \boldsymbol{d_{c}}) +||\boldsymbol{d_{c}}||^{2} \hspace{25pt} (3.3.2) \\ \end{aligned}$$

となります。 $||\boldsymbol{d_{dl}}||^{2}=r_{k}^{2}$ より解の公式から、$\beta$ は以下のようになります。

$$\begin{aligned} \beta = \dfrac{-\boldsymbol{d_{c}}^{T}( \boldsymbol{d_{gn}}- \boldsymbol{d_{c}}) \pm \sqrt{\left(\boldsymbol{{d_{c}}}^{T}( \boldsymbol{d_{gn}}- \boldsymbol{d_{c}}) \right)^{2} +||\boldsymbol{d}_{gn}-\boldsymbol{d}_{c}||^{2} \left(r_{k}^{2}-\boldsymbol{|| \boldsymbol{d}_{c}||}^{2} \right)}} {||\boldsymbol{d}_{gn}-\boldsymbol{d}_{c}||^{2}} \hspace{25pt} (3.3.3) \end{aligned}$$

解は2つ存在しますが、$0\lt \beta \lt 1$ の範囲の解を選びます。

コーシー点ではなくドッグレッグステップを使う理由はドッグレッグステップにおける二次近似関数 $S_{(2)}$ の値が コーシー点における $S_{(2)}$ の値以下となるためです。証明は 6 章の補足をご覧ください。

4. 信頼半径・パラメータの更新

信頼半径の更新には前章にて計算した探索ベクトル $\boldsymbol{d^{(k)}}$ を用いて以下の評価関数 $\rho_{k}$ を使用します。評価関数 $\rho_{k}$ は残差平方和 $S(\boldsymbol{a})$ とその二次近似関数 $S_{(2)}(\boldsymbol{a})$ の減少量の比を表します。

$$\begin{aligned} \rho_{k} &= \dfrac{S(\boldsymbol{a^{(k)}+d^{(k)}})-S(\boldsymbol{a^{(k)}})} {S_{(2)}(\boldsymbol{a^{(k)}+d^{(k)}})-S_{(2)}(\boldsymbol{a^{(k)}})} \\ &= \dfrac{S(\boldsymbol{a^{(k)}+d^{(k)}})-S(\boldsymbol{a^{(k)}})} {S_{(2)}(\boldsymbol{a^{(k)}+d^{(k)}})-S(\boldsymbol{a^{(k)}})} \\ &= \dfrac{S(\boldsymbol{a^{(k)}+d^{(k)}})-S(\boldsymbol{a^{(k)}})} {\nabla S(\boldsymbol{a^{(k)}})^{T} \boldsymbol{d^{(k)}} +\dfrac{1}{2} (\boldsymbol{d^{(k)}})^T J_{f}^{T}J_f \boldsymbol{d^{(k)}} } \\ &= \dfrac{S(\boldsymbol{a^{(k)}+d^{(k)}})-S(\boldsymbol{a^{(k)}})} {\nabla S(\boldsymbol{a^{(k)}})^{T} \boldsymbol{d^{(k)}} +\dfrac{1}{2} ||J_f \boldsymbol{d^{(k)}} ||^{2} } \hspace{25pt} (4.1)\\ \end{aligned}$$

上記 $\rho_{k}$ の計算結果によって、以下のように信頼半径およびパラメータの更新をします。

[1] 信頼半径の更新

1. $\rho_{k} \lt 0.25$ のとき
   $r_{k+1} = 0.25 r_{k}$

2. $\rho_{k} \gt 0.75$ かつ $||\boldsymbol{d^{(k)}}||=r_{k}$
   $r_{k+1} = 2r_{k}$
   ※信頼半径の上限値 $r_{max}$ 超える場合は上限値を使用

3. 1 および 2 の条件を満たさないとき
   $r_{k+1} = r_{k}$

目的関数とその近似関数のズレが大きい場合には信頼半径を小さく、ズレが小さくガウス・ニュートンステップが信頼領域外にある場合に信頼半径を大きくします。

[2] パラメータの更新
$0 \sim 0.25$ の判定で設定した閾値 $\eta$ を用います。

1. $\rho_{k} \geq \eta$ のとき
   $\boldsymbol{a^{(k+1)}} = \boldsymbol{a^{(k)}} +\boldsymbol{d^{(k)}}$

2. $\rho_{k} \lt \eta$ のとき
   $\boldsymbol{a^{(k+1)}} = \boldsymbol{a^{(k)}}$

目的関数とその近似関数のズレが大きい場合には、パラメータ $\boldsymbol{a}$ は更新せず、信頼半径を小さくして再計算します。

5. 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-6 #許容値
    
    
    r_tr = 10 #信頼半径の初期値
    eta = 0.2 #閾値 η
    r_max = 100 #信頼半径の上限値の設定
    lamb_min = 1e-10 #ガウス・ニュートン点が求められない（ヘッセ行列の近似列が正則でない）場合に加える単位行列の係数

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

        # ドッグレッグ法による計算

        #ヤコビ行列の計算
        Jf = np.array([f_a(a,x),f_b(a,x),f_c(a,x)]).T
        #print(f'ヤコビ行列Jf\n{Jf}')

        # ヘッセ行列 =Jf^{T}Jf の計算 
        Hs = Jf.T @ Jf
        #print(f'ヘッセ近似行列Hs\n{Hs}')

        lamb = lamb_min
        M = Hs

        # ガウス・ニュートン点の計算
        while(1):
            try:
                d_gn = np.linalg.solve(M,-nabla_S)
                break
                
            except Exception as e:
                #ヘッセ近似行列が正則でない場合、単位行列を足して正則化
                M = Hs + lamb * np.eye(3)
                lamb = lamb * 2
                print(f'ヘッセ近似行列が正則でない可能性があります。正則化して再計算します。\n' +\
                      f'ERROR:{e}')
            

        # ガウス・ニュートンステップの大きさと信頼半径を比較
        dgn_norm = np.linalg.norm(d_gn)

        if dgn_norm <= r_tr:
            dk = d_gn
            d_c = [0,0,0]
            dc_norm=0

            # 探索ベクトルの大きさ＝信頼半径の場合に1となるフラグ
            tr_flag = 0

        # ガウス・ニュートンステップが信頼領域外の場合、コーシー点を計算
        else:

            #コーシー点の計算
            alpha0 = np.sum(nabla_S**2)/ np.sum( (Jf @ nabla_S)**2)
            d_c = -alpha0 * nabla_S

            #コーシー点の大きさと信頼半径を比較
            dc_norm = np.linalg.norm(d_c)

            # コーシー点が信頼領域内の場合
            if dc_norm <= r_tr:
                
                # 係数βの計算
                a_quadeq = np.sum((d_gn-d_c)**2)
                b_quadeq = np.dot(d_c,(d_gn-d_c))
                c_quadeq = dc_norm**2 -r_tr**2
                D_quadeq = b_quadeq **2 - a_quadeq * c_quadeq

                beta = (-b_quadeq+np.sqrt(D_quadeq))/a_quadeq

                if (beta<0) or (beta>1):
                    beta = (-b_quadeq-np.sqrt(D_quadeq))/a_quadeq

                # ドッグレッグステップの計算
                d_dl = beta * d_gn + (1-beta)*d_c
                ddl_norm = np.linalg.norm(d_dl)

                # 探索ベクトルをドッグレッグステップに設定
                dk = d_dl

            else:

                #コーシー点が信頼領域外の場合は、コーシー点と同じ方向で大きさ＝信頼半径の点を探索ベクトルとする
                dk = r_tr/dc_norm * d_c

            # 探索ベクトルの大きさ＝信頼半径の場合1となるフラグ （ガウス・ニュートンステップが信頼領域外の場合）
            tr_flag = 1 

        # ρk の計算　残差平方和S(a) の二次近似関数
        rho = (S(a+dk,x,y)-S(a,x,y))/(np.dot(nabla_S,dk) + 1/2*np.sum((Jf @ dk)**2))

        if rho < 0.25:
            r_tr = 0.25*r_tr

        elif (rho > 0.75) and (tr_flag == 1):
            r_tr = np.amin([2*r_tr,r_max])


        #パラメータの更新 評価関数ρkが閾値ηを超えた場合のみ更新
        if rho >= eta:
            a = a + dk
        
        #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)

        #途中経過の観察用 
        ''' 
        print(f'\n----{rep}回目計算途中結果-----\n'+\
              f'パラメータ[a,b,c]={a}\n'+\
              f'∇S(a,b,c)={nabla_S}\n'+ \
              f'二乗誤差:{sq_err}\n'+\
              f'ガウス・ニュートン点{d_gn}\n大きさ:{dgn_norm}\n' +\
              f'コーシー点{d_c}\n大きさ:{dc_norm}\n' +\
              f'信頼半径:{r_tr}\n'+\
              f'dk大きさ:{np.linalg.norm(dk)}\n'+\
              f'評価関数ρ:{rho}\n'+\
              f'--------------------------------')
        time.sleep(0.3)#'''

        # 反復回数の更新
        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キーで最初に戻ります。')

6.【補足】ドッグレッグステップとコーシー点の大小

本章では、ドッグレッグステップにおける二次近似関数 $S_{(2)}$ の値がコーシー点における $S_{(2)}$ の値以下となることを証明します。
$S_{(2)}$ を $\beta$ の関数 $g(\beta)$ として以下のように定義します。

$$\begin{aligned} g(\beta) &=\nabla S(\boldsymbol{a^{(k)}}) + \nabla S(\boldsymbol{a^{(k)}})^{T}\{\beta\boldsymbol{d_{gn}}+(1-\beta)\boldsymbol{d_{c}}\} +\dfrac{1}{2}\{\beta\boldsymbol{d_{gn}}+(1-\beta)\boldsymbol{d_{c}}\}^{T}J_{f}^{T}J_{f}\{\beta\boldsymbol{d_{gn}}+(1-\beta)\boldsymbol{d_{c}}\} \\ &= S_{(2)}(\boldsymbol{a^{(k)}}+\beta\boldsymbol{d_{gn}}+(1-\beta)\boldsymbol{d_{c}}) \hspace{25pt} (6.1) \end{aligned}$$

方針は $\beta_{2}=\beta_{1}+\Delta \beta\ (1 \gt \beta_{2} \gt \beta_{1} \gt 0)$ として、$g(\beta)$ を以下の形に変形します。

$$\begin{aligned} g(\beta_{2}) &=g(\beta_{1})+\delta \hspace{25pt} (6.2) \end{aligned}$$

$\delta$ が 0 以下であることを示せば、題意を満足することになります。

まず $g(\beta)$ を変形します。$(3.3.1)$ 式より

$$\begin{aligned} g(\beta) &=\nabla S(\boldsymbol{a^{(k)}}) - \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ \beta (J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta) I \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ & \hspace{15pt} +\dfrac{1}{2} \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ \beta (J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta) I \right\}^{T} J_{f}^{T}J_{f} \left\{ \beta (J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta) I \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ &=\nabla S(\boldsymbol{a^{(k)}}) - \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ \beta (J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta) I \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ & \hspace{15pt} +\dfrac{1}{2} \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ \beta (J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta) I \right\} \left\{ \beta I + \alpha_{0}(1-\beta) J_{f}^{T}J_{f} \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ &=\nabla S(\boldsymbol{a^{(k)}}) - \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ \beta (J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta) I \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ & \hspace{15pt} +\dfrac{1}{2} \nabla S(\boldsymbol{a^{(k)}})^{T} \left[ \beta^{2} (J_{f}^{T}J_{f})^{-1} + 2\alpha_{0}(1-\beta)\beta I +\{ \alpha_{0}(1-\beta)\}^{2} J_{f}^{T}J_{f} \right] \nabla S(\boldsymbol{a^{(k)}}) \hspace{25pt} (6.3)\\ \end{aligned}$$

よって $g(\beta_{2})=g(\beta_{1}+\Delta\beta)$ については

$$\begin{aligned} g(\beta_{2}) &=g(\beta_{1}+\Delta\beta) \\ &=\nabla S(\boldsymbol{a^{(k)}}) - \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ (\beta_{1}+\Delta\beta) (J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta_{1}-\Delta\beta) I \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ & \hspace{15pt} +\dfrac{1}{2} \nabla S(\boldsymbol{a^{(k)}})^{T} \left[ (\beta_{1}+\Delta\beta)^{2} (J_{f}^{T}J_{f})^{-1} + 2\alpha_{0}(1-\beta_{1}-\Delta\beta)(\beta_{1}+\Delta\beta) I +\{ \alpha_{0}(1-\beta_{1}-\Delta\beta)\}^{2} J_{f}^{T}J_{f} \right] \nabla S(\boldsymbol{a^{(k)}}) \\ &=\nabla S(\boldsymbol{a^{(k)}}) - \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ \beta_{1}(J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta_{1})I +\Delta\beta (J_{f}^{T}J_{f})^{-1}-\alpha_{0}\Delta\beta I \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ & \hspace{15pt} +\dfrac{1}{2} \nabla S(\boldsymbol{a^{(k)}})^{T} \Big[ \beta_{1}^{2} (J_{f}^{T}J_{f})^{-1} + 2\alpha_{0}(1-\beta_{1})\beta_{1}I +\{ \alpha_{0}(1-\beta_{1})\}^{2} J_{f}^{T}J_{f} \\ & \hspace{83pt} + (2\beta_{1}\Delta\beta+\Delta\beta^{2})(J_{f}^{T}J_{f})^{-1} +2\alpha_{0} \{-\Delta\beta(\beta_{1}+\Delta\beta)+(1-\beta_{1})\Delta\beta \}I \\ & \hspace{83pt} + \alpha_{0}^{2}\{-2\Delta\beta(1-\beta_{1})+ \Delta\beta^{2} \} J_{f}^{T}J_{f} \Big] \nabla S(\boldsymbol{a^{(k)}}) \\ &=\nabla S(\boldsymbol{a^{(k)}}) - \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ \beta_{1}(J_{f}^{T}J_{f})^{-1} + \alpha_{0}(1-\beta_{1})I \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ & \hspace{15pt} +\dfrac{1}{2} \nabla S(\boldsymbol{a^{(k)}})^{T} \Big[ \beta_{1}^{2} (J_{f}^{T}J_{f})^{-1} + 2\alpha_{0}(1-\beta_{1})\beta_{1}I +\{ \alpha_{0}(1-\beta_{1})\}^{2} J_{f}^{T}J_{f} \Big]\nabla S(\boldsymbol{a^{(k)}}) \\ & \hspace{15pt} +\dfrac{1}{2} \nabla S(\boldsymbol{a^{(k)}})^{T} \Big[ -2\Delta\beta (J_{f}^{T}J_{f})^{-1}+2\alpha_{0}\Delta\beta I +\Delta\beta(2\beta_{1}+\Delta\beta)(J_{f}^{T}J_{f})^{-1} +2\alpha_{0}\Delta\beta (1-2\beta_{1}-\Delta\beta )I \\ & \hspace{83pt} +\alpha_{0}^{2}\Delta\beta(2\beta_{1} + \Delta\beta -2) J_{f}^{T}J_{f} \Big] \nabla S(\boldsymbol{a^{(k)}}) \\ &=g(\beta_{1}) +\dfrac{1}{2} \nabla S(\boldsymbol{a^{(k)}})^{T} \Big[ \Delta\beta(2\beta_{1}+\Delta\beta-2)(J_{f}^{T}J_{f})^{-1} +2\alpha_{0}\Delta\beta \{2-2\beta_{1}-\Delta\beta\}I \\ & \hspace{105pt} -\alpha_{0}^{2}\Delta\beta(2-\beta_{1} - \beta_{2}) J_{f}^{T}J_{f} \Big] \nabla S(\boldsymbol{a^{(k)}}) \\ &=g(\beta_{1}) -\dfrac{1}{2} \nabla S(\boldsymbol{a^{(k)}})^{T} \Big[ \Delta\beta(2-\beta_{1} - \beta_{2})(J_{f}^{T}J_{f})^{-1} -2\alpha_{0}\Delta\beta \{2-\beta_{1} - \beta_{2}\}I \\ & \hspace{105pt} +\alpha_{0}^{2}\Delta\beta(2-\beta_{1} - \beta_{2}) J_{f}^{T}J_{f} \Big] \nabla S(\boldsymbol{a^{(k)}}) \\ &=g(\beta_{1}) -\dfrac{1}{2}\Delta\beta \{2-\beta_{1} - \beta_{2}\} \nabla S(\boldsymbol{a^{(k)}})^{T} \Big[ (J_{f}^{T}J_{f})^{-1} -2\alpha_{0}I +\alpha_{0}^{2}J_{f}^{T}J_{f} \Big] \nabla S(\boldsymbol{a^{(k)}}) \\ &=g(\beta_{1}) -\dfrac{1}{2}\Delta\beta \{2-\beta_{1} - \beta_{2}\} \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ (J_{f}^{T}J_{f})^{-1} -\alpha_{0}I \right\} \left\{I -\alpha_{0}J_{f}^{T}J_{f} \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ &=g(\beta_{1}) -\dfrac{1}{2}\Delta\beta \{2-\beta_{1} - \beta_{2}\} \nabla S(\boldsymbol{a^{(k)}})^{T} \left\{ (J_{f}^{T}J_{f})^{-1} -\alpha_{0}I \right\} J_{f}^{T}J_{f} \left\{(J_{f}^{T}J_{f})^{-1} -\alpha_{0}I \right\} \nabla S(\boldsymbol{a^{(k)}}) \\ &=g(\beta_{1}) -\dfrac{1}{2}\Delta\beta \{2-\beta_{1} - \beta_{2}\} \left[J_{f} \left\{(J_{f}^{T}J_{f})^{-1} -\alpha_{0}I \right\} \nabla S(\boldsymbol{a^{(k)}}) \right]^{T} \left[ J_{f} \left\{(J_{f}^{T}J_{f})^{-1} -\alpha_{0}I \right\} \nabla S(\boldsymbol{a^{(k)}}) \right]\\ &=g(\beta_{1}) -\dfrac{1}{2}\Delta\beta \{2-\beta_{1} - \beta_{2}\} \left|\left| J_{f} \left\{(J_{f}^{T}J_{f})^{-1} -\alpha_{0}I \right\} \nabla S(\boldsymbol{a^{(k)}}) \right|\right|^{2} \hspace{25pt} (6.4) \end{aligned}$$

となります。$0 \lt \beta_{1} \lt \beta_{2} \lt1$ より $2-\beta_{1}-\beta_{2} \gt 0$ であることから、 $g(\beta_{2})\leq g(\beta_{1})$ が成立します。以上より、$\beta$ が大きくなるほど二次近似関数 $S_{(2)}$ は小さくなり、コーシー点よりもドッグレッグステップの方が適していることが分かります。