最適化の精度とパフォーマンスの改善

この例では、非線形方程式系を解く場合に、Symbolic Math Toolbox™ が誤差の最小化に役立つことを示します。

この例では、次の Symbolic Math Toolbox 機能を使用します。

gradientを使用した解析的勾配の計算
jacobianを使用した解析ヤコビアンの計算
matlabFunctionを使用したシンボリックな結果の数値関数への変換
fplot を使用した解析の可視化

目標は、ローゼンブロック関数 (Rosenbrock function)、いわゆる "バナナ" 関数の最小値を求めることです。

$f (x) = \sum_{i = 1}^{n / 2} 100 (x_{2 i} - x_{2 i - 1}^{2})^{2} + (1 - x_{2 i - 1})^{2}$

fplot を使用して、2 変数のローゼンブロック関数を即座に可視化します。

syms x y
a=1; b=100;
f(x,y)=(a-x)^2+b*(y-x^2)^2

f(x, y) =  ${(x - 1)}^{2} + 100 {(y - x^{2})}^{2}$

fsurf(f,[-2 2 -1 3])
view([-45 50])
colormap jet
colorbar

Figure contains an axes object. The axes object contains an object of type functionsurface.

目的

Optimization Toolbox™ のヤコビアンを使用した場合と使用しない場合の非線形系の解法 (Optimization Toolbox)の例では、関数fsolve (Optimization Toolbox)を使用して同じ問題を解いています。その例に示されたワークフローには、潜在的誤差要因が 2 つあります。次が必要となります。

ローゼンブロック関数とヤコビアンの方程式の、テキスト形式の数式から数値コードへの変換。
ヤコビアンの陽的な計算。複雑な方程式の場合、この作業は誤差が発生しやすくなります。

この例では、シンボリック数学を使用して問題のステートメントを記述し、以降の手順でこのような誤差を最小化する、あるいは取り除くことを説明します。

ローゼンブロック関数の非線形方程式系への書き換え

最初に、ローゼンブロック関数 f を非線形方程式系 F に変換します。ここで、F は f の勾配です。

clear
n = 64;
x = sym('x',[n,1]);
f = sum(100*(x(2:2:n)-x(1:2:n).^2).^2 + (1-x(1:2:n)).^2);
F = gradient(f,x);

最初の 10 個の方程式を表示します。

F(1:10)

ans = 
 $(\begin{array}{c} 2 x_{1} - 400 x_{1} (x_{2} - {x_{1}}^{2}) - 2 \\ 200 x_{2} - 200 {x_{1}}^{2} \\ 2 x_{3} - 400 x_{3} (x_{4} - {x_{3}}^{2}) - 2 \\ 200 x_{4} - 200 {x_{3}}^{2} \\ 2 x_{5} - 400 x_{5} (x_{6} - {x_{5}}^{2}) - 2 \\ 200 x_{6} - 200 {x_{5}}^{2} \\ 2 x_{7} - 400 x_{7} (x_{8} - {x_{7}}^{2}) - 2 \\ 200 x_{8} - 200 {x_{7}}^{2} \\ 2 x_{9} - 400 x_{9} (x_{10} - {x_{9}}^{2}) - 2 \\ 200 x_{10} - 200 {x_{9}}^{2} \end{array})$

ヤコビアンを使用した場合と使用しない場合の非線形系の解法 (Optimization Toolbox)に示された中間結果と一致させるため、同じ形式の F を使用します。

F(1:2:n) = simplify(F(1:2:n) + 2*x(1:2:n).*F(2:2:n));
F(1:2:n) = -F(1:2:n)/2;
F(2:2:n) = F(2:2:n)/20;

これで方程式系は同一です。

F(1:10)

ans = 
 $(\begin{array}{c} 1 - x_{1} \\ 10 x_{2} - 10 {x_{1}}^{2} \\ 1 - x_{3} \\ 10 x_{4} - 10 {x_{3}}^{2} \\ 1 - x_{5} \\ 10 x_{6} - 10 {x_{5}}^{2} \\ 1 - x_{7} \\ 10 x_{8} - 10 {x_{7}}^{2} \\ 1 - x_{9} \\ 10 x_{10} - 10 {x_{9}}^{2} \end{array})$

方程式系を表す行列のヤコビアンの計算

jacobianを使用して、F のヤコビアンを計算します。この関数は、ヤコビアンをシンボリックに計算するため、導関数の数値近似に伴う誤差を防ぎます。

JF = jacobian(F,x);

ヤコビ行列の初めの 10 行 10 列を表示します。

JF(1:10,1:10)

ans = 
 $(\begin{array}{c} - 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ - 20 x_{1} & 10 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - 20 x_{3} & 10 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - 20 x_{5} & 10 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 20 x_{7} & 10 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 20 x_{9} & 10 \end{array})$

ヤコビ行列 JF のほとんどの要素は 0 です。ただし、この行列をmatlabFunctionに変換する場合、結果は数値密行列になります。多くの場合、スパース行列での演算は、密行列での同じ演算よりも効率的です。

