MATLAB: Comparing accuracy of symbolic computation with matlabFunction

accuracyfminconMATLABmatlabfunctionOptimization ToolboxsymbolicSymbolic Math Toolbox

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.

Best Answer

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.
Related Question