Initial version of Lie-Euler

[mateuszbaran]

I've stared working on a Lie-Euler solver. A few things to work out:

• [x] I don't exactly understand the correct use of mutable caches. Examples I've found are very confusing, filling it with things that aren't used by perform_step! or duplicating information from integrator.
• [x] Whether ODEProblem is OK for this or do we actually need ExplicitManifoldODEProblemType? I don't know yet.
• [x] DiffEq returns a solution that performs interpolation when asked for a value outside of the timesteps. The default interpolator isn't correct for manifolds and we may have to replace it with our own.
• [x] More testing.
• [x] Links and badges.

[mateuszbaran]

An accidental cool looking spiral :slightly_smiling_face:

using GLMakie, Makie, LinearAlgebra

n = 10

θ = [0;(0.5:n-0.5)/n;1]
φ = [(0:2n-2)*2/(2n-1);2]
x = [cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
y = [sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
z = [cospi(θ) for θ in θ, φ in φ]

function f2(x, y, z)
    return exp(x) * cross([x, y, z], [1.0, 0.0, 0.0]) + exp(y) * cross([x, y, z], [0.0, 1.0, 0.0])
end

tans = f2.(vec(x), vec(y), vec(z))
u = [a[1] for a in tans]
v = [a[2] for a in tans]
w = [a[3] for a in tans]

scene = Scene();

#surface(x, y, z)
arr = Makie.arrows(
    vec(x), vec(y), vec(z), u, v, w;
    arrowsize = 0.1, linecolor = (:gray, 0.7), linewidth = 0.02, lengthscale = 0.1
)

using ManifoldDiffEq, OrdinaryDiffEq

A = ManifoldDiffEq.ManifoldDiffEqOperator{Float64}() do u, p, t
    return f2(u...)
end
prob = ODEProblem(A, [0.0, 1.0, 0.0], (0, 20.0))
alg = ManifoldDiffEq.ManifoldLieEuler(Sphere(2), ExponentialRetraction())
sol1 = solve(prob, alg, dt = 0.01)

Makie.lines!([u[1] for u in sol1.u], [u[2] for u in sol1.u], [u[3] for u in sol1.u]; linewidth = 10, color=:red)

[mateuszbaran]

Interesting, it should probably be a closed loop and the spiral is a solver error.

[kellertuer]

Do we have a way to get this more stable?

[mateuszbaran]

Smaller timestep helps but probably a higher-order method will be a better choice. I'll add some RKMK variants once I'm happy with the functionality of Lie-Euler.

[mateuszbaran]

CG2 and CG3 should be even easier and work on this out-of-the box.

[mateuszbaran]

@kellertuer I need Hermite interpolation on manifolds to interpolate solutions. Where should it be put? I've found this paper about it: https://arxiv.org/abs/1908.05875 but it requires some variant of differential of Riemannian logarithm, which is currently a part of Manopt. Should I contribute it to Manopt, or maybe we should have another package for that?

[kellertuer]

Oh, that gets interesting. I am not sure where to put that, but yes, currently I have it in Manopt to compute differentials.

[kellertuer]

We nearly used up our private repo testing, just as a side note :)

[mateuszbaran]

Then we can just open the repository :slightly_smiling_face: . I think it's OK to open even now?

[kellertuer]

Sure, we can open that now; we are far from registering this, but still we can do this openly, no problem.

[codecov[bot]]

Codecov Report

:exclamation: No coverage uploaded for pull request base (main@97530a1). Click here to learn what that means. The diff coverage is n/a.

@@           Coverage Diff           @@
##             main       #1   +/-   ##
=======================================
  Coverage        ?   92.18%           
=======================================
  Files           ?        7           
  Lines           ?      192           
  Branches        ?        0           
=======================================
  Hits            ?      177           
  Misses          ?       15           
  Partials        ?        0           

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 97530a1...9484367. Read the comment docs.

[mateuszbaran]

I think I've got CG2 right finally and it seems to work great (the green curve; the red one is Lie-Euler with the same time step).

[kellertuer]

Nice! I still did not yet have the time to look at the code.

[mateuszbaran]

Thanks! The code isn't particularly finished yet, I've just today figured out the essential parts, so probably it's better that you didn't spend time reading the old versions :slightly_smiling_face: .

[mateuszbaran]

I tried implementing CG3 but it turns out to be very unstable on my test problem. @kellertuer could you check if my approach here is correct? I'm thinking mostly about this: https://github.com/JuliaManifolds/ManifoldDiffEq.jl/pull/1/files#diff-27108084e34e3b6b2e3fd70d737b0dccc6517a57b7bd1b792dbcf3d2bb37aacdR204-R227 . Vector transports are supposed to express the flow along frozen coefficient fields but I'm not sure if I transport the intermediate stages correctly.

[mateuszbaran]

RKMK4 seems to work really well though.

[mateuszbaran]

TODO: update citation keys to include all authors.

[kellertuer]

I like the new writing with M makes it more readable :)

