最適化問題の応用例

Posted Nov 3, 2024 Updated Nov 3, 2024

By Guang-Ze Yang 27 min read

3.5 最適化問題の応用例

実際、これらの例の最終目標は、一見無関係に見える問題を最適化問題に変換する方法を示す。最適化問題をモデル化して解決するプロセスを通じて、元の問題の結果を見つけようとする。

3.5.1 線形回帰問題

Solutions of linear regression problems

ある関数の線形結合が次のように定義されていると仮定する。

\[g\left(x\right)=c_1 f_1 \left(x\right)+\ldotp \ldotp \ldotp +c_n f_n \left(x\right)\]

$f$ ：既知関数
$c$ ：未知係数

\[\textrm{Ac}=\mathit{\mathbf{y}}\]

ただし、 $x,y$ は観測データである。

したがって、最小二乗解は次のように求められる。

\[\mathit{\mathbf{c}}=\mathit{\mathbf{A}}\backslash \mathit{\mathbf{y}}\]

\\: 行列の左除算(left matrix divide) $=A^{-1} y$

例題 3.27

  
x=[0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5]';
y=[2.88;2.2576;1.9683;1.9258;2.0862;2.109;2.198;
2.541;2.9627;3.155;3.2052]; % samples input
A=[ones(size(x)) exp(-3*x),cos(-2*x).*exp(-4*x) x.^2];
c=A\y; c1=c' % least squares solution

c1 = 1x4    
1.2202    2.3390   -0.6792    0.8699

  
x0=(0:0.01:1.5)';
B=[ones(size(x0)) exp(-3*x0) cos(-2*x0).*exp(-4*x0) x0.^2];
y1=B*c; plot(x0,y1,x,y,'x') % draw fitting curve

3.5.2 最小二乗曲線近似

Least-squares curve fitting

3.5.1の線形回帰問題では、関数 f(x) は与えられた関数であると想定されていたため、その問題は線形方程式に変換できた。しかし、実際には、関数 f(x) には未定の係数が含まれることもある。したがって、線形方程式から解くことはできない。最適化問題に変換する必要がある。

定義 3.5

与えられたデータセット $\left(x_i ,y_i \right),i=1,2,\ldots,m$ が、あるモデル関数 $\hat{y} \left(x\right)=f\left(\mathit{\mathbf{a}},x\right)$ に基づいていると仮定する。ここで、 $\mathit{\mathbf{a}}$ は未定の係数のベクトルである。最小二乗法を用いた曲線近似の目的は、この係数を求めることにある。

この最適化問題を低レベルの方法（Low-level solution, 勾配降下法やニュートン法などの基本的なアルゴリズム）で解くことも可能であるが、より簡便に解く方法として、MATLABの関数 lsqcurvefit() を直接使用することもできる。lsqcurvefit() は最小二乗近似問題(Least Squares Fitting problem )に特化した関数で、与えられたモデル関数とデータに基づいて、係数ベクトル a の最適値を計算する。

Fun: MATLAB の関数または無名関数（anonymous function）を指定できる
a0: 未定の係数の初期値を含むベクトルで、初期推定値として与える。
x と y: 元の入力ベクトル xと出力ベクトル yである。多変数関数のフィッティングの場合、行列として指定できる。
options: Optimization Toolbox の設定である。

lsqcurvefit() 関数を使用すると、以下の出力が得られる。

a: 最適化された未定係数ベクトルで、データに最も適合する係数である。
Jm: 目的関数の値（残差平方和の最小値）で、フィッティングの精度を示す。
flag: 解法の成功を示し、正の値であれば最適化が成功して収束したことを示す。

アルゴリズム

By MATLAB

レーベンバーグ・マルカート法および信頼領域 Reflective 法は、fsolve にも使用される非線形最小二乗アルゴリズムに基づいています。

既定の信頼領域 Reflective 法アルゴリズムは部分空間の信頼領域法であり、[1] および [2] で説明する interior-reflective ニュートン法に基づいています。各反復は、前処理付き共役勾配 (PCG) 法を使用する大型線形システムの近似解を伴います。詳細については、信頼領域 Reflective 法の最小二乗を参照してください。
レーベンバーグ・マルカート法は [4]、[5]、[6] の参考文献で説明されています。詳細については、レーベンバーグ・マルカート法を参照してください。

'interior-point' アルゴリズムは fmincon 'interior-point' アルゴリズムにいくつかの変更を加えて使用します。詳細については、制約付き最小二乗法に合わせて変更された fmincon アルゴリズムを参照してください。

By ChatGPT:

1. 信頼領域 Reflective 法 ( Trust Region Reflective Algorithm )

アルゴリズムの概要:

