ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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

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

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

  • gradient を使用した解析的勾配の計算

  • jacobian を使用した解析ヤコビアンの計算

  • matlabFunction を使用したシンボリックな結果の数値関数への変換

  • fplot を使用した解析の可視化

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

f(x)=i=1n/2100(x2i-x2i-12)2+(1-x2i-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-12+100y-x22
fsurf(f,[-2 2 -1 3])
view([-45 50])
colormap jet
caxis([0,1000])
colorbar

目的

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

  1. ローゼンブロック関数とヤコビアンの方程式の、テキスト形式の数式から数値コードへの変換。

  2. ヤコビアンの陽的な計算。複雑な方程式の場合、この作業は誤差が発生しやすくなります。

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

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

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

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 = 

(2x1-400x1x2-x12-2200x2-200x122x3-400x3x4-x32-2200x4-200x322x5-400x5x6-x52-2200x6-200x522x7-400x7x8-x72-2200x8-200x722x9-400x9x10-x92-2200x10-200x92)

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 = 

(1-x110x2-10x121-x310x4-10x321-x510x6-10x521-x710x8-10x721-x910x10-10x92)

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

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

JF = jacobian(F,x);

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

JF(1:10,1:10)
ans = 

(-1000000000-20x1100000000000-1000000000-20x3100000000000-1000000000-20x5100000000000-1000000000-20x7100000000000-1000000000-20x910)

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

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

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

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

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

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

FJFfun はリストで変数を受け取りますが、fsolve は変数のベクトルを要求します。関数 genericObj はベクトルをリストに変換します。無名関数 bananaObj は、genericObj の引数の値を確定するために定義されています。関数 genericObj のスパース パターンの使用方法に注目してください。また、このように FJFfungenericObj を共に使用するアプローチは、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)          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         102.771         1.5625           6.39            1.56
     6          7         102.771        3.90625           6.39            3.91
     7          8         87.7443       0.976563           2.19           0.977
     8          9         74.1426        2.44141           6.27            2.44
     9         10         74.1426        2.44141           6.27            2.44
    10         11          52.497       0.610352           1.52            0.61
    11         12         41.3297        1.52588           4.63            1.53
    12         13         34.5115        1.52588           6.97            1.53
    13         14         16.9716        1.52588           4.69            1.53
    14         15         8.16797        1.52588           3.77            1.53
    15         16         3.55178        1.52588           3.56            1.53
    16         17         1.38476        1.52588           3.31            1.53
    17         18        0.219553        1.16206           1.66            1.53
    18         19               0      0.0468565              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 でも解くことができます。 vpasolvefsolve の大きな違いは、fsolve が倍精度の演算で方程式を解くのに対して、vpasolve は可変精度の演算を使用するため、計算精度を制御できることです。

sol2 = vpasolve(F);

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

sol2.x1
ans = 1.0

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

for k=(1:64)
  sol2_array(k) = sol2.(char(x(k)));
end
sol2_array(5:25)
ans = (1.01.01.01.01.01.01.01.01.01.01.01.01.01.01.01.01.01.01.01.01.0)

補助関数

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

Copyright © 2016 The MathWorks, Inc.