Optimization of the function and Jacobian matrix

Structural knowledge and functional properties should be used by the symbolic environment in order to improve the efficiency of the resulting compiled code [21]. Optimisation can be used for more efficient evaluation of the function and its Jacobian matrix in Newton's method, in much the same way as in the Runge-Kutta example of section 6.4.4. For example, the optimization process can resolve matrix symmetries.

Illustration of these issues is given by a simple example. Consider a Jacobian matrix derived from a function of two variables f(x,y).

In[1]:= f[x_,y_] := {x + y + 2 x Sin[y]^2, x + y + x^2 Sin[2 y]};

In[2]:= jac[f_List,vars_List] := Outer[D,f,vars];

In[3]:= matrix = jac[f[x,y],{x,y}];

In[4]:= matrix //MatrixForm

Out[4]//MatrixForm=
            2
1 + 2 Sin[y]            1 + 4 x Cos[y] Sin[y]

                               2
1 + 2 x Sin[2 y]        1 + 2 x  Cos[2 y]

In fact this matrix turns out to be symmetric:

In[5]:= simpmat = Simplify[matrix];
In[6]:= simpmat //MatrixForm

Out[6]//MatrixForm=
2 - Cos[2 y]        1 + 2 x Sin[2 y]

                           2
1 + 2 x Sin[2 y]    1 + 2 x  Cos[2 y]

A considerable amount of computation can be saved if the property m₁₂=m₂₁ is utilised. Moreover, the functions Cos[2 y] and Sin[2 y] need only be evaluated once if the values are stored as temporary variables. This information is again resolved by syntactic optimization. An optimized computational sequence in FORTRAN may be obtained as:

In[7]:= FortranAssign[m, simpmat, AssignOptimize->True,
          OptimizePower->True,OptimizeTimes->False]
Out[7]//OutputForm=
        o1=cos(2.d0*y)
        o2=sin(2.d0*y)
        o3=1.d0+2.d0*o2*x
        m(1,1)=2.d0-o1
        m(1,2)=o3
        m(2,1)=o3
        m(2,2)=1.d0+2.d0*o1*x**2

Notice that the matrix components are optimized, but the encapsulating lists themselves are not, since the default setting OptimizeNull->{List} has been used to ignore such objects. Also we use the option OptimizeTimes->False to ignore sub-expressions with head Times. Such operations do not yield additional optimizations in this example and ignoring them speeds up the optimization process.