信頼領域 Reflective 法は、部分空間の信頼領域法**をベースとし、**内点反射（interior\-reflective）ニュートン法を適用して非線形最小二乗問題を解きます。この方法は、大規模な最適化問題にも適応できるよう設計されており、fsolve関数にも使用されています。

初期化：

初期のパラメータ値 $x_0$ を設定し、初期信頼領域の半径 $Delta_0$ を設定します。

信頼領域内の解の探索：

各反復ステップで、目的関数の周辺を局所的に2次近似（通常は二次関数）し、信頼領域内でこの近似モデルの最適解を探索します。
このとき、制約がある場合には Reflective（反射） のアプローチを用いて制約範囲を外れないように解を調整します。

更新：

求めた解が実際の目的関数値の改善に寄与しているかを確認します。
もし改善が十分であれば、信頼領域の半径を拡大して次の反復へ進みます。
改善が不十分であれば、信頼領域の半径を縮小し、再度解を探索します。

収束条件の確認：

解が収束条件を満たした場合、アルゴリズムを終了します。収束条件には、勾配の大きさや信頼領域の半径の小ささなどが含まれます。

適用分野: 大規模で制約のある非線形最小二乗問題に適しています。

2. レーベンバーグ・マルカート法( Levenberg–Marquardt algorithm )

アルゴリズムの概要:

レーベンバーグ・マルカート法は、**勾配降下法とニュートン法**を組み合わせた手法で、目的関数の凹凸に応じて更新方向を調整します。このアプローチにより、非線形最小二乗問題の解を安定的に収束させることが可能です。

初期設定：

パラメータの初期値 x0 とダンピングパラメータ λ を設定します。
ダンピングパラメータ λ は初期状態では大きな値に設定され、勾配降下法に近い動作をします。

反復処理：

各ステップで、目的関数の近似を行い、更新されたパラメータ $x_{k+1}$ を次のように求めます。

ここで、

$J$ は残差ベクトル $r\left(x\right)$ のヤコビ行列
$r\left(x\right)$ は観測データと予測モデルの残差ベクトル
$I$ は単位行列
λ はダンピングパラメータ

このステップでは、更新が勾配降下法とニュートン法の間で切り替わります。

パラメータ λ の調整 ：

更新したパラメータ $x_{k+1}$ によって目的関数が改善された場合、つまり誤差が減少した場合は λ を小さくし、ニュートン法に近い動作を行うようにします。
逆に、改善が見られない場合は λ を大きくして、より保守的な勾配降下法に近い更新を行います。

収束判定：

目的関数の値やパラメータの変動がある基準を下回った場合に収束したとみなし、計算を終了します。

特徴:

勾配方向により安定性が向上し、収束速度が速くなります。特に小規模で滑らかな非線形問題に効果的です。

適用例:

MATLABのlsqcurvefitやlsqnonlinで、レーベンバーグ・マルカート法を使用した非線形最小二乗問題の解法が可能です。

3. 内点法（ interior-point ）

アルゴリズムの概要:

制約付き最適化のためのアルゴリズムであり、fminconのinterior-pointアルゴリズムを修正したものです。これにより、制約付き最小二乗問題にも対応できます。

初期化：

実行可能領域の内部にある初期点を選びます。
ペナルティパラメータやステップサイズなどのアルゴリズムに必要な初期設定を行います。

障害関数（バリア関数）の導入：

制約条件を考慮するために、目的関数に「障害関数」を加え、制約境界に近づくほど値が大きくなるように調整します。
例えば、制約付きの最適化問題を次のようにバリア項を追加した形に変換します。

ここで、 $g_{i\left(x\right)}$ は制約条件、 $\mu \;$ はペナルティパラメータです.

反復更新：

各反復ステップで、バリア関数を最適化しながら目的関数の最小化を進めます。
反復ごとにペナルティパラメータ μ を減少させ、障害関数の影響を弱めて制約境界に近づけます。

収束判定：

収束条件（例: 勾配が一定値以下になる、目的関数値の変動が小さくなるなど）が満たされたら反復を終了し、最適解を得ます。

適用分野: 制約がある非線形最小二乗問題に適しています。

例題 3.28

  
x=0:0.4:10; % generate samples
y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);

f=@(a,x)a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x);
a0=[1,1,1,1,1]; [xx,res]=lsqcurvefit(f,a0,x,y); % fitting

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>

xx

xx = 1x5    
    0.1200    0.2130    0.5400    0.1700    1.2300

  
x1=0:0.01:10;y1=f(xx,x1);plot(x1,y1,x,y,'o') % comparison

赤色点：データ
青色線：近似した曲線

実際、低レベルの方法で最適化問題を解くのは難しくない。新しい無名関数を定義することで、目的関数を表現し、これを使って最適化を実行できる。

  
F=@(a)norm(f(a,x)-y); x1=fminunc(F,a0)

