WIP Analytical Derivation of Sensitivity Equations

[FedericoV]

There are three major ways of calculating the time-dependent sensitivity of a system of differential equations:

• Automatic Differentiation
• Finite Differences
• Sensitivity Equations

Unlike FD or AD, sensitivity equations can only be used when each of the equations that maps the state variables to the differential equations can be differentiated with respect to the parameters. Here is a gentle introduction to the approach: http://www.cs.toronto.edu/~hzp/index_files/presentations/SONAD08hossein.pdf

In my work using ODE, I've found that Sensitivity Equations are particularly useful when trying to estimate the parameters of a system of ODE using gradient based methods like LM.

I've written some code for my own use in Python to calculate the expanded system of Sensitivity Equations for a given model using SymPy:

https://github.com/FedericoV/Sensitivity_Equations_Julia

The file that does all the work is this:

https://github.com/FedericoV/Sensitivity_Equations_Julia/blob/master/sympy_tools.py

I also started some very very preliminary work in Julia to extract the equations from a function using Julia's introspection abilities, but I didn't really get anywhere yet. Is this something more people are interested in?

[stevengj]

I assume you mean "finite differences", and not "finite elements"?

(If the equations are not differentiable, every method is going to run into trouble at some point.)

Note that one also should differentiate (no pun intended) between forward and reverse (aka adjoint methods) differentiation techniques. The latter is more efficient if you are differentiating with respect to a large number of parameters, but is more storage-intensive without checkpointing or similar complications.

[FedericoV]

Yes, oops. Just the doing (f(x+delta_x) - f(x)) / delta_x

On Wed, Mar 5, 2014 at 7:32 PM, Steven G. Johnson [email protected]:

I assume you mean "finite differences", and not "finite elements"?

Reply to this email directly or view it on GitHubhttps://github.com/JuliaLang/ODE.jl/issues/22#issuecomment-36776243 .

[ChrisRackauckas]

I'm interested in adding something like this in DifferentialEquations.jl. Have you made any more progress?

[FedericoV]

Hi Christopher,

I actually went ahead and developed a package to do this in pure python using sympy. It's very rough around the edges and definitely not fit for public consumption (yet) but I'm happy to put it on github if you want to take a look.

Currently, it can only handle simple differential equations - but it works well for my own use case.

Federico

On Sun, 7 Aug 2016 at 22:19 Christopher Rackauckas [email protected] wrote:

I'm interested in adding something like this in DifferentialEquations.jl. Have you made any more progress?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaLang/ODE.jl/issues/22#issuecomment-238105855, or mute the thread https://github.com/notifications/unsubscribe-auth/AAmddrH6Ty8QBSDnXDKtWeszXUatDXzpks5qdj3MgaJpZM4BnCoG .

[mauro3]

@FedericoV: this would definitely be cool to have in ODE.jl, if you ever find the time to look into it again.

@ChrisRackauckas: as you probably know, Sundials can do this. There is this example (probably non-functioning): https://github.com/JuliaLang/Sundials.jl/pull/50/files#diff-bcbf0efa3881f3441672006bef42ace7R1

[FedericoV]

Mauro: I might be mistaken but I think Sundials approximates the sensitivity equations/jacobian using a finite difference scheme? As far as I know, they don't have access to the symbolic form of the equations so I'm not sure how they could calculate the forward or backwards sensitivity equations without symbolic manipulation.

On Mon, 8 Aug 2016 at 09:50 Mauro [email protected] wrote:

@FedericoV https://github.com/FedericoV: this would definitely be cool to have in ODE.jl, if you ever find the time to look into it again.

@ChrisRackauckas https://github.com/ChrisRackauckas: as you probably know, Sundials can do this. There is this example (probably non-functioning): https://github.com/JuliaLang/Sundials.jl/pull/50/files#diff-bcbf0efa3881f3441672006bef42ace7R1

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/JuliaLang/ODE.jl/issues/22#issuecomment-238165066, or mute the thread https://github.com/notifications/unsubscribe-auth/AAmddiZynxtcJdh3OcEEUIsZCmnDYcofks5qdt_NgaJpZM4BnCoG .

