Looking for ways to speed up the numeric evaluation of a symbolic expression in Matlab

{Summary: I have a symbolic expression DCritnF expressed in terms of two variables x1 and x2. I need to find it's numeric value and I used combination of double and subs as given below.

FuncVal = double(subs(DCritnF,[x1,x2],[x(1),x(2)]));


Challenge is that since DCritnF is a very convoluted expression (generated by symbolic toolbox using combination of matrix algebra and some calculus), evaluating FuncVal takes atleast 0.1 seconds (used timeit function). This is unacceptable as I'm performing optimization of a function and there are multiple function calls to above line of code and a single optimization run takes anywhere from 4 seconds to 40 seconds depending on type of algorithm used. I already determined using profiler that the one and only culprit is the above line of code. How to speed up the above function evaluation?}

Details: Im working on global optimization of nonlinear and nonconvex functions by supplying different start points for the decision variables through a For Loop (only have the optimization toolbox in matlab and dont have the global optimization toolbox so simulating multistart functionality).

The objective function I have is quite complex, based on a bunch of matrix computations (transpose, multiplication, inverse and finally determinant) using the gradient of function in equation 11 of http://www.amstat.org/publications/jse/v12n3/goos.html . The objective function (DCritn) is listed in the end of this question for reference.

Here q1 and q2 are parameters and x1/x2 are the decision variables.

For any given value of q1 and q2, like let's say 0.6 and 0.3, I try to find the values of x1/x2 for which DCritn is minimum. In order to do the optimization, I supply different start points for x and use different set of algorithms like interior point, SQP etc to see which one generates best possible local mimima.

syms x1 x2 q1 q2;
q=[q1,q2];
x=[x1,x2];
qvalue=[0.6,0.3];


%bunch of symbolic computations to arrive at symbolic form of DCritn which is the objective function (see end result after signature)

DCRGrad = jacobian(DCritn,x).';%Gradient of the the objective function
DCritnF = subs(DCritn,q,qvalue.');%DCritn objective function specific to qvalue

xstart=[5,10];% supplying start points for the optimization

%fmincon code to optimize DCritn using different start values of x. Options for using gradient specified (fmincon code skipped here)


I have a seperate function file "Objfungrad", provided below, which fmincon can tap in to for calculating function value and gradient for any given value of x

function [FuncVal,gradVal] = Objfungrad(x)

global x1 x2;
syms x1 x2;

try
% Objective function
FuncVal = double(subs(DCritnF,[x1,x2],[x(1),x(2)]));
catch exception
FuncVal = NaN;
end
% Gradient of the objective function
if nargout  > 1
end


As per profiling tool, more than 95% of optimization time is spent in evaluating lines of code pertaining to FuncVal and gradVal.

If instead of double and subs functions used above, If I had hard-coded the complete symbolic expression for FuncVal and gradVal using qvalue as fixed parameters, then optimization gets completed in less than 10% of the time. Unfortunately this hard-coding is not possible, as parameters q can change.

So question is, how can I speed up evaluation of a symbolic expression based on double and subs?

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)))

• If you hard-code your expression, try factorisation of your expression. I see many common terms in the sums. – AlexE Jun 23 '14 at 7:54
• Welcome to SciComp! If you indent your code by four spaces (or use the {} button), it will become much more readable. – Christian Clason Jun 23 '14 at 13:32
• Thanks Christian, I see that you have already edited it. Looks very good. Will remember to use {} button. – Hari Jun 23 '14 at 14:06
• AlexE, I tried factoring DCritn using factor function and then doing subs, double etc but this only increases the time taken to compute. – Hari Jun 23 '14 at 14:30
• I did a little more research on this and realized that the bulk of the time is actually taken by subs function. Double takes less than 20% of the time and rest is on subs. Infact if I replace the double by vpa it becomes faster. Still, how do I overcome the slowness with subs since it represents the majority of the evaluation time – Hari Jun 23 '14 at 15:39

Running your expression through Mathematica's FullSimplify produced this:

(exp(2*(q1 + q2)*(x1 + x2))*(q1 - q2)^6)/(q1^2*(exp(q2*(x1 + x2))*q1*(x1 - x2) +
exp(q1*(x1 + x2))*q2*(x1 - x2) +
exp(q1*x1 + q2*x2)*(-q2*x1 + q1*x2 + q1*(-q1 + q2)*x1*x2) +
exp(q2*x1 + q1*x2)*(q2*x2 + q1*x1 (-1 + q1*x2 - q2*x2)))^2)