[kellertuer]

Thanks for switching for example to retraction_method. I closed and opened it so that the documenter CI will hopefully pick up the label.

[mateuszbaran]

• Do we want to do a short general intro on the setup of differential equations on manifolds? Might be good to establish the notation

I'd like to do that but maybe not necessarily in this PR.

• Do we want to provide a list/overview of available solvers? In the long run maybe even a page per solver would be nice?

For now I think one page for Lie group solvers and one for frozen coefficient solvers should be enough, we can always add more details later. I'll add a list of solvers at the top of each page.

If I see correctly we currently have 3 (or 4?) solvers already – a plain forward Euler. CG2, CG3, RKMK4? For all it would be nice to know where they are applicable; for example whether on all manifolds or only on Lie.

We have 5: frozen Euler, CG2, CG3, Lie Euler and RKMK4. I'll add some details regarding applicability to docs but in short frozen coefficients solvers work on any manifold and Lie group solvers need an explicit action of a Lie group on a manifold in question. There are some conditions this action needs to satisfy but I'm not entirely sure in what cases these solvers would stop working.

[kellertuer]

Ok, so just for the last point – if even I do not see that in the docs – maybe we should at least give the overview somewhere.

[mateuszbaran]

I have no idea what the problem with CG3 is, I've been checking it and it seems according to papers but it's performing very poorly on my test problem (blue curve is CG3, green is CG2, red is Lie-Euler, dt is 0.05: ). It's bad even on R^3 so I'll probably just skip this one.

[mateuszbaran]

There might be a problem renaming things like p or u: the integrator interface has fixed names and if we change them, things are going to not work that well:

function perform_step!(integrator, cache::CG2Cache, repeat_step = false)
    @unpack t, dt, uprev, u, f, p, alg = integrator

[kellertuer]

Hm sad. Then we have to stick with them; but that means that unpack takes integrator.u and puts it into u? Ok. So we stick with them if that's what they want.

[kellertuer]

Can we still use X and p in our caches?

[mateuszbaran]

that means that unpack takes integrator.u and puts it into u

Yes, that's how I understand it.

Can we still use X and p in our caches?

Yes, but I think for p it may be more confusing than helpful. Using X instead of k should be fine.

[kellertuer]

Oh, I should maybe have done that here in changes and not at the commit, sorry.

What about the first and fourth task?

• For the first we should look whether the constructor can be shortened?
• Which tests are missing?
• We should add a Readme with a few shields/badges and links to ManifoldsBase.jl/Manifolds.jl/DiffEq.jl.

[mateuszbaran]

Oh, I should maybe have done that here in changes and not at the commit, sorry.

No problem, commenting on a commit is also fine.

• For the first we should look whether the constructor can be shortened?

I think the key thing of the first point is solved, I've started using these caches. Which constructor would you like to shorten?

• Which tests are missing?

I think I'd at least like to check that on R^n our solutions match DiffEq.jl methods with the same tableau. That would probably be enough for now, though we could also for example check the empirical order of convergence on some examples.

• We should add a Readme with a few shields/badges and links to ManifoldsBase.jl/Manifolds.jl/DiffEq.jl.

Right, good point.

[kellertuer]

• For the first we should look whether the constructor can be shortened?

I think the key thing of the first point is solved, I've started using these caches. Which constructor would you like to shorten?

I marked one occurrence above which I found a little long (but maybe there is no way around?)

• Which tests are missing?

I think I'd at least like to check that on R^n our solutions match DiffEq.jl methods with the same tableau. That would probably be enough for now, though we could also for example check the empirical order of convergence on some examples.

for now the Rn checks should be fine.

• We should add a Readme with a few shields/badges and links to ManifoldsBase.jl/Manifolds.jl/DiffEq.jl.

Right, good point.

I added such a Readme now.

[mateuszbaran]

I'll review this once again, I haven't added yet test for Lie group action solvers -- but indeed that should be more or less enough :slightly_smiling_face: .

[kellertuer]

I refined the docs a little bit, but some docstrings are still a little minimalistic. bFor a first version this might be fine.

[mateuszbaran]

I've added tests for Lie group methods and I think merging this should be fine now.

[kellertuer]

Well, then. Go for it :)