What is considered good practice for coding up a function's derivatives, e.g. its Jacobian and Hessian matrices?

Hi,

I'm currently practicing numerical root-finding, using simple sets of nonlinear equations and writing my own solvers -- with an eye towards using Matlab's fsolve for my real problem that's in many more variables. This way, I'll be using fsolve while having a good sense of what root-finding methods generally do, what the pros and cons of various methods are, etc.

Since my practice problems are easy for now, I can differentiate the functions to compute the Jacobian matrices myself, and then code it up in the script files. However, I imagine it's not best practice to continue doing this, especially when I start considering many more variables, and the Jacobian / Hessian matrices get larger.

What's considered "best practice" for computing derivatives to use for, say, a root-finding method? Should I purchase and use the Symbolic toolbox, or is there another way to approach differentiation, without having to do it by hand?

I think I'll eventually get into some "finite-differencing" methods, but I'm not there yet and know nothing about them -- I'm maybe a few weeks away. So, any thoughts and recommendations are welcome.

Thanks,

Answer 1

[Matt J]

However, I imagine it's not best practice to continue doing this, especially when I start considering many more variables, and the Jacobian / Hessian matrices get larger.

On the contrary, if you have many variables, it is best to differentiate the functions yourself analytically. However, the efficiency of doing it this way only manifests if you express your Jacobians in vectorized terms. For example, you should recognize that a quadratic form x.'*Q*x has a Jacobian efficiently implemented as (2*x.')*Q. The Symbolic Toolbox is not capable of differentiating in matrix-vector form and will give you very unwieldy expressions expanded into a large number of scalar variables.

There are tools on the File Exchange that will help with numerical differentiation if you want to go that route,

https://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation

but the most efficient code will always come from an analytical differentiation customized to the objective function at hand.

[Noob]

Hi Matt,

Efficient because if I differentiate by hand, and supply the Jacobians / Hessians in the script file, then Matlab doesn't have to spend the time differentiating in several variables for me. Is that what you're saying? If so, that's interesting -- and I'll stick with analytical differentiation for now, and look for efficient ways of doing so, like in your quadratic form example.

Please let me know if I misunderstood you, and thanks for your answer,

[Matt J]

No, there's never any reason you would have Matlab do the calculus for you in real time. If you were to use the Symbolic Toolbox for this, you would just use matlabFunction() offline to do the calculus and write you a function for the Jacobian evaluation before you would ever even invoke FSOLVE. The problem is that what it will give you is not very efficient code. Here's a 10-variable example, for the quadratic form case that I outlined above,

>> Q=randi(4,10);Q=Q*Q.';
>> x=sym('x',[10,1])
 
x =
 
  x1
  x2
  x3
  x4
  x5
  x6
  x7
  x8
  x9
 x10

As we've said, the efficiently coded expression for the Jacobian would be J=@(x) 2*x.'*Q. But using the Symbolic Toolbox's commands matlabFunction() and jacobian() gives us the messier, non-vectorized expression below.

 >> J=matlabFunction(jacobian(x.'*Q*x))

J =

  function_handle with value:

    @(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10)[x1.*8.8e+1+x2.*1.04e+2+x3.*1.34e+2+x4.*9.0e+1+x5.*1.08e+2+x6.*9.4e+1+x7.*1.14e+2+x8.*1.12e+2+x9.*9.6e+1+x10.*8.4e+1,x1.*8.4e+1+x2.*1.14e+2+x3.*1.46e+2+x4.*1.16e+2+x5.*1.32e+2+x6.*1.16e+2+x7.*1.36e+2+x8.*1.1e+2+x9.*9.6e+1+x10.*1.3e+2,x1.*1.04e+2+x2.*1.38e+2+x3.*1.7e+2+x4.*1.2e+2+x5.*1.36e+2+x6.*1.22e+2+x7.*1.46e+2+x8.*1.36e+2+x9.*1.2e+2+x10.*1.14e+2,x1.*1.34e+2+x2.*1.7e+2+x3.*2.24e+2+x4.*1.54e+2+x5.*1.78e+2+x6.*1.42e+2+x7.*1.82e+2+x8.*1.7e+2+x9.*1.52e+2+x10.*1.46e+2,x1.*9.0e+1+x2.*1.2e+2+x3.*1.54e+2+x4.*1.34e+2+x5.*1.26e+2+x6.*1.18e+2+x7.*1.26e+2+x8.*1.24e+2+x9.*1.0e+2+x10.*1.16e+2,x1.*1.08e+2+x2.*1.36e+2+x3.*1.78e+2+x4.*1.26e+2+x5.*1.62e+2+x6.*1.28e+2+x7.*1.5e+2+x8.*1.4e+2+x9.*1.24e+2+x10.*1.32e+2,x1.*9.4e+1+x2.*1.22e+2+x3.*1.42e+2+x4.*1.18e+2+x5.*1.28e+2+x6.*1.44e+2+x7.*1.36e+2+x8.*1.38e+2+x9.*1.08e+2+x10.*1.16e+2,x1.*1.14e+2+x2.*1.46e+2+x3.*1.82e+2+x4.*1.26e+2+x5.*1.5e+2+x6.*1.36e+2+x7.*1.86e+2+x8.*1.58e+2+x9.*1.38e+2+x10.*1.36e+2,x1.*1.12e+2+x2.*1.36e+2+x3.*1.7e+2+x4.*1.24e+2+x5.*1.4e+2+x6.*1.38e+2+x7.*1.58e+2+x8.*1.7e+2+x9.*1.32e+2+x10.*1.1e+2,x1.*9.6e+1+x2.*1.2e+2+x3.*1.52e+2+x4.*1.0e+2+x5.*1.24e+2+x6.*1.08e+2+x7.*1.38e+2+x8.*1.32e+2+x9.*1.32e+2+x10.*9.6e+1]

The code is correct, but it has written everything out as scalar expressions - it doesn't implement the calculation using Matlab's fast matrix multiplication capabilities. Imagine now what you would get with 100 unknown variables instead of just 10, and how slow it would be to execute..

[Noob]

Hi Matt,

Ok, got it -- thanks again!

Take care,

[Matt J]

You are quite welcome, but please Accept-click the answer since we seem to have resolved your question.

[Noob]

Hi Matt,

Quick question: what about having Matlab invert my matrices? I'm currently using inv( ) to invert my Jacobian matrices and then left-multiplying with a function (multivariable Netwon's method).

Should I invert the matrices by hand, and then code it in the script file without using the inv( ) function?

Thanks,

[Matt J]

You should never form the inverse of a matrix to solve an equation J*x=-g. You should be solving it with mldivide():

x=J\(-g)