Comparing accuracy of symbolic computation with matlabFunction
显示 更早的评论
Hi,
I have an expression called as DCritn (dynamically generated based on user input) which has 2 parameters and 2 decision variables.
I'm trying to minimize this expression using fmincon for different combinations of parameters (and different start points)
I have created a handle for the expression using matlabFunction and it works fine (well for most part).
For certain combination of parameters and start points, the function handle evaluation does not match the symbolic evaluation of the same expression by using double(subs(expression)). Why? Which is more correct?
%Generating the list of decision/explanatory variables;
t1 = sym(zeros(1, 2));
for k=1:2
t1(k) = sym(sprintf('x%d', k));
end
%Generating the list of parameters;
q = sym(zeros(1, 2));
for k=1:2
q(k) = sym(sprintf('q%d', k));
end
%Below is the value for expression/objective function DCritn. Since it is dynamically generated so providing here an example;
%Note, for below example, "simplify" function could be used to make it shorter. Consciously not using it, as in my real-life example I might have 6 or even 10 decision variables with even longer expressions for which "simplify" takes forever to run and never able to shorten it;
DCritn = ((q1 - q2)^4*(q1^2*exp(2*q1*x1)*exp(2*q1*x2)*exp(2*q2*x1)*exp(2*q2*x2) + q2^2*exp(2*q1*x1)*exp(2*q1*x2)*exp(2*q2*x1)*exp(2*q2*x2) - 2*q1*q2*exp(2*q1*x1)*exp(2*q1*x2)*exp(2*q2*x1)*exp(2*q2*x2)))/(q1^2*(q1^2*x1^2*exp(2*q1*x2)*exp(2*q2*x1) + q2^2*x1^2*exp(2*q1*x1)*exp(2*q1*x2) + q1^2*x1^2*exp(2*q2*x1)*exp(2*q2*x2) + q1^2*x2^2*exp(2*q1*x1)*exp(2*q2*x2) + q2^2*x1^2*exp(2*q1*x1)*exp(2*q2*x2) + q2^2*x2^2*exp(2*q1*x1)*exp(2*q1*x2) + q1^2*x2^2*exp(2*q2*x1)*exp(2*q2*x2) + q2^2*x2^2*exp(2*q1*x2)*exp(2*q2*x1) - 2*q1^2*x1^2*exp(q1*x2)*exp(q2*x2)*exp(2*q2*x1) - 2*q1^2*x2^2*exp(q1*x1)*exp(q2*x1)*exp(2*q2*x2) - 2*q2^2*x1^2*exp(q1*x2)*exp(q2*x2)*exp(2*q1*x1) - 2*q2^2*x2^2*exp(q1*x1)*exp(q2*x1)*exp(2*q1*x2) + q1^4*x1^2*x2^2*exp(2*q1*x1)*exp(2*q2*x2) + q1^4*x1^2*x2^2*exp(2*q1*x2)*exp(2*q2*x1) - 2*q2^2*x1*x2*exp(2*q1*x1)*exp(2*q1*x2) - 2*q1^2*x1*x2*exp(2*q2*x1)*exp(2*q2*x2) - 2*q1^3*x1*x2^2*exp(2*q1*x1)*exp(2*q2*x2) - 2*q1^3*x1^2*x2*exp(2*q1*x2)*exp(2*q2*x1) + q1^2*q2^2*x1^2*x2^2*exp(2*q1*x1)*exp(2*q2*x2) + q1^2*q2^2*x1^2*x2^2*exp(2*q1*x2)*exp(2*q2*x1) - 2*q1*q2*x1^2*exp(q1*x1)*exp(q2*x1)*exp(2*q1*x2) - 2*q1*q2*x1^2*exp(q1*x1)*exp(q2*x1)*exp(2*q2*x2) - 2*q1*q2*x2^2*exp(q1*x2)*exp(q2*x2)*exp(2*q1*x1) - 2*q1*q2*x2^2*exp(q1*x2)*exp(q2*x2)*exp(2*q2*x1) + 2*q1^2*x1*x2*exp(q1*x1)*exp(q2*x1)*exp(2*q2*x2) + 2*q2^2*x1*x2*exp(q1*x1)*exp(q2*x1)*exp(2*q1*x2) + 2*q1^2*x1*x2*exp(q1*x2)*exp(q2*x2)*exp(2*q2*x1) + 2*q2^2*x1*x2*exp(q1*x2)*exp(q2*x2)*exp(2*q1*x1) - 2*q1*q2^2*x1*x2^2*exp(2*q1*x2)*exp(2*q2*x1) - 2*q1*q2^2*x1^2*x2*exp(2*q1*x1)*exp(2*q2*x2) + 2*q1^2*q2*x1*x2^2*exp(2*q1*x1)*exp(2*q2*x2) + 2*q1^2*q2*x1*x2^2*exp(2*q1*x2)*exp(2*q2*x1) + 2*q1^2*q2*x1^2*x2*exp(2*q1*x1)*exp(2*q2*x2) + 2*q1^2*q2*x1^2*x2*exp(2*q1*x2)*exp(2*q2*x1) + 2*q1^3*x1*x2^2*exp(q1*x1)*exp(q2*x1)*exp(2*q2*x2) - 2*q1^3*x1^2*x2*exp(q1*x1)*exp(q2*x1)*exp(2*q2*x2) - 2*q1^3*x1*x2^2*exp(q1*x2)*exp(q2*x2)*exp(2*q2*x1) + 2*q1^3*x1^2*x2*exp(q1*x2)*exp(q2*x2)*exp(2*q2*x1) - 2*q1*q2*x1*x2*exp(2*q1*x1)*exp(2*q2*x2) - 2*q1*q2*x1*x2*exp(2*q1*x2)*exp(2*q2*x1) - 2*q1^3*q2*x1^2*x2^2*exp(2*q1*x1)*exp(2*q2*x2) - 2*q1^3*q2*x1^2*x2^2*exp(2*q1*x2)*exp(2*q2*x1) + 2*q1^3*x1*x2^2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) + 2*q1^3*x1^2*x2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) + 2*q1*q2*x1*x2*exp(q1*x1)*exp(q2*x1)*exp(2*q1*x2) + 2*q1*q2*x1*x2*exp(q1*x1)*exp(q2*x1)*exp(2*q2*x2) + 2*q1*q2*x1*x2*exp(q1*x2)*exp(q2*x2)*exp(2*q1*x1) + 2*q1*q2*x1*x2*exp(q1*x2)*exp(q2*x2)*exp(2*q2*x1) - 2*q1^4*x1^2*x2^2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) + 4*q1*q2*x1^2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) + 4*q1*q2*x2^2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) - 2*q1^2*x1*x2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) - 2*q2^2*x1*x2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) + 2*q1*q2^2*x1*x2^2*exp(q1*x1)*exp(q2*x1)*exp(2*q1*x2) - 2*q1*q2^2*x1^2*x2*exp(q1*x1)*exp(q2*x1)*exp(2*q1*x2) - 2*q1^2*q2*x1*x2^2*exp(q1*x1)*exp(q2*x1)*exp(2*q1*x2) + 2*q1^2*q2*x1^2*x2*exp(q1*x1)*exp(q2*x1)*exp(2*q1*x2) - 2*q1*q2^2*x1*x2^2*exp(q1*x2)*exp(q2*x2)*exp(2*q1*x1) + 2*q1*q2^2*x1^2*x2*exp(q1*x2)*exp(q2*x2)*exp(2*q1*x1) - 2*q1^2*q2*x1*x2^2*exp(q1*x1)*exp(q2*x1)*exp(2*q2*x2) + 2*q1^2*q2*x1*x2^2*exp(q1*x2)*exp(q2*x2)*exp(2*q1*x1) + 2*q1^2*q2*x1^2*x2*exp(q1*x1)*exp(q2*x1)*exp(2*q2*x2) - 2*q1^2*q2*x1^2*x2*exp(q1*x2)*exp(q2*x2)*exp(2*q1*x1) + 2*q1^2*q2*x1*x2^2*exp(q1*x2)*exp(q2*x2)*exp(2*q2*x1) - 2*q1^2*q2*x1^2*x2*exp(q1*x2)*exp(q2*x2)*exp(2*q2*x1) + 4*q1^3*q2*x1^2*x2^2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) - 2*q1^2*q2^2*x1^2*x2^2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) + 2*q1*q2^2*x1*x2^2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) + 2*q1*q2^2*x1^2*x2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) - 4*q1^2*q2*x1*x2^2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) - 4*q1^2*q2*x1^2*x2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2) - 4*q1*q2*x1*x2*exp(q1*x1)*exp(q1*x2)*exp(q2*x1)*exp(q2*x2)));
%Evaluating gradient of objective function for supplying to fmincon;
DCritnGrad= jacobian(DCritn ,t1).';
%Defining matlab function handle for objective func;
DCritnF = matlabFunction(DCritn,'vars',{t1.',q.'},'file',FileNmForDCritn);
%Defining matlab function handles for objective func and gradient together;
%This will be useful when fmincon needs both and saves time in evaluating them
%independently;
DCritnObjFObjGradF = matlabFunction(DCritn,DCritnGrad,'vars',{t1.',q.'},'file',FileNmForDCritnObj_Grad);
Now I try to evaluate the value of DCritn in three different ways. As seen below the answers don't match. The differences might seem less but since I'm running fmincon for many different parameter combinations and different start points, wrong results get returned for considerable instances.
>> DCritnF([.0001,85].',[0.24875,0.07].')
ans =
1.206125170348200e+09
>> DCritnObjFObjGradF([.0001,85].',[0.24875,0.07].')
ans =
1.206125876160960e+09
>> double(subs(DCritn,[t1,q],{[.0001,85,0.24875,0.07]}))
ans =
1.206125421571699e+09
The problem is acute as for some start points the fmincon stops with 1 iteration with exitflag -3. On the other hand if I had simplified DCritn using "simplify" function and then fed to matlabFunction and then run fmincon with same start point as before it terminates with exitflag 1 which is good.
So I guess it comes down to accuracy of function evaluation.
采纳的回答
Christopher Berry
2014-8-8
1 个投票
Hari,
It seems you have answered your own question here, but just to clarify on the issue of accuracy, I will add a little to it.
Functions subs(s) vs matlabFunction(s) Accuracy
Using the functions|double(subs(s))| will definitely give you a more accurate evaluation of the symbolic expression, s, than matlabFunction(s) will.
This is because matlabFunction(s) breaks the symbolic expression s up into pieces (just look at your FileNmForDCritn for proof of this) which are evaluated and stored into double precision variables. Then the complete expression is evaluated (summed,multiplied,etc..) using these stored values. Each evaluation and storage is a chance for floating point error to occur and accumulate with subsequent steps.
The double(subs(s)) method of evaluating s, however, does not convert to double precision until the final step. Each intermediary step is computed with a "working" precision specified by digits (just like the function vpa). The default value for the the working precision is 32.
This means that every exp calculated in your symbolic expression is calculated to 32 digits of accuracy, roughly twice as much as when using a double.
Functions matlabFunction(s1) vs matlabFunction(s1,s2) Accuracy
This difference is harder to predict and compare, since it does not seem like there should be any reason these two methods of evaluating s1 arrive at different values. When you look inside the two function created (like your FileNmForDCritn and FileNmForDCritnObj_Grad), however, you will probably see vastly different expressions for s1 in each file. Before breaking the expressions up, matlabFunction first tries to simplify and optimize the pieces it stores as doubles taking into consideration both outputs the function will be producing simultaneously, not one at a time. This means that the order of operations can be different for each function.
Again, what it ultimately boils down to is that floating point arithmetic does not necessarily satisfy associative, distributive or commutative laws, so how you evaluate an expression matters.
更多回答（1 个）
Christopher Creutzig
2014-8-29
Small correction to Christopher Berry's answer: double(subs(s)) does not compute to 32 digits (unless the values substituted already have been created with the vpa function), it computes exact values – for floating point input, exact rationals:
>> syms x, s = x^100; x = 0.123; subs(s)
ans =
97838805977257474352566705351629014033137938449734350966526074342064414099511156930426773522415958061389200997320437636836296142253482249885877442849062074323416253749444792245426920843456133929113701176246001/1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
To get better performance, thus, explicitly substitute with vpa values:
>> syms x, s = x^100; subs(s, x, vpa(0.123))
ans =
9.7838805977257474352566705351629e-92
Of course, that implies all the computations take place with floating-point numbers and will accumulate floating-point errors – just with higher precision (which you can set as high or low as you want using the digits function), which often leads to higher accuracy.
类别
在 帮助中心 和 File Exchange 中查找有关 Conversion Between Symbolic and Numeric 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
