このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
最適化の精度とパフォーマンスの改善
この例では、非線形方程式系を解く場合に、Symbolic Math Toolbox™ が誤差の最小化に役立つことを示します。
この例では、次の Symbolic Math Toolbox 機能を使用します。
gradient
を使用した解析的勾配の計算jacobian
を使用した解析ヤコビアンの計算matlabFunction
を使用したシンボリックな結果の数値関数への変換fplot
を使用した解析の可視化
目標は、ローゼンブロック関数 (Rosenbrock function)、いわゆる "バナナ" 関数の最小値を求めることです。
fplot
を使用して、2 変数のローゼンブロック関数を即座に可視化します。
syms x y a=1; b=100; f(x,y)=(a-x)^2+b*(y-x^2)^2
f(x, y) =
fsurf(f,[-2 2 -1 3])
view([-45 50])
colormap jet
colorbar
目的
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 =
ヤコビアンを使用した場合と使用しない場合の非線形系の解法 (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 =
方程式系を表す行列のヤコビアンの計算
jacobian
を使用して、F
のヤコビアンを計算します。この関数は、ヤコビアンをシンボリックに計算するため、導関数の数値近似に伴う誤差を防ぎます。
JF = jacobian(F,x);
ヤコビ行列の初めの 10 行 10 列を表示します。
JF(1:10,1:10)
ans =
ヤコビ行列 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 =
sol2
をシンボリック ベクトルに変換した後は、一定の範囲の解にアクセスすることもできます。たとえば、10 番目から 20 番目までの解を表示します。
for k=(1:64) sol2_array(k) = sol2.(char(x(k))); end sol2_array(10:20)
ans =
補助関数
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