Solver stopped prematurely.

fminunc stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 5.000000e+02.
x1 = 1x5    
    0.1200    0.2130    0.5400    0.1700    1.2300

例題 3.29

$y\left(x\right)=\textrm{ax}+{\textrm{bx}}^2 e^{-\textrm{cx}} +d$

  
x=0.1:0.1:1; % input samples
y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090, 4.2147,4.5191,4.8232,5.1275];

  
f=@(a,x)a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);
a=[1;2;2;3]; % prototype function and initial values
while (1)
    [a,f0,cc,flag]=lsqcurvefit(f,a,x,y); 
    if flag>0 
        break;
    end
end % least squares solutions

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>

  
y1=f(a,x); plot(x,y,x,y1,'o')

例題 3.30

独立変数 x, y, z がデータである。このデータを使って最小二乗法を用いて係数 $a_i$ を求める。

  
f=@(a,X)a(1)*X(:,1).^(a(2)*X(:,1))+...
a(3)*X(:,2).^(a(4)*(X(:,1)+X(:,2)))+...
a(5)*X(:,3).^(a(6)*(X(:,1)+X(:,2)+X(:,3))); % prototype
XX=load('c3data1.dat'); X=XX(:,1:3); v=XX(:,4); % load samples
a0=[2 3 2 1 2 3]; 
[a,f0,err,flag] = lsqcurvefit(f,a0,X,v)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>
a = 1x6    
    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000

f0 = 1.0904e-07
err = 150x1    
1.0e-04 *

   -0.3596
    0.1326
    0.1289
   -0.3689
   -0.0006
    0.3872
   -0.2540
    0.3401
   -0.1031
    0.1154

flag = 1

3.5.3 微分方程式の境界値問題におけるシューティング法

Shooting method in boundary value differential equations

シューティング法を使った最適化の利用例として、微分方程式の問題を最適化問題に変換する方法を示す。常微分方程式（ODE, ordinary differential equations）の初期値問題は、数値解析で広く研究されている。微分方程式の数学モデルが既知であれば、初期値問題を解くための手法が用意されている。初期値問題は、次のように数学的に記述される。

ここで、 $x\prime \left(t\right)$ は x の t に関する導関数、 $f\left(t,x\right)$ は x と t に依存する関数である。このような形の初期値問題を解くためには、MATLAB の ode45 関数がよく使用される。

実際の応用では、初期時刻 $t_0$ でのいくつかの成分 $x\left(t_0 \right)$ と、ある時刻 $t_0$ でのいくつかの成分 $x\left(t_n \right)$ が与えられている場合に、このような微分方程式を解くには、シューティング法(shooting methodshooting method)がよく使われる。

> By ChatGPT

シューティング法

未知の初期値を含む $x\left(t_0 \right)$ に対し適切な値を割り当てることで、元の境界値問題(BVP: Boundary Value Problem)を初期値問題(IVP: Initial Value Problem)に変換する。

シューティング法の手順

未知の初期値の仮定: 初期時刻 $t_0$ における未知の初期値 $x\left(t_0 \right)$ を仮定する。この仮定により、境界値問題を初期値問題として取り扱えるようになる。
初期値問題の数値解法: ode45() などの数値解法を用いて、初期値問題を解き、方程式の最終時刻 $t_n$ における近似解 $\hat{x} \left(t_n \right)$ を得る。
誤差の測定: $\hat{x} \left(t_n \right)$ と与えられた境界条件 $x\left(t_n \right)$ との間の誤差を計算する。この誤差は、目的関数として最小化される。
初期値の調整: 得られた誤差に基づき、初期値 $x\left(t_0 \right)$ を調整し、再び微分方程式を解く。このプロセスを、目標の境界条件に達するまで繰り返す。

このループ構造により、問題は初期値の最適化問題に変換され、最適化アルゴリズムで境界条件に適合する初期値が求められる。

初期値問題と境界値問題の違いまとめ: 微分方程式を解くため


特徴	初期値問題 (IVP)	境界値問題 (BVP)
条件の指定方法	初期時刻の値	範囲の両端（境界）の値
解の一意性	通常一意	一意とは限らない
数値解法の例	Runge-Kutta法, Euler法	シューティング法, 有限差分法

例題 3.31

境界条件 $x_1 \left(0\right)=1,x_2 \left(0\right)=2,x_3 \left(10\right)=-0\ldotp 021677,x_4 \left(10\right)=0\ldotp 15797$ が与えられる。 $t\in \left(0,10\right)$ 内の常微分方程式を解く。

まず、変数 $x_3 \left(0\right),x_4 \left(0\right)$ を決定係数として、元の問題を最適化問題に変換する。