[mauro3]

@FedericoV yes, you're right! I guess I was just referring to sensitivity analysis.

[pwl]

As far as I know, sensitivity analysis can be done via some sort of validated numerics, @lbenet or @PerezHz should know more about it. They were working on a high order Taylor method to solve ODEs in Julia, hopefully a package will be ready soon.

[ChrisRackauckas]

Validated numerics and interval can help. I am hoping the use of ArbReals can get that. But a standard sensitivity analysis could be implemented with something like ForwardDiff as well.

Supporting symbolic is hard when there really isn't a symbolic CAS for Julia. Right now there's SJulia, SymPy, and Nemo. SymPy is just a Python bridge and doesn't support the whole type system. SJulia needs to grow. Nemo isn't focused on these issues.

[lbenet]

I am not familiar with sensitivity analysis, but from what I read above, I guess it is related to some (maybe higher-order) stability-like analysis that considers small changes in the initial conditions or on the parameters of the equations. So I think the issue is if it is possible to obtain those derivatives from the equations of motion. Is this correct?

If so, TaylorSeries.jl may be useful.

julia> using TaylorSeries

julia> set_variables("x y μ") # This changes the internal parameters to deal with 3 variables
3-element Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 x + 𝒪(‖x‖⁷)
  1.0 y + 𝒪(‖x‖⁷)
  1.0 μ + 𝒪(‖x‖⁷)

julia> f1(x1,x2,mu) = x2 # RHS of 1st equation of motion
f1 (generic function with 1 method)

julia> f2(x1,x2,mu) = mu * (1-x1^2) * x2 - x1 # RHS of 2nd equation of motion
f2 (generic function with 1 method)

julia> # Create three `TaylorN` variables, named x, y, and μ.
       # Each one is expanded around zero.
       x = TaylorN(1)
       y = TaylorN(2)
       μ = TaylorN(3)
 1.0 μ + 𝒪(‖x‖⁷)

julia> # Now expand the equations of motion
       f1(x,y,μ)
 1.0 y + 𝒪(‖x‖⁷)

julia> f2(x,y,μ)
 - 1.0 x + 1.0 y μ - 1.0 x² y μ + 𝒪(‖x‖⁷)

julia> # You can expand around other values
       f2(x,y,3+μ)
 - 1.0 x + 3.0 y + 1.0 y μ - 3.0 x² y - 1.0 x² y μ + 𝒪(‖x‖⁷)

I hope this helps.

[ChrisRackauckas]

While an interesting project, in the docs I see it only supports:

exp, log, sqrt, sin, cos and tan

Thus it may be more fruitful to go directly from autodifferentiation (ForwardDiff).

[FedericoV]

If you use autodifferentiation directly, you don't really need to explicitly use sensitivity equations. The downside (and reason I didn't opt for that approach) is that it means that you can only use julia code, and are stuck with julia-lang only integrators.

Since there's a lot of incredibly high quality integrator libraries written in C, it's hard to use dual numbers there.

On Wed, 10 Aug 2016 at 06:03 Christopher Rackauckas < [email protected]> wrote:

While an interesting project, in the docs I see it only supports:

exp, log, sqrt, sin, cos and tan

Thus it may be more fruitful to go directly from autodifferentiation (ForwardDiff).

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaLang/ODE.jl/issues/22#issuecomment-238760431, or mute the thread https://github.com/notifications/unsubscribe-auth/AAmddngESM0KiYdn_1Pj29Cm1d3OW4fnks5qeU2egaJpZM4BnCoG .

[ChrisRackauckas]

It would add sensitivity analysis to the many Julia integrators we already have. Sundials already has its own. If ODE.jl did the same, then the only package this would really be missing is ODEInterface (which could be implemented using finite differences in a callback function). It would be easy in DifferentialEquations.jl to format all of these similarly in a sensitivity output. I think it's fine that not all integrators will support all functionality, it would just have to be noted in the docs.