Disable AD for Continuation

[fh9g12]

Is it possible to disable automatic differentiation and use finite difference for the calculation of period orbits and normal forms?

I am just starting with Julia and BifurcationKit and as a first step I was trying to replicate results I have generated using AUTO in another programming language. To do this, I wrapped an external call to a function in a different langauge

function derivfunc(X::Vector{Float64},params::Vector{Float64})
    # call to deriv function in other software 
end

derivfunc is not suitable for Automatic differentiation and as such the inputs of the function can only be vectors of real numbers. ( I know this is a slightly contrived example, but the generalised case where automatic differentiation cannot be used in the call to the derivative function seems generic enough)

In terms of the continuation of equilibria, my question is related to #68 . By setting the J = (x,p)->BifurcationKit.finite_differences(z->derivfunc(z,p), x)) in the BifurcationProblem I was able to solve for equilibrium positions. (I have a quick side note here: both #68 and the documentation for BifurcationProblem suggest using BifurcationKit.finiteDifferences. I guess the naming was changed at some point, but the doc string wasn't updated?)

However, I cannot solve for periodic orbits (the continuation of a simple Hopf). Regardless of which method I pick the code will always error on line 243 of Problems.jl

if isnothing(d2F)
    d2F = (x, p, dx1, dx2) -> ForwardDiff.derivative(t -> d1Fad(x .+ t .* dx2, p, dx1), zero(eltype(dx1)))
    d2Fc = (x, p, dx1, dx2) -> BilinearMap((_dx1, _dx2) -> d2F(x, p, _dx1, _dx2))(dx1, dx2)
else
    d2Fc = d2F
end

e.g. it still attempts to calculate the second (and third) derivatives using automatic differentiation. I first tried to 'cheat' this by making d2F not nothing, e.g. (d2F = (x,p,dx1,dx2)->throw(DomainError("Second derivatives not implemented"))) but alas, all the periodic methods I tried require the calculation of d2F. Additionally, I would like to estimate the normal form, which also requires the derivatives.

So my questions:

1. do any periodic methods not require the calculation of higher-order derivatives?
2. can you suggest a function to use to calculate d2F and d3F numerically? (I know AUTO / COCO have internal functions for this, but I'll be honest, I've always found it hard to grasp the proper implementation of the higher order derivatives)

Thank you for any help you can provide.

Cheers,
Fintan

[rveltz]

Hi,