$y_1 =x_3 \left(0\right),y_2 =x_4 \left(0\right)$ とすると、以下の式が得られる。

  
function z=c3mode(y,f,x0,xn)
x0=[x0(1:2); y]; 
[t,x]=ode45(f,[0,10],x0);
z=norm([x(end,3:4)]'-xn);
end

f=@(t,x)[x(2); -x(1)-3*x(2)+exp(-5*t); x(4); 2*x(1)-3*x(2)-3*x(3)-4*x(4)-sin(t)];
x0=[1; 2]; 
xn=[-0.021677; 0.15797];
g=@(x)c3mode(x,f,x0,xn); % convert the problem with extra parameters
x2=rand(2,1); %初期値ランダム
x3=fminunc(g,x2) % solve optimization problem

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>
x3 = 2x1    
   -0.7093
    0.0757

  
[t,x]=ode45(f,[0,10],[x0; x3]); 
plot(t,x), legend("x_1","x_2","x_3","x_4")

  
x(end,:)

ans = 1x4    
    0.0474   -0.0181   -0.0217    0.1580

  
[t,x]=ode45(f,[0,10],[x0; -1; 2]); %境界条件が異なる場合
plot(t,x), legend("x_1","x_2","x_3","x_4")

  
x(end,:)

ans = 1x4    
    0.0474   -0.0181   -0.0217    0.1579

ode45: 初期値で微分方程式を数値的に解き、値を求める。
fminuc: 目的関数 g(x) = c3mode(x, f, x0, xn) を最小化し、未知の初期条件 x2 を最適化する。

3.5.4 代数方程式を最適化問題への変換

Converting algebraic equations into optimization problems

非線形代数方程式の数値解法は、第2章で紹介された。実際には、方程式の解を求める問題も無制約最適化問題に変換できる。方程式 $f\left(\mathit{\mathbf{x}}\right)=0$ が与えられているとする。このとき、 $\mathit{\mathbf{x}}$ が満たすべき条件は $f_i \left(\mathit{\mathbf{x}}\right)=0$ でなければならない。そこで、問題を最適化問題に変換でき、目的関数として各方程式の平方和を取ることができる。このようにして目的関数を最小化する x を決定ベクトルとする：

例題 3.32

代数方程式の問題を制約なし最適化問題に変換する：

無名関数で代数方程式を定義し、fimplicit()で陰関数をプロットする。

  
f1=@(x1,x2)4*x1.^3+4*x1.*x2+2*x2.^2-42*x1-14;
f2=@(x1,x2)4*x2.^3+2*x1.^2+4*x1.*x2-26*x2-22;
fimplicit(f1), hold on, fimplicit(f2), hold off, legend("f_1","f_2"), xlabel("x_1"),ylabel("x_2")

曲線の交点が代数方程式の解であり、たくさんあるので、最適化問題で解けると、初期値により結果が変わる。

  
f=@(x)(f1(x(1),x(2)))^2+(f2(x(1),x(2)))^2;
[x,f0]=fminunc(f,rand(2,1))

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>
x = 2x1    
   -0.2708
   -0.9230

f0 = 8.6138e-12

まとめ

たくさん問題は、最適化問題に変換することができる。今回では以下の問題を取り扱った

線形回帰問題: 予測値と実際のデータ値の差の平方和を最小化する
最小二乗曲線近似: 最小二乗法を用いて誤差の平方和を最小化する
微分方程式の境界値問題: 初期値問題に変換し、境界条件を満たすようにパラメータを調整するため、誤差を最小化する
代数方程式の問題: 目的関数として方程式の誤差の平方和を最小化する。

また、人工知能のニューラルネットワークでは、Loss関数を最小化するために、最適化問題として設定される。勾配降下法や確率的勾配降下法などが用いられる。強化学習では、エージェントが環境と相互作用して報酬を最大化する方策を学ぶ際に、報酬関数の期待値を最大化する。ハイパーパラメータのチューニングでは、モデルの性能を最適化するために、ハイパーパラメータの最適な組み合わせを探索する最適化が行われる。

参考資料

Math, Optimization

This post is licensed under CC BY 4.0 by the author.

3.5 最適化問題の応用例

3.5.1 線形回帰問題

例題 3.27

3.5.2 最小二乗曲線近似

定義 3.5

アルゴリズム

1. 信頼領域 Reflective 法 ( Trust Region Reflective Algorithm )

2. レーベンバーグ・マルカート法( Levenberg–Marquardt algorithm )

3. 内点法（ interior-point ）

例題 3.28

例題 3.29

例題 3.30

3.5.3 微分方程式の境界値問題におけるシューティング法

シューティング法

例題 3.31

3.5.4 代数方程式を最適化問題への変換

例題 3.32

まとめ

参考資料

Trending Tags