したがって、効率的なコードを作成するには、JF の非ゼロの要素のみをシンボリックベクトルに変換してから、matlabFunction を呼び出します。is および js は JF のスパースパターンを表します。

[is,js] = find(JF);
JF = JF(JF~=0);

方程式系とヤコビアンの MATLAB 関数への変換

ローゼンブロック関数を表す方程式系 F は、シンボリック式で構成されるシンボリック行列です。これを関数 fsolve で解くには、この系をmatlabFunctionに変換します。この手順では、F とヤコビアン JF の両方を、1 つのファイルベースの MATLAB 関数 FJFfun に変換すると便利です。基本的には、これにより F と JF は変数を再利用できるようになります。または、個々の MATLAB 関数に変換することもできます。

matlabFunction([F;JF],'var',x,'file','FJFfun');

FJFfun はリストで変数を受け取りますが、fsolve は変数のベクトルを要求します。関数 genericObj はベクトルをリストに変換します。無名関数 bananaObj は、genericObj の引数の値を確定するために定義されています。関数 genericObj のスパースパターンの使用方法に注目してください。また、このように FJFfun と genericObj を共に使用するアプローチは、F で表される特定の式に限らず、任意の F にそのまま使用できます。

bananaobj = @(y) genericObj(y,@FJFfun,is,js)

bananaobj = function_handle with value:
    @(y)genericObj(y,@FJFfun,is,js)

`fsolve` を使用した非線形方程式系の数値的求解

MATLAB 関数に変換された方程式系に fsolve を使用します。初期の点として、奇数インデックスの場合は x(i) = -1.9、偶数インデックスの場合は x(i) = 2 を使用します。ソルバーの進行状況を確認するため、'Display' を 'iter' に設定します。 bananaobj に定義されたヤコビアンを使用するために、'Jacobian' を 'on' に設定します。

x0(1:n,1) = -1.9;
x0(2:2:n,1) = 2;
options = optimoptions(@fsolve,'Display','iter','Jacobian','on');
[sol1,Fsol,exitflag,output,JAC] = fsolve(bananaobj,x0,options);

                                             Norm of      First-order   Trust-region
 Iteration  Func-count     ||f(x)||^2           step       optimality         radius
     0          1             8563.84                             615              1
     1          2             3093.71              1              329              1
     2          3             225.104            2.5             34.8            2.5
     3          4              212.48           6.25             34.1           6.25
     4          5              212.48           6.25             34.1           6.25
     5          6              212.48         1.5625             34.1           1.56
     6          7             116.793       0.390625             5.79          0.391
     7          8             109.947       0.390625            0.753          0.391
     8          9             99.4696       0.976563              1.2          0.977
     9         10             83.6416        2.44141             7.13           2.44
    10         11             77.7663        2.44141             9.94           2.44
    11         12             77.7663        2.44141             9.94           2.44
    12         13              43.013       0.610352             1.38           0.61
    13         14             36.4334       0.610352             1.58           0.61
    14         15             34.1448        1.52588             6.71           1.53
    15         16             18.0108        1.52588             4.91           1.53
    16         17             8.48336        1.52588             3.74           1.53
    17         18             3.74566        1.52588             3.58           1.53
    18         19             1.46166        1.52588             3.32           1.53
    19         20             0.29997        1.24265             1.94           1.53
    20         21                   0      0.0547695                0           1.53

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

vpasolve を使用した非線形方程式系の数値的求解

非線形方程式系 F は、Symbolic Math Toolbox で提供される数値ソルバーの関数 vpasolve でも解くことができます。 vpasolve と fsolve の大きな違いは、fsolve が倍精度の演算で方程式を解くのに対して、vpasolve は可変精度の演算を使用するため、計算精度を制御できることです。

sol2 = vpasolve(F);

方程式系の解を単一の変数 sol2 に代入すると、vpasolve は構造体配列の形式で解を返します。個々の解 (構造体配列の各フィールド) にアクセスできます。

sol2.x1

ans =  $1.0$

sol2 をシンボリックベクトルに変換した後は、一定の範囲の解にアクセスすることもできます。たとえば、10 番目から 20 番目までの解を表示します。

for k=(1:64)
  sol2_array(k) = sol2.(char(x(k)));
end
sol2_array(10:20)

ans =  $(\begin{array}{c} 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 \end{array})$

補助関数

function [F,J] = genericObj(x,FJFfun,is,js)
% genericObj takes as input, vector x, Function and Jacobian evaluation 
% function handle FJFfun, and sparsity pattern is,js and returns as output,
% the numeric values of the Function and Jacobian: F and Jfunction 
% FJs(1:length(x)) is the numeric vector for the Function
% FJs(length(x)+1:end) is the numeric vector for the non-zero entries of the Jacobian

xcell = num2cell(x);
FJs = FJFfun(xcell{:});
F = FJs(1:length(x));
J = sparse(is,js,FJs(length(x)+1:end));

end