2
$\begingroup$

I have a RKF45 numerical integrator that simulates polymerization of proteins using CUDA. It does so by tracking the populations of discrete length polymers, e.g. monomers, dimers, trimers, etc. all the way up to 8192-mers using the equations of the type that follow.

$$\frac{dc_r}{dt}= k_nc_r^{n_c}+2k_a(c_{r-1}-c_r) + k_m(2\sum^{max}_{s=r+2}c_s-(r-1)c_r) $$

where r represents the length of the polymer and c its concentration at that given time. The k's are the parameters that I'd like to minimize an objective function with respect to. These equations are integrated over about 20,000-30,000 timesteps.

My goal is to start fitting the model to experimental data. The first type of fit I'm trying is fitting a value derived from these polymer populations at arbitrary timesteps. For instance, I want to calculate the average length of polymers at five different times in the simulation/integration and compare them to externally supplied data. The average length is derived for each timestep from the concentrations like this :

$$L(t) = \frac{M(t)}{N(t)}$$ $$M(t) = \sum^{max}_{s=n_c}s*c_s(t)$$ $$N(t) = \sum^{max}_{s=n_c}c_s(t)$$

The objective function I'm trying to minimize thusly looks like this :

$$y(t)=\sqrt{\sum_{fit data}(L_{calculated}-L_{experimental})^2}$$

I've already implemented a solution to do this and run it through the following norm and minimize said norm/objective function using Scipy's Nelder-Mead optimize method. It seems to work, but if it's possible to make it perform better, I'd like to make it happen. I want to try and utilize more robust and swifter converging methods like L-BFGS-B, especially because they also allow me to impose bounds on the parameters, that need to remain positive.

Now, the Runge-Kutta algorithm, for the unfamiliar, is here : http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

Due to the recursive nature of the Runge Kutta calculations such as $$k_1=f(y_n)$$ $$k_2=f(y_n+a_{21}k_1)$$ $$k_3=f(y_n+a_{31}k_1+a_{32}k_2)$$

I'm having issues figuring out how to calculate the Jacobian, which I believe would be this.

$$J(t) = [\frac{dy}{dk_a} \frac{dy}{dk_m} \frac{dy}{dk_n}]$$

So, if I framed this question properly, is it possible at some point during the RKF calculation to construct a Jacobian to pass out to Scipy along with the objective function evaluation? I'd like to get an optimized optimization script up and running ASAP.

Improve this question
$\endgroup$
1
  • $\begingroup$ An alternative (if this is too hard to do directly) would be to solve another set of ODEs satisfied by $\frac{\partial c_r}{\partial k_a}$ (and the other partial derivatives), whose r.h.s. depends on the approximate solution $c_r$. $\endgroup$
    – Kirill
    Sep 16, 2014 at 5:23

1 Answer 1

Reset to default
3
$\begingroup$

So, if I framed this question properly, is it possible at some point during the RKF calculation to construct a Jacobian to pass out to Scipy along with the objective function evaluation?

The Runge-Kutta-Fehlberg 4(5) method you describe is an explicit numerical method for solving an ODE. As such, it does not require the solution of a nonlinear algebraic equation; Jacobian matrices are not required for that algorithm.

Is it possible to construct a Jacobian? Sure; the RKF4(5) method has nothing to do with it.

What you want are parameter sensitivities, namely $\partial{c_{r}}/\partial{k_{a}}$, $\partial{c_{r}}/\partial{k_{m}}$, etc., so that you can calculate the derivatives you need.

These can be calculated numerically in a few ways, to my knowledge; all of them involve solving a set of auxiliary differential equations:

  1. the adjoint method; see Steven G. Johnson's notes, and this paper by Petzold, et al., among others by Petzold. Adjoint methods work best when you have fewer state variables than parameters. (That is, the dimension of your $c_{r}$ is smaller than the number of parameters you have.)
  2. the simultaneous direct method, where you solve the sensitivity ODE system and your original state ODE system at the same time as a larger coupled ODE system
  3. the staggered direct method (see M. Caracotsios, W.E. Stewart, Sensitivity analysis of initial value problems with mixed ODEs and algebraic equations, Comput. Chem. Engrg. 9 (4) (1985) 359–365)
  4. the simultaneous corrector method, where you solve for your state variables and sensitivities at the same time (see T. Maly, L.R. Petzold, Numerical methods and software for sensitivity analysis of differential–algebraic systems, Appl. Numer. Math. 20 (1997) 57–79)
  5. the staggered corrector method (by Barton and Feehery; disclaimer: Barton was one of my advisors); see Feehery's thesis, Chapter 4, also see W.F. Feehery, J.E. Tolsma, P.I. Barton, Efficient sensitivity analysis of large-scale differential–algebraic systems, Appl. Numer. Math. 25 (1997) 41–54

In chemistry and chemical engineering, there are a few other types of sensitivity analysis available, but I don't think of them will give you the Jacobian matrix information you need for parameter fitting.

One final warning: ODE-constrained optimization problems are notoriously tricky, and generally nonconvex. It is possible to apply a standard nonlinear programming solver to these problems and obtain local minima that are globally suboptimal. See the thesis of Joseph K. Scott for examples. (Disclaimer: Dr. Scott was a former colleague.)

Improve this answer
$\endgroup$
4
  • $\begingroup$ (+1) "all of them involve solving a set of auxiliary differential equations" Is this true? The ode solver is, in the end, a complicated rational function of the initial conditions and other inputs. So I think it makes sense to ask how one can easily differentiate that specific function with respect to its inputs. $\endgroup$
    – Kirill
    Sep 16, 2014 at 22:56
  • $\begingroup$ @Kirill: You could discretize in time and then take derivatives with respect to parameters. That approach would be collocation and Lorenz Biegler has a bunch of papers that use that approach in optimization. Normally, when you do that, though, you give up adaptive time-stepping (because that would cause nonsmoothness in what are effectively the sensitivity derivatives), so you probably wouldn't use RKF4(5). If you take derivatives of the continuous equation, adaptive time-stepping can still be used sometimes on the sensitivity ODE system. $\endgroup$
    – Geoff Oxberry
    Sep 16, 2014 at 23:23
  • $\begingroup$ Geoff, thanks for the great answer, as usual. As I'm parsing through the proffered papers, am I reading in the second link you posted that the forward method, not the account method, is the one I'm looking for? I think you might have gotten them confused, but I know nothing about sensitivity analysis. $\endgroup$
    – Hair of Slytherin
    Sep 17, 2014 at 15:26
  • $\begingroup$ @KarstenChu: Yes, you're right. I was thinking about PDE-constrained optimization and mixed up the applications in my head; adjoint is better when you have more parameters than state variables. $\endgroup$
    – Geoff Oxberry
    Sep 17, 2014 at 17:57

Your Answer

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

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