Improving Accuracy and Performance in Optimization

Open Live Script

This example demonstrates that the Symbolic Math Toolbox™ helps minimize errors when solving a nonlinear system of equations.

This example uses the following Symbolic Math Toolbox capabilities:

Computing the analytical gradient using gradient
Computing the analytical jacobian using jacobian
Converting symbolic results into numeric functions using matlabFunction
Analytical visualization with fplot

The goal is to find the minimum of the Rosenbrock function, or the so-called "banana" function.

$f (x) = \sum_{i = 1}^{n / 2} 100 (x_{2 i} - x_{2 i - 1}^{2})^{2} + (1 - x_{2 i - 1})^{2}$

Use fplot to create a quick visualization of a Rosenbrock function of two variables.

syms x y
a=1; b=100;
f(x,y)=(a-x)^2+b*(y-x^2)^2

f(x, y) =  ${(x - 1)}^{2} + 100 {(y - x^{2})}^{2}$

fsurf(f,[-2 2 -1 3])
view([-45 50])
colormap jet
colorbar

Figure contains an axes object. The axes object contains an object of type functionsurface.

Motivation

The Solve Nonlinear System Without and Including Jacobian (Optimization Toolbox) example in Optimization Toolbox™ solves the same problem by using the fsolve (Optimization Toolbox) function. The workflow shown in that example has two potential sources of errors. You would need to

Translate the equations for the Rosenbrock function and Jacobian from the math form in the text to numeric code.
Explicitly calculate the Jacobian. For complicated equations, this task is prone to errors.

This example shows that using symbolic math to describe the problem statement and subsequent steps minimizes or even eliminates such errors.

Rewrite Rosenbrock Function as a System of Nonlinear Equations

First, convert the Rosenbrock function f to a system of nonlinear equations F, where F is a gradient of 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);

Show the first 10 equations:

F(1:10)

ans = 
 $(\begin{array}{c} 2 x_{1} - 400 x_{1} (x_{2} - {x_{1}}^{2}) - 2 \\ 200 x_{2} - 200 {x_{1}}^{2} \\ 2 x_{3} - 400 x_{3} (x_{4} - {x_{3}}^{2}) - 2 \\ 200 x_{4} - 200 {x_{3}}^{2} \\ 2 x_{5} - 400 x_{5} (x_{6} - {x_{5}}^{2}) - 2 \\ 200 x_{6} - 200 {x_{5}}^{2} \\ 2 x_{7} - 400 x_{7} (x_{8} - {x_{7}}^{2}) - 2 \\ 200 x_{8} - 200 {x_{7}}^{2} \\ 2 x_{9} - 400 x_{9} (x_{10} - {x_{9}}^{2}) - 2 \\ 200 x_{10} - 200 {x_{9}}^{2} \end{array})$

To match the intermediate results shown in the Solve Nonlinear System Without and Including Jacobian (Optimization Toolbox), use the same form of 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;

Now the systems of equations are identical:

F(1:10)

ans = 
 $(\begin{array}{c} 1 - x_{1} \\ 10 x_{2} - 10 {x_{1}}^{2} \\ 1 - x_{3} \\ 10 x_{4} - 10 {x_{3}}^{2} \\ 1 - x_{5} \\ 10 x_{6} - 10 {x_{5}}^{2} \\ 1 - x_{7} \\ 10 x_{8} - 10 {x_{7}}^{2} \\ 1 - x_{9} \\ 10 x_{10} - 10 {x_{9}}^{2} \end{array})$

Calculate Jacobian of Matrix Representing the System of Equations

Use jacobian to calculate the Jacobian of F. This function calculates the Jacobian symbolically, thus avoiding errors associated with numerical approximations of derivatives.

JF = jacobian(F,x);

Show the first 10 rows and columns of the Jacobian matrix.

JF(1:10,1:10)

ans = 
 $(\begin{array}{c} - 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ - 20 x_{1} & 10 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - 20 x_{3} & 10 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - 20 x_{5} & 10 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 20 x_{7} & 10 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 20 x_{9} & 10 \end{array})$

Most of the elements of the Jacobian matrix JF are zeros. Nevertheless, when you convert this matrix to a matlabFunction, the result is a dense numeric matrix. Often, operations on sparse matrices are more efficient than the same operations on dense matrices.

Therefore, to create efficient code, convert only the nonzero elements of JF to a symbolic vector before invoking matlabFunction. is and js represent the sparsity pattern of JF.

[is,js] = find(JF);
JF = JF(JF~=0);

Convert System of Equations and Jacobian to a MATLAB Function

The system of equations F, representing the Rosenbrock function, is a symbolic matrix that consists of symbolic expressions. To be able to solve it with the fsolve function, convert this system to a matlabFunction. At this step, it is convenient to convert both F and its Jacobian, JF, to a single file-based MATLAB function, FJFfun. In principle, this allows F and JF to reuse variables. Alternatively, you can convert them to individual MATLAB functions.

matlabFunction([F;JF],'var',x,'file','FJFfun');

FJFfun accepts variables as a list, but fsolve expects a vector of variables. The function genericObj converts the vector to a list. Anonymous function bananaObj is defined to fix argument values for genericObj. Note the use of the sparsity pattern in the function genericObj. Further note that this approach of using FJFfun and genericObj together is not tied to the particular expression represented by F and can be used as is for any F.

bananaobj = @(y) genericObj(y,@FJFfun,is,js)

bananaobj = function_handle with value:
    @(y)genericObj(y,@FJFfun,is,js)

Use `fsolve` to Numerically Solve the System of Nonlinear Equations

Use fsolve for the system of equations, converted to MATLAB function. Use the starting point x(i) = –1.9 for the odd indices, and x(i) = 2 for the even indices. Set 'Display' to'iter' to see the solver's progress. Set 'Jacobian' to 'on' to use the Jacobian defined in bananaobj.

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.

Use vpasolve to Numerically Solve the System of Nonlinear Equations

The system of nonlinear equations, F, can be alternatively solved with the vpasolve function, which is a numeric solver available in Symbolic Math Toolbox. A key difference between vpasolve and fsolve is that fsolve uses double-precision arithmetic to solve the equations, whereas vpasolve uses variable-precision arithmetic and therefore, lets you control the accuracy of computations.

sol2 = vpasolve(F);

If you assign the solution of a system of equations to one variable, sol2, then vpasolve returns the solutions in form of a structure array. You can access each solution (each field of the structure array):

sol2.x1

ans =  $1.0$

You also can convert sol2 to a symbolic vector, and then access a range of the solutions. For example, display the 10th to 20th solutions:

for k=(1:64)
  sol2_array(k) = sol2.(char(x(k)));
end
sol2_array(10:20)

ans =  $(\begin{array}{c} 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 \end{array})$

Helper Functions

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