This is much shorter at least.

• Thanks for the tip on using simplify. I used the MATLAB equivalent and it came reasonably close to what you provided above. But the problem is still not resolved. The reason is that at some point I need to apply the subs (and double) function and this is where the major overhead is present (infact it is in subs). So it does not matter much in terms of time if the expression is big or small, and the subsequent subs function takes up majority of time..Any thoughts on how I could overcome it? – Hari Jun 24 '14 at 5:36
• I have one thought, but not sure if possible. Since I see the issue being with sub function, I was thinking if I could avoid it somehow or use proxies. Going back to the expression you provided above, is it possible to replace x1 with x(1) and x2 with x(2) in above using code. If yes (how?), then I could feed this modified expression in to the function evaluation file that fmincon calls and in this new approach using double alone is enough – Hari Jun 24 '14 at 6:20

Have you tried matlabFunction? For example,

syms x y real
f = matlabFunction(x^2 + y^2, 'vars', {[x, y]})


will create a (non-symbolic) function f that takes a two-component vector as input (rather than two different input arguments). The function is vectorized automatically, which potentially allows MATLAB to optimize its evaluation.

• matlabFunction was surprisingly worse for me as compared to using combination of subs/double so thats why I used later. Ofcourse I'm still not proficient in matlab and came to know about vectorization only recently (later not implemented in my code). Since Im running fmincon from many different start points (points generated by me), can this be run in a vectorized manner? – Hari Jun 24 '14 at 5:56
• I'm very surprised about it not working as fast as subs/double, I'd encourage you to try it again taking advantage of vectorization (which your Objfungrad function doesn't currently use). I expect vectorization to be useful in your case inside for fmincon for a particular starting point. The other scenario (many start points) seems more suited to parallel evaluation. – Juan M. Bello-Rivas Jun 24 '14 at 12:58
• @Hari, make sure that you're not calling matlabFunction on every iteration. it only needs to be run once. The the output function can be used. You should also combine this answer with that from @VictorLiu to make the resultant function more efficient. Or you can just take the equation from @VictorLiu's answer and make your own vectorized function. There's nothing magic about matlabFunction. Your objective function should be doing everything in double precisions and nothing in symbolic math when you're done. – horchler Jun 24 '14 at 15:58
• Thanks @horchler , @VictorLou , @Juanmbellorivas . By using following in main code DCritnF = matlabFunction(simplify(subs(DCritn,q,qvalue.'),'IgnoreAnalyticConstraints', true,'steps',500),'vars',t1); it is 10 times faster and makes me super happy. Kudos!. Two related ques: Q1. Within "Objfungrad" function unless I put global DCritnF; the "DCritnF" created using matlabFunction is not recognized. I read global should be avoided? Q2. @Juanmbellorivas mentions that with matlabFunction code automatically gets vectorized. I understood vectoring to mean replacing for loops with matrix opns? – Hari Jun 24 '14 at 19:01
• One admin question. Should I clicking on the yes/check mark against both VictorLou and Juanmbellorivas answers ?. What is the protocol – Hari Jun 24 '14 at 19:04

If the expression is as "simple" (I realize that's a relative term) as VictorLiu shows, then you might try hand-coding it.

If hand-coding an equivalent MATLAB expression is unwieldy (for instance, if you have other, more complicated expressions to work with), you could try using MATLAB's code generation facilities. MATLAB can convert MuPAD symbolic expressions into MATLAB, Fortran, or C code. If speed is truly critical, converting the symbolic expression into C code and then writing a MEX wrapper around the C code to call it from MATLAB might be a viable option.

• Thanks for the comment Geoff. You are indeed correct that I would have other expressions to work with. For the present exercise, starting function I used was from equation 11 of amstat.org/publications/jse/v12n3/goos.html and it went through a couple of matrix computations to arrive at the present form of DCritnF. But in reality, my starting function could be any equation from physical/chemical/biological world. Additionally I took only 2 unknown but in real problems number of unknowns would be higher. I guess if am not able to solve speed issues, I will have to look in to C code – Hari Jun 24 '14 at 6:04