Chaostools tangent_integrator stiff solvers

Are you creating any cache variables? If you are, are they dual cached?

I tested some more algorithms, and so far none of Kvaerno5, `KenCarp4, Rodas4P worked.

cds = ContinuousDynamicalSystem(lorenz_iip,x_start,param)
  • leaving out the Jacobian (I assume its auto diffed then), does not work either.

I am just running the code above in a Jupyter notebook. What do you mean by cached/dual cached in this context?

To the best of my understanding, tangent_integrator propagates everything to __init and ODEProblem. See here: https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/continuous.jl#L50-L68 :

function tangent_integrator(ds::CDS{IIP}, Q0::AbstractMatrix;
    u0 = ds.u0, diffeq...) where {IIP}

    t0 = ds.t0
    Q = safe_matrix_type(Val{IIP}(), Q0)
    u = safe_state_type(Val{IIP}(), u0)

    k = size(Q)[2]
    k > dimension(ds) && throw(ArgumentError(
    "It is not possible to evolve more tangent vectors than the system's dimension!"
    ))

    tangentf = create_tangent(ds.f, ds.jacobian, ds.J, Val{IIP}(), Val{k}())
    tanprob = ODEProblem{IIP}(tangentf, hcat(u, Q), (t0, typeof(t0)(Inf)), ds.p)

    solver = _get_solver(diffeq)
    return __init(tanprob, solver; DEFAULT_DIFFEQ_KWARGS..., internalnorm = _tannorm,
                  save_everystep = false, diffeq...)
end

(also, since the problem comes from using different integration algorithms, it cannot stem from DynamicalSystems.jl since the step! functionality is beyond our control; we only initialize the states and functions and give everything to init)

Can you point to the implementation of create_tangent? How do you tie together the “f” and “jacobian”? Could this be something, the stiff solvers do not like?

The tangent space function is created generally in this part of the code: https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/dynamicalsystem.jl#L342

e.g. for the case of in-place system with provided jacobian function one has (https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/dynamicalsystem.jl#L343-L355 )

function create_tangent(@nospecialize(f::F), @nospecialize(jacobian::JAC), J::JM,
    ::Val{true}, ::Val{k}) where {F, JAC, JM, k}
    J = deepcopy(J)
    tangentf = (du, u, p, t) -> begin
        uv = @view u[:, 1]
        f(view(du, :, 1), uv, p, t)
        jacobian(J, uv, p, t)
        mul!((@view du[:, 2:(k+1)]), J, (@view u[:, 2:(k+1)]))
        nothing
    end
    return tangentf
end

I can not post the whole stack trace here, as it is too long (over 32000 characters). I can post it somewhere else, if someone is interested :slight_smile:

@Datseris So if i switch to Rodas5(autodiff=false) this works. So it seems that differentiating tangentf fails for some reason.

1 Like

Hm, nice that you found the root of the problem. I’ve cited above how the tangent function is created, but I don’t really have a clue of what Rodas does (or what it is in general), so I can’t really help on why this fails.

Maybe Chris has some ideas?

Hmm, I do know way too little to track down this error further. Are the non stiff algorithms not using ForwardDiff ? Because they work flawlessly here. Maybe @ChrisRackauckas can help?

Yup, you’re creating cache variables which aren’t dual cached which violates the assumptions of autodiff.

http://docs.juliadiffeq.org/latest/basics/faq.html#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods…?-1

1 Like

Ah cool, this seems like something that can be fixed!

@jamblejoe if you want to open a PR for this on DynamicalSystemsBase you are welcomed! (I won’t have time to fix it myself)

Yeah, it just takes some manpower. And it really is fundamental to autodiff, because dual numbers take up (N+1)*bitsize memory. So for example if you have a dual number of Float64s with 3 partials, that dual number is actually a 4 dimensional number of 64 bits each, which is why there’s no way this is ever going to work without caching a bit differently. But the trick for how to use dispatch + reinterpret is documented in the FAQ there, which essentially builds a cache for the two possible sizes, chooses between them, and fixes the tag so ForwardDiff doesn’t error.

1 Like

@ChrisRackauckas now I got what you mean, Chris, thanks! I assume, as the ForwardDiff package only errors, when using solvers for stiff problems, that the non-stiff solvers are not using ForwardDiff? Are they still auto-differentiating by default or do they approximate by finite differences?

So yes, until DynamicalSystems.jl makes its caches compatible with ForwardDiff, it cannot use automatic differentiation for the stiff ODE solvers. That is fine since they all have a switch to use numerical differencing (Rodas5(autodiff=false)). The non-stiff ODE solvers aren’t building Jacobian so they don’t use any differencing.

@ChrisRackauckas I played around with stiff solvers and autodiff=false a bit and, compared to the non-stiff solvers in parameter regimes of my systems, where they work, they seem to be very slow. Providing a Jacobian by hand to avoid finite differences is not an option as they are simply too big.

I had a (short) look at ParametrizedFunctions and ModellingToolkit and I really like the idea of symbolically calculating Jacobians, etc. It would be really nice, and less error prone for the user, if one can simply define the actual system, say f, via e.g. something like @ode_def, and not only have the Jacobian Df, but also the Jacobian of the coupled system hcat(f,Df) automatically calculated, which then can be used to solve the coupled system without the need of auto-differentiating or finite differences :slight_smile:

Is there a way to achieve this atm?

ModelingToolkit.jl is a symbolic system, so there is a way to get it once you’d symbolicicized to it. Explain what you’re looking to generate and maybe I can help you get there. Will you be at the JuliaCon workshop? We plan to go over this.

Also look into SparseDiffTools.jl

Automatic detection of sparsity + Jacobian coloring will get you pretty close to analytical on a large sparse system.

The whole idea of using the tangent_integrator is to calculate Lyapunov exponents. Say we have a differential equation

d/dt x = f(x)

with Jacobian Df, then we need to solve the coupled system

d/dt x = f(x)
d/dt W = Df(x)*W

This is what create_tangent does - it ties both differential equations together in one differential equation, say g, which then can be solved. The solver then might need to compute Dg. the Jacobian of g, for solving the for some initial values.

So if we have symbolically defined f, we can calculate the Jacobian Df of f, couple the system and get g and the symbolically calculate Dg. So we only need the define f and the rest could be done by the symbolic engine. This would save a lot of code writing as the dimension of Dg grows asymptotically to the power of 4 of the dimension of f, e.g. for the Lorenz system
dim(f) = 3
dim(Df) = 3x3 = 9
dim(g) = 3 + 9 = 12
dim(Dg) = 12x12 =144.

SparseDiffTools would be very handy here, as the dimensions grow quickly.

I will not be at JuliaCon this year, unfortunately. I picked up Julia two months ago and am just scratching the surface right now :slight_smile:

JuliaCon is going on, so remind me to do this, but a good Hackathon thing would be to implement the Jacobian function. Then you could literally write ModelingToolkit.jacobian(f)*W and it would generate that equation (matrix multiplication already works). Right now, I’d just calculate the Jacobian using the primitive expend_derivative(ex), looping through and writing the derivative operator (which is how the Jacobian function would be created).

I will have a look at this! Thanks Chris for all your help!