(I have a quick side note here: both #68 and the documentation for BifurcationProblem suggest using BifurcationKit.finiteDifferences. I guess the naming was changed at some point, but the doc string wasn't updated?)

Thank you for pointing this. The issue is based on an older version of BK so that's why the name difference. Where do you see this in the docs?

[rveltz]

However, I cannot solve for periodic orbits (the continuation of a simple Hopf). Regardless of which method I pick the code will always error on line 243 of Problems.jl

Yes, contrary to Auto and Matcont, BK uses the Hopf normal form to have a robust branch switching from a simple Hopf point. That answers your point 1. I could try to add a Basic branch switching ala Auto but that will take some time.

[rveltz]

For the other derivatives, you can perhaps use https://github.com/JuliaDiff/FiniteDifferences.jl

Otherwise, this could work (not tested)

δ = 1e-8
d1Fad(x,p,dx1) = (F(x .+ δ .* dx1, p) .- F(x .- δ .* dx1, p)) /(2δ)
d2F(x, p, dx1, dx2) = (d1Fad(x .+ δ .* dx2, p, dx1) - d1Fad(x .- δ .* dx2, p, dx1) )/(2δ)

[fh9g12]

Thank you for pointing this. The issue is based on an older version of BK so that's why the name difference. Where do you see this in the docs?

I saw it in the doc string of BifurcationProblem. Line 129 of Problems.jl, near the end of the line.

[fh9g12]

Thanks for your suggestions. I did a bit of testing with the problem outlined in https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/tutorials/ode/tutorialsODE/#nmepo

dx1 = [1.0,0,0]
# AD third third derivative
d1Fad(x,p,dx1) = ForwardDiff.derivative(t -> TMvf(x .+ t .* dx1, p), 0.0)
d2F_ad = (x, p, dx1, dx2) -> ForwardDiff.derivative(t -> d1Fad(x .+ t .* dx2, p, dx1), 0.0)
d3F_ad = (x, p, dx1, dx2, dx3) -> ForwardDiff.derivative(t -> d2F_ad(x .+ t .* dx3, p, dx1, dx2), 0.0)

J_ad = ForwardDiff.jacobian(x -> TMvf(x, par_tm), z0)
d2_ad = d2F_ad(z0, par_tm, dx1, dx1)
d3_ad = d3F_ad(z0, par_tm, dx1, dx1, dx1)
res_ad = [J_ad d2_ad d3_ad]

# FD third derivative
δ = 1e-4
J =   (x, p) -> BifurcationKit.finite_differences(z->TMvf(z,p), x)
d2F_fd = (x, p, dx1, dx2) -> 1/(4δ^2)*(TMvf(x .+ δ*dx1 .+ δ*dx2, p) - TMvf(x .+ δ*dx1 .- δ*dx2, p) - TMvf(x .- δ*dx1 .+ δ*dx2, p) + TMvf(x .- δ*dx1 .- δ*dx2, p))
d3F_fd = (x, p, dx1, dx2, dx3) -> 1/(2δ)*(d2F_fd(x .+ δ*dx1, p, dx2, dx3) - d2F_fd(x .- δ*dx1, p, dx2, dx3))

J_fd = J(z0, par_tm)
d2_fd = d2F_fd(z0, par_tm, dx1, dx1)
d3_fd = d3F_fd(z0, par_tm, dx1, dx1, dx1)
res_fd = [J_fd d2_fd d3_fd]

Somewhat unsurprisingly, the results are highly dependent on δ. Anything smaller than 1e-6 and the second derivative starts significantly losing accuracy, and anything smaller than 1e-4 and the third derivative starts losing accuracy. Using δ=1e-2, the normal form at the second hopf matches well for b but not a.

AD

┌─ a = -856.549747865979 - 524.9763996706965im
└─ b = 0.04946305648728988 + 0.1911290780993767im

FD

┌─ a = -7.75693297424067 - 6.581284199549626im
└─ b = 0.04946532901777781 + 0.1911302664401131im

And I can get somewhere close with Periodic orbit simulations.

I'm sure I could get closer with a bit more fiddling (e.g., higher-accuracy finite difference schemes / tweaking continuation parameters), and perhaps I chose a poor model to try this on. But in general, I wouldn't call the process robust.

Overall, I think my conclusion is that if the code is not suitable for AD, stick to other tools with less reliance on the higher-order derivatives (you wouldn't know of an AUTO julia wrapper would you? :) ). However, I shall continue to explore using BifurcationKit for code written purely in Julia.

Cheers for your help,
Fintan

[rveltz]

Hi,

You probably need to scale δ = 1e-8 for first order, δ = 1e-4 for second and δ=1e-2 for third. But this is not used for the continuation of periodic orbits. It is only used to compute the normal form.

Once the initial guess is computed, the code is pretty similar to Auto and Matcont.

I also wrote a branch switching method that does not rely on AD.

For your plot, I am surprised. The Floquet coefficients should only derail on the way to the homoclinic point. You may be using an old version of the code.

I will add a backend for full FD when I have some time

[fh9g12]

Once the initial guess is computed, the code is pretty similar to Auto and Matcont.

Ah ok, perhaps I should hold back on making any judgement yet.

For your plot, I am surprised. The Floquet coefficients should only derail on the way to the homoclinic point. You may be using an old version of the code.

I'm using version 0.4.6, which I think is the latest release.

I also wrote a branch switching method that does not rely on AD.
...
I will add a backend for full FD when I have some time

I'll be interested in trying these updates in the future. My current research goal is to run continuation on higher-dimensional aeroelastic models, which may not always play nicely with AD...

This wasn't a rabbit hole I thought I would go down (this was meant to be an easy first problem), but it's been a good learning exercise for Julia. If it's of any use as a test case, my adaptation of the example is below. The 2nd order derivative uses δ = 1e-4 and the third order δ = 1e-2.

For a nominal test case, the finite difference at δ = 1e-2 compares well to AD, and the estimation of the normal form also matches well with AD

┌─ a = -856.5497693138566 - 524.9764134777083im
└─ b = 0.049469181253591295 + 0.19113570713773032im

but the PO continuation currently fails when compared to the AD solution

Once again, cheers for your active communication!

using Revise, Plots
using BifurcationKit
using ForwardDiff
const BK = BifurcationKit

function TMvf(z, p, t = 0)
	(;J, α, E0, τ, τD, τF, U0) = p
	E, x, u = z
	SS0 = J * u * x * E + E0
	SS1 = α * log(1 + exp(SS0 / α))
	dz = zeros(typeof(z[1]),3)
	dz[1] = (-E + SS1) / τ
	dz[2] =	(1.0 - x) / τD - u * x * E
	dz[3] = (U0 - u) / τF +  U0 * (1.0 - u) * E
	dz
end
# parameter values
par_tm = (α = 1.4, τ = 0.013, J = 3.07, E0 = -2.0, τD = 0.20, U0 = 0.3, τF = 1.5)

# initial condition
z0 = [0.238616, 0.982747, 0.367876]

dx = [1.0,0,0]
# estimate 3rd deriv with AD
d1Fad(x,p,dx1) = ForwardDiff.derivative(t -> TMvf(x .+ t .* dx1, p), 0.0)
d2F_ad = (x, p, dx1, dx2) -> ForwardDiff.derivative(t -> d1Fad(x .+ t .* dx2, p, dx1), 0.0)
d3F_ad = (x, p, dx1, dx2, dx3) -> ForwardDiff.derivative(t -> d2F_ad(x .+ t .* dx3, p, dx1, dx2), 0.0)
println(d3F_ad(z0, par_tm, dx, dx,dx))

# estimate with finite differnece
J =   (x, p) -> BifurcationKit.finite_differences(z->TMvf(z,p), x)
d2F_fd = (x, p, dx1, dx2, δ) -> 1/(4δ^2)*(TMvf(x .+ δ*dx1 .+ δ*dx2, p) - TMvf(x .+ δ*dx1 .- δ*dx2, p) - TMvf(x .- δ*dx1 .+ δ*dx2, p) + TMvf(x .- δ*dx1 .- δ*dx2, p))
d3F_fd = (x, p, dx1, dx2, dx3, δ) -> 1/(2δ)*(d2F_fd(x .+ δ*dx1, p, dx2, dx3, δ) - d2F_fd(x .- δ*dx1, p, dx2, dx3, δ))
println(d3F_fd(z0, par_tm, dx, dx,dx,1e-2))

# Bifurcation Problem
prob_fd = BifurcationProblem(TMvf, z0, par_tm, (@optic _.E0);
	record_from_solution = (x, p; k...) -> (E = x[1], x = x[2], u = x[3]),
	J =J, 
	d2F = (x, p, dx1, dx2) -> d2F_fd(x, p, dx1, dx2, 1e-4), 
	d3F = (x, p, dx1, dx2, dx3) -> d3F_fd(x, p, dx1, dx2, dx3, 1e-2))

# continuation options
opts_br = ContinuationPar(p_min = -2.0, p_max = -1.)

# continuation of equilibria
br_fd= continuation(prob_fd, PALC(), opts_br)

scene = plot(br_fd, plotfold=false, markersize=4, legend=:topleft)

#estimate normal form of 2nd Hopf
nf_fd = get_normal_form(br_fd,4)
println(nf_fd)

# run PO continuation
opts_po_cont = ContinuationPar(opts_br, dsmax = 0.1, ds = 0.004, dsmin = 1e-4,
	max_steps = 80, tol_stability = 1e-6)

args_po = (	record_from_solution = (x, p; k...) -> begin
	xtt = get_periodic_orbit(p.prob, x, p.p)
	return (max = maximum(xtt[1,:]),
			min = minimum(xtt[1,:]),
			period = getperiod(p.prob, x, p.p))
end,
# we use the supremum norm
normC = norminf)
    
br_fd_po = continuation(br_fd, 4, opts_po_cont,
PeriodicOrbitOCollProblem(50, 4; meshadapt = true);args_po...)
plot(br_fd, br_fd_po)

[rveltz]

I'll be interested in trying these updates in the future. My current research goal is to run continuation on higher-dimensional aeroelastic models, which may not always play nicely with AD...

It is very difficult to write code that has no bug in particular when computing the normal forms. Using AD removes a lot of difficulties but now that it has settled I should relax that need.

In BK you really only need higher order derivatives to compute predictors for branch switching. But I will soon allow to use a very simple predictor (ala Auto - Matcont)

$pred = x_0+amp * eigenvec$

Now for your high dimensional aeroelastic models, you will be able to do the same as the TMModel but with Jacobian-Free methods (which are not in Auto or MatCont).

[rveltz]

Ah ok, perhaps I should hold back on making any judgement yet.

No problem. Your example shows me that my treatment of stability of the periodic orbits may be too simple. Maybe Auto or Matcont is doing something clever with the Floquet coefficients. In your plot, you see a lot of spurious bifurcations. This is because the trivial multiplier $1$ is not exactly one numerically (say 1 + 1e-6), and when it is above one, a (spurious) bifurcation is reported.

I should note that the (numerical) value of the Floquet coeffs is the same as Auto's or Matcont's.

[fh9g12]

Now for your high dimensional aeroelastic models, you will be able to do the same as the TMModel but with Jacobian-Free methods (which are not in Auto or MatCont).

Ah, interesting. I've not used Jacobian free methods before, so I'll definitely give them a go at some point!

Maybe Auto or Matcont is doing something clever with the Floquet coefficients

I'm not 100% sure I'm in the right part of the source code, but I wouldn't call it clever. it seems to be floq < 1+tol

https://github.com/auto-07p/auto-07p/blob/8b6cf/src/periodic.f90#L2029-L2039

There's also a bit of a discussion in the floquet multiplier section about methods to minimise roundoff errors

https://github.com/auto-07p/auto-07p/blob/beb2fb79c9bb/src/floquet.f90#L144-L167