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.

[Hari]

I realize now (after experimenting with code and reading it on matlab site) that double(subs(symbolic expression)) is more accurate (less chances of floating point issues) than using matlabFunction/anonymous function. But then, former is very very slow while later is blazing fast. For my problem, fmincon runs 20 times faster (on average) when using matlabFunction/anonymous function approach as opposed to using double(subs)).

What I decided to do is use matlabFunction approach as my primary input for fmincon and re-run fmincon with double(subs) only for those situations when fmincon + matlabFunction does not provide "right output". The rule that I use for determining not having "right output" is when either of the following is true: -

1. Exitflag non-positive

2. Exit flag positive but gradient and hessian re-evaluated using optimal x using double(subs) do not satisfy the optimality measures (gradient not zero or hessian not positive).

3. Constraints not satisfied if optimal x is substituted in to constraint equations and evaluated using double(subs)

Would appreciate any feedback.

Answer 1

[Christopher Berry]

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.

[Hari]

Thanks for the answer Christopher. I had guessed what you wrote but now it is formally confirmed.

Really like your statement "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." as this is what was happening

Answer 2

[Christopher Creutzig]

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.

[Hari]

Thanks for the tip Christopher. This helps. One question. I come to realize that double(symfun) is faster than double(subs) even though both have EXACT precision. Now for later u have recommended to modify double(subs) to double(subs(s,x,vpa(h))) so we get increased precision than matlabFunction (though not exact) and still faster than without vpa. Can I adopt the same approach in double(symfun) also. That is using double(symfun(vpa(arguments))) rather than double(symfun(arguments)) to get increased precision over matlabFunction (though not exact) and have better peformance over later.