0
$\begingroup$

I am trying to estimate a few unknown parameters of my continuous non-linear PDAE model (simulated through finite-volume method spatial discretisation, and time-stepping through method-of-lines). I am planning to do this by minimising the least-squares error of the model output to experimental data.

I read about the various optimisation solvers that can be used. However, I am having a weird situation. I do not have the exact analytical form of my function (since its a system of PDAES, which after discretisation yields a lot of DAEs). While, this is fairly common and the recommended method is the method of lines, I have a slightly different situation. I do have the Jacobian matrix (evaluated through Automatic Differentiation package 'CASADI').

I'd like to really take advantage of the power of the jacobian matrix, i.e. directional information towards solution. Furthermore, my function evaluation is very cheap, since my simulation model is hand-optimised and takes advantage of various intrinsic CPU instructions.

What would be the recommended software package for implementing a global optimisation for parameter estimation, given this situation ?

Improve this question
$\endgroup$

3 Answers 3

Reset to default
2
$\begingroup$

Since you have the Jacobian matrix, you can apply it within a Gauss-Newton or Levenberg-Marquardt method to effectively approximate the Hessian and gradient of your least squares objective function (the gradient is $J^{T}f$, and the Hessian is to first order $J^{T}J$.)

You could also use the Jacobian to compute the gradient of your least squares objective functions and then solve the problem with any gradient based scheme (e.g. BFGS or a limited memory BFGS method.) This might be a better approach if the number of parameters to be estimated is very large and you don't have sufficient space to deal with solving a system of equations involving $J^{T}J$. An alternative is to use an iterative method to solve the normal equations that arise in the GN or LM methods.

All of these methods can (with appropriate attention to detail) be guaranteed to converge to a local minimum of the least squares objective from any starting solution. This is sometimes confusingly called "global convergence", even though it doesn't mean convergence to a global minimum.

Assuming that your least squares objective is non-convex (most likely), then you'll want to consider combining a local nonlinear least squares solver with some kind of stochastic search in hopes of finding a global minimum. The simplest approach would be to use "multi-start", in which you run a local search from lots of randomly chosen solutions, and pick the best local minimum that you see.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Hi Brian, I think I didn't word the problem very clearly. I have the Jacobian only for the simulation model. But here I am trying to minimise the sum of squares of errors between the model's output and experimental data, in order to perform a parameter estimation of the model's parameters. I.e. the objective function is the sum of squares of difference of the experimental data and model output, whose Jacobian is not available. I am now thinking that even having the analytical Jacobian is a waste , and we might have to resort to derivative-free optimisation methods. Am I thinking on the right $\endgroup$
    – Dr Krishnakumar Gopalakrishnan
    Commented Oct 19, 2016 at 8:43
  • $\begingroup$ It's easy to get the gradient and approximate Hessian of the sum of squares from the Jacobian of the model predictions. Let $h(p)$ be the vector valued function of the model predictions with Jacobian $J(p)$. Let $f_{i}(p)=(h_{i}(p)-y_{i})$ be the difference between the model prediction and observation. Let $F(p)=\sum_{i=1}^{n} f_{i}(p)^{2}$. Then $\nabla F(p)=J(p)^{T}f(p)$ and $\nabla^{2} F(p) \approx J(p)^{T}J(p)$. $\endgroup$
    – Brian Borchers
    Commented Oct 19, 2016 at 12:13
0
$\begingroup$

You can either use some derivative free optimization, for which there are several packages available, or: If in addition to the jacobian you can obtain the hessian then there are many methods that you can write up. If you don't have constraints the newton method and steepest descent are very easy to write, and work rather well if you have a good method of estimating step size ( armijo conditions etc), constrained methods take more work to program. I recommend you the reference Numerical Optimization by J Nocedal, it contains pseudo code and references to software that would help you with your problem.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Yes, The CASADI package helps in obtaining the Hessian's as well. I thought Newton's method was a local solver ? Thanks for the reference book. $\endgroup$
    – Dr Krishnakumar Gopalakrishnan
    Commented Oct 18, 2016 at 0:35
  • $\begingroup$ If your function is convex a local optimum is a global optimum. There are global solvers for non convex functions ( e.g. By partitioning the space and using local solvers) but they are very computationally intensive and may not be able to guarantee optimality. A least square problem is convex though, unless you have some weird constraints. $\endgroup$
    – Septimus G
    Commented Oct 18, 2016 at 0:41
  • $\begingroup$ Ah. Thank you for the explanation. In that case, maybe I should be able to try the 'CVX' package? $\endgroup$
    – Dr Krishnakumar Gopalakrishnan
    Commented Oct 18, 2016 at 0:48
0
$\begingroup$

Let me clarify from my previous comment that I had linear least squares problems in mind, non linear least squares need to meet several conditions in order to be convex. I would look for some simple code for non linear least square problems (which are basically non linear methods that take advantage of the structure of the problem to reduce in computations, e.g. Gauss-newton and levenberg-marquardt) and modify it to add the jacobian and hessian from your program. You may need some understanding of the methods to do this. Most packages will expect that you provide the objective function.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.