Polish the Differential Interfaces

[kellertuer]

This PR fixes #428, should also resolve #88. To some extend it might also resolve #42 together with the last PR? At least in the embedding Zygote can be used easily now; for intrinsic methods, there is even in theory still more work to do.

• with the introduction of properly solving the Jacobians, there is a few get_vector! functions missing, that are now required (yields a few ErrorExceptions where MethodErrors were exempted and a few missing implementation errors
• Test Coverage has to be checked.

This PR also makes the dependency on FiniteDifferences optional, since we now have enough to choose from.

[kellertuer]

edit: To avoid confusion with git commit comments from other branches I cleared this one by squashing the work until now into one commit.

edit2: I then had to do the same strange merge conflict from yesterday again. Squashing and merging is really strange sometimes. I will not merge another dev branch ever again ;)

[mateuszbaran]

It looks like merging branches can be not only confusing for humans but for github as well 😉 .

I think after this we can mark #42 as resolved; the intrinsic approach is a somewhat separate issue, tracked in JuliaManifolds/ManifoldDiff.jl#28 IIUC.

[kellertuer]

Cool, I think so, too. Then this PR fixes even 3 issues. Nice!

[kellertuer]

I noticed one challenge when I tried fixing the first errors.

Currently we compute the new local metric jacobian here

Manifolds.jl/src/manifolds/MetricManifold.jl

Lines 471 to 480 in b2f2674

	∂g = reshape(
	_jacobian(
	c -> local_metric(M, retract(M, p, get_vector(M, p, c, B), retraction), B),
	p,
	backend,
	),
	d,
	d,
	d,
	)

using a default retraction.

This has a limit for the case where we solve the ODE for the geodesic/exp at

Manifolds.jl/src/manifolds/ConnectionManifold.jl

Lines 164 to 172 in b2f2674

	@decorator_transparent_fallback function exp!(M::AbstractConnectionManifold, q, p, X)
	tspan = (0.0, 1.0)
	A = get_default_atlas(M)
	i = get_chart_index(M, A, p)
	B = induced_basis(M, A, i, TangentSpace)
	sol = solve_exp_ode(M, p, X, tspan, B; dense=false, saveat=[1.0])
	n = length(p)
	return copyto!(q, sol.u[1][(n + 1):end])
	end

which internally calls solve_exp_ode which internally requires local_jacobian which .... uses the default retraction (exp) and the circle repeats finishing in s stack overflow.

We have to avoid, that (for the case the local jacbian is used to compute the exp problem) that the ExponentialRetraction is used.

[mateuszbaran]

Hm, that's an interesting problem. My guess would be that the only way to make this work is working in a chart, that is local_metric_jacobian can would only be defined for bases induced by charts, where you don't have to use retract.

[kellertuer]

uff. I would prefer checking (before calling local_jacobian that the (explicit, non-default) retraction is not the exponential retraction in solve_ode_exp. That would be maybe be the easier fix?

[mateuszbaran]

uff. I would prefer checking (before calling local_jacobian that the (explicit, non-default) retraction is not the exponential retraction in solve_ode_exp. That would be maybe be the easier fix?

I don't see any retraction in solve_exp_ode, it just has christoffel_symbols_second which calls the problematic Jacobian calculation. You could check, though, that the retraction local_metric_jacobian gets is not the fallback case -- at least after refactoring the code of exp!(M::AbstractConnectionManifold, q, p, X) a bit...

[mateuszbaran]

That is, exp!(M::AbstractConnectionManifold, q, p, X) would be an ODERetractionMethod, which would be the default retraction on AbstractConnectionManifold, and local_metric_jacobian could check if its retraction is different from ODERetractionMethod.

[kellertuer]

uff. I would prefer checking (before calling local_jacobian that the (explicit, non-default) retraction is not the exponential retraction in solve_ode_exp. That would be maybe be the easier fix?

I don't see any retraction in solve_exp_ode, it just has christoffel_symbols_second which calls the problematic Jacobian calculation. You could check, though, that the retraction local_metric_jacobian gets is not the fallback case -- at least after refactoring the code of exp!(M::AbstractConnectionManifold, q, p, X) a bit...

Of course not. Until now it used q+ ... , that is what we are changing, so we are introducing this difficulty here. solve_ode_exp will have a retraction parameter and yes, we have to refactor here also the exp! .

[kellertuer]

That is, exp!(M::AbstractConnectionManifold, q, p, X) would be an ODERetractionMethod, which would be the default retraction on AbstractConnectionManifold, and local_metric_jacobian could check if its retraction is different from ODERetractionMethod.

That is a very good idea since it only approximates exp, that solves the problem to that extend.

[mateuszbaran]

uff. I would prefer checking (before calling local_jacobian that the (explicit, non-default) retraction is not the exponential retraction in solve_ode_exp. That would be maybe be the easier fix?

I don't see any retraction in solve_exp_ode, it just has christoffel_symbols_second which calls the problematic Jacobian calculation. You could check, though, that the retraction local_metric_jacobian gets is not the fallback case -- at least after refactoring the code of exp!(M::AbstractConnectionManifold, q, p, X) a bit...

Of course not. Until now it used q+ ... , that is what we are changing, so we are introducing this difficulty here. solve_ode_exp will have a retraction parameter and yes, we have to refactor here also the exp! .

OK, I thought that's not the issue we are solving here. It seems that exp!(M::AbstractConnectionManifold, q, p, X) could either work in a chart (and then no retraction is needed as it just does +) or use ManifoldDiffEq.jl, and both method could even coexist. Working in charts means a wider selection of ODE solvers, while using ManifoldDiffEq.jl means no problems caused by charts.

[kellertuer]

uff. I would prefer checking (before calling local_jacobian that the (explicit, non-default) retraction is not the exponential retraction in solve_ode_exp. That would be maybe be the easier fix?

I don't see any retraction in solve_exp_ode, it just has christoffel_symbols_second which calls the problematic Jacobian calculation. You could check, though, that the retraction local_metric_jacobian gets is not the fallback case -- at least after refactoring the code of exp!(M::AbstractConnectionManifold, q, p, X) a bit...

Of course not. Until now it used q+ ... , that is what we are changing, so we are introducing this difficulty here. solve_ode_exp will have a retraction parameter and yes, we have to refactor here also the exp! .

OK, I thought that's not the issue we are solving here. It seems that exp!(M::AbstractConnectionManifold, q, p, X) could either work in a chart (and then no retraction is needed as it just does +) or use ManifoldDiffEq.jl, and both method could even coexist. Working in charts means a wider selection of ODE solvers, while using ManifoldDiffEq.jl means no problems caused by charts.

Hm, then I though about a completely different way of solving this. Remove theexp! you mentioned and do that on a Retraction with an ExponentialODEApproximationRetraction that has a retraction internally (and with the right basis that retraction could be +).

but for your two approaches I have not yet understood how you mean to solve those, neither in theory for the charts with + not in practice for both.

[mateuszbaran]

Could you explain how would you like to use a retraction in ExponentialODEApproximationRetraction?

[kellertuer]

I would remove the default fallback of exp! to the ODE. Reason: The ODE is only a numerical approximation and we avoid the ambiguity we are currently running into.

A retraction of type ExponentialODEApproximationRetraction would have a field, which retraction (and maybe also with basis?) to use and this would (in the constructor) be checked that it can not be the ExponentialRetraction (since that is what is introducing the current infinite loop).

[mateuszbaran]

I would remove the default fallback of exp! to the ODE. Reason: The ODE is only a numerical approximation and we avoid the ambiguity we are currently running into.

Yes, I agree that this is a good idea.

A retraction of type ExponentialODEApproximationRetraction would have a field, which retraction (and maybe also with basis?) to use and this would (in the constructor) be checked that it can not be the ExponentialRetraction (since that is what is introducing the current infinite loop).

I don't quite see yet how you want to rewrite solve_exp_ode. Anyway, I think the in-chart method could be said to use the following retraction:

struct ChartRetraction{TA<:AbstractAtlas,TI} <: AbstractRetractionMethod
    A::TA
    i::TI
end

function retract(M, p, X, cr::ChartRetraction)
    pc = get_parameters(M, cr.A, cr.i, p)
    Xc = get_coordinates(M, p, X, InducedBasis(TangentSpaceType, A, i))
    return get_point(M, cr.A, cr.i, pc + Xc)
end

The problem is that it's currently only implemented for StereographicAtlas.

[kellertuer]

Ah that looks also nice and would be the inner retraction then, sure.

[mateuszbaran]

Notice though that when you write a solver for geodesic using that retraction you can eliminate calling get_parameters and get_point on every step (and only do so at the very beginning at at the very end), and in doing so you effectively end up with my "geodesic in chart" idea.

[kellertuer]

Yes, I see. One could implement a special solver for this case (dispatching on ExponentialODEApproximationRetraction{ChartRetraction}.

[mateuszbaran]

Yes, and I think that solver would be the easiest one to write (and the closest to the current implementation).

[kellertuer]

Could you try to propose that? Because I think the generic one should be relatively close to what is done now actually (with the retraction inside instead of the +, sure your approach would really just use + in _jacobian call).

[mateuszbaran]

I think I could write that.

Also, which + do you want to replace with retraction? Currently the addition (for performing a step) is done inside OrdinaryDiffEq.jl.

[kellertuer]

Before we started this PR, there was basically a q+... (yes within OrdinaryDiffEq), which we now shifted accordingly into the coordinates of a vector; this in principle (since the function we work on in f \circ retr (circ get_vector) is now replaced by the retraction. So I do not want to, I already did that; just due to the loop that the retr causes an inner soilve_ode – again and again – it is just failing horribly, but hour combined solution ideas should fix that :)

[kellertuer]

ON my way of fixing the ODE thingy, I noticed that both connection and metric decorator implement the jacobian-things (local metric, first, second), but while Metric just calls an ODE for the local metric jacobian and continues from there, the ConnectionManifold computes the second jacobians with an ODE (but neither first nor local metric jacobina9: Is that an inconsistency or intended?

[mateuszbaran]

The test was expecting local_metric to be defined for DefaultOrthonormalBasis. Removing that method of inverse_local_metric didn't break anything so I removed it (I would expect the default implementation to work correctly).

[kellertuer]

I ran into a next problem when starting the own ODE tests:

In the version before this version, the ode-parameters were stored as points in the embedding, but now I thought it is necessary to reduce the size to the manifold_dimension? This is what the code currently does, but then there is a size mismatch when initialising the ODE

Manifolds.jl/src/differentiation/ode.jl

Lines 15 to 16 in a4a2716

	u0[iv] .= X
	u0[ix] .= p

so we still might habe to revise this to fit.

Edit: Taking the representation size should work for now, but I again get stuck on local metrics. They are really something to hunt me with. Implementing it for the metric manifold did not work, implementing it for the manifold itself did not work, tried both with less or more parametrised bases – never tries to “find/catch” - I seem generally unable to “track” down the right signature for bases to be found in dispatch on them with all their parameters.

[mateuszbaran]

In the version before this version, the ode-parameters were stored as points in the embedding, but now I thought it is necessary to reduce the size to the manifold_dimension?

I'm not sure. I only know how to do this ODE solving either in a chart or through ManifoldDiffEq.jl.

tried both with less or more parametrised bases – never tries to “find/catch” - I seem generally unable to “track” down the right signature for bases to be found in dispatch on them with all their parameters.

I'm fixing it, the main problem I think was that you forgot Manifolds. to add methods to the function instead of making a new function.

[mateuszbaran]

There are still errors but much further along (local_metric_jacobian is broken).

[kellertuer]

Thanks for the help, I think in my many tries I had that variant, too, but I really seem to struggle a lot with things involving metric/basis combined here, especially which parameters to specify for the basis (I feel like bases really have many parameters, though it is also just 3 or something?).

Sure there are more errors down the road, I will slowly work towards that.

Is there a reason I can not ask the sphere for the metric? I had hoped it just has a metric tensor (namely the one from the embedding restricted)?

[mateuszbaran]

Thanks for the help, I think in my many tries I had that variant, too, but I really seem to struggle a lot with things involving metric/basis combined here, especially which parameters to specify for the basis (I feel like bases really have many parameters, though it is also just 3 or something?).

Bases really have only two meaningful parameters: whether it's a real or complex basis and whether it's a basis of the tangent space or some other space. Probably the length of <:ManifoldsBase.TangentSpaceType makes it seem like something much more involved and it would be nice to figure out a shorter notation for bases of tangent spaces.

Is there a reason I can not ask the sphere for the metric? I had hoped it just has a metric tensor (namely the one from the embedding restricted)?

The only reason is that no one has implemented that function yet 😄 . Metric tensor from the embedding can't work because:

1. On an N-sphere it's an NxN matrix, not an (N+1)x(N+1).
2. There is by default no relation established between bases of embedded manifolds and bases of the embedding.
   Both problems could be solved but that's relatively complicated. And it's much easier to just write a special case for sphere.

[kellertuer]

Bases really have only two meaningful parameters: whether it's a real or complex basis and whether it's a basis of the tangent space or some other space. Probably the length of <:ManifoldsBase.TangentSpaceType makes it seem like something much more involved and it would be nice to figure out a shorter notation for bases of tangent spaces.

My problem is I tried (a) TestSphere with the Basis (b) TestSphere with basis and one Parameter (c) two parameters (d) MetricManifold with basis (no parameters) e) one parameter, f) two parameters – ... got annoyed and left. That is my usual workflow. 6-8 tries, none works -> you fix it and then it looks like one that I already tried. I do not yet know what I am doing wrong in my near-infinite number of tries.

The only reason is that no one has implemented that function yet 😄

Ah these lazy developers again 😄

Metric tensor from the embedding can't work because:

1. On an N-sphere it's an NxN matrix, not an (N+1)x(N+1).
2. There is by default no relation established between bases of embedded manifolds and bases of the embedding.
   Both problems could be solved but that's relatively complicated. And it's much easier to just write a special case for sphere.

no in general that is too complicated but I might then add the local metric for the sphere just to reuse it here for the test sphere.

[mateuszbaran]

My problem is I tried (a) TestSphere with the Basis (b) TestSphere with basis and one Parameter (c) two parameters (d) MetricManifold with basis (no parameters) e) one parameter, f) two parameters – ... got annoyed and left. That is my usual workflow. 6-8 tries, none works -> you fix it and then it looks like one that I already tried. I do not yet know what I am doing wrong in my near-infinite number of tries.

I think the reasonable workflow would be:

1. When doing metric stuff, always add methods for the manifold wrapped in MetricManifold -- unless there is a clear reason for writing a method for the wrapper-less case.
2. When doubt, add more parameters to the basis.
3. Check simple cases (i.e. check that you are actually able to directly call the method you want to be called and it produces the desired result).
4. If problems persist, let me know 🙂 .

I might then add the local metric for the sphere just to reuse it here for the test sphere.

That's a good idea 👍

[kellertuer]

Actually, we have the local metric already, since the Sphere has the EuclideanMetric inherited from its embedding and for any manifold with the Euclidean metric, we already return the diagonal matrix accordingly. So we (both) were just missing to call the local Metric with the MetricManifold(Sphere(n), EuclideanMetric())...edit: Ah a.ain only for the induced Basis. Then I will add it still.

[kellertuer]

Somewhere we still have some dimension mixup (i.e. whether we are in coordinates at a place or in the representation). Currently I somehow end up in a get_vector on the sphere S2 with a coordinate vector of length 3, which is definetlz wrong. We might have to check that thoroughly again (until now – also with the testSphere over in Metric coordinate and vector dimensions always agreed I think?)...

Edit computed the jacobian in a basis, so we have to convert dx as well, now all dimensions are correct, just the result is still wrong.

Edit2: I know why: The extrinsic (embedded) Jacobian computation is now wrong in the generic setting. I think I have an Idea how to fix this, but I need some time to think about it. The idea is (a) in the embedded case: embed -> jacobian -> project/coordinates and (b) implement explicit formula in other cases.

[kellertuer]

I am a little stuck (maybe I thought about it too much or something) – somehow I can not manage to write the solve_exp_ode correctly, currently I either end up with two things . Let d be the manifold dimension and n be the representation size (for a matrix represented in vectors, otherwise reshape to vectors so d and n are numbers) Then I am able to

a) write a variant that tries to work in a too large space (n) where the @einsum complains (which is totally correct, since that is of size d\times d and the vectors I would then multiply it with in the ensue are n and n

b) write a variant that works in the right space, then I am stuck in how to initialise the p (just zero?) and how to fill that part of the u anyways.

And in general I am not sure in which of the variants to work, since we should have a version (i.e. in some chart) where we can do a “Euclidean” ode, i.e. one where we have a plus right?

I also just tried to write the fallback for exp! directly on the AbstractManifold level, since for the default metric (or default connection, too) the Jacobian would just be implemented for Sphere for example, right?

[mateuszbaran]

a) write a variant that tries to work in a too large space (n) where the @einsum complains (which is totally correct, since that is of size d\times d and the vectors I would then multiply it with in the ensue are n and n

Yes, in this variant you need a different way of calculating the tangent vector for ODE and, unfortunately, I don't currently see a good way to do this.

b) write a variant that works in the right space, then I am stuck in how to initialise the p (just zero?) and how to fill that part of the u anyways.
And in general I am not sure in which of the variants to work, since we should have a version (i.e. in some chart) where we can do a “Euclidean” ode, i.e. one where we have a plus right?

The variant I've written in #431 does exactly that: p is initialized based on its position in a chart (it's 0 if you take for example a retraction chart at a point but you can take a different one and it's not zero then) and a normal "Euclidean" ODE is solved by OrdinaryDiffEq.jl.

I also just tried to write the fallback for exp! directly on the AbstractManifold level, since for the default metric (or default connection, too) the Jacobian would just be implemented for Sphere for example, right?

I think we may wait with that, currently the preferred way of working in Manifolds.jl is to implement explicit retractions, and the support for ODE-based retractions is very experimental. I'd wait with making it the default for non-metric manifolds until we have it at least somewhat tested and benchmarked on a simple problem (such as geodesics on an ellipsoid).

[kellertuer]

Thanks for the help. I will take the approach (a) back to what we had, just generalise it to (isometrically) embedded sub manifolds where I think I have an idea (and then getting to finish this PR and just rewriting/reorganising the tests a little bit).

For the generic fallback my idea stems more from the fact that for the default metrics we “strip away” the metric part anyways (in most cases), since we assume to have the implicit (first/default) metric to be implemented without the metric decorator.

[mateuszbaran]

Thanks for the help. I will take the approach (a) back to what we had, just generalise it to (isometrically) embedded sub manifolds where I think I have an idea (and then getting to finish this PR and just rewriting/reorganising the tests a little bit).

OK, I'll wait until you try your idea 🙂 .

For the generic fallback my idea stems more from the fact that for the default metrics we “strip away” the metric part anyways (in most cases), since we assume to have the implicit (first/default) metric to be implemented without the metric decorator.

I think that could work but you'd still need ManifoldDiffEq.jl instead of OrdinaryDiffEq.jl, right?

[kellertuer]

Currently 26 tests fail in metric, and I have no clue why, because I do not understand what they test, since there is no explanations around and code without comments is a little hard to reverse-engineer.

So I am not sure how to continue here – revert all the ODE stuff and just keep the renaming of the embedded AD stuff? Omit this PR totally? The last part would be a little sad since I spent really a lot of time on this. Even the first is a little sad, but I fear I do not understand enough of the ODE and how to do it in a basis or the embedding to fix it here.

[mateuszbaran]

Currently 26 tests fail in metric, and I have no clue why, because I do not understand what they test, since there is no explanations around and code without comments is a little hard to reverse-engineer.

Looking at that code, I would be surprised if it made sense to you 😅 . It has a few things done very much not right.

OK, so this is what these tests supposed to do:

1. We take a scaled sphere with the default atlas (IIRC it's the exponential retraction atlas):

g = TestSphericalMetric()
M = MetricManifold(Sr, g)
A = Manifolds.get_default_atlas(M)

2. We take a point on that sphere and a chart that contains that point:

p = vtype([θ, ϕ])
chart_p = get_chart_index(M, A, p)
B_p = induced_basis(M, A, chart_p, TangentSpace)

This is somewhat incorrect as p is not a point in the correct representation, I must've updated it incorrectly. p should be just a normal point, like [1.0, 0.0, 0.0]. For some reason local_metric for the scaled sphere is written in a way to take this into account but that's not how local_metric should work. So it needs to be rewritten to a) accept point in the embedding representation and b) we should make sure that local_metric actually returns the coefficients of metric in the basis indicated in the argument.

This leads us to reconsider what point (1) does, and thus we should use the chart corresponding to spherical coordinates, not a chart from the default atlas. It's almost as if someone only cared about making these tests pass 😅 . Sorry.

3. After the chart confusion is cleared, we encounter the failing tests. They are supposed to test if the coefficients of these various tensors are correct (expressed in the spherical coordinates chart).

Perhaps the mix between these two charts breaks your new code but I'm not entirely sure. I can at least fix the old tests so that they make sense but you need to give me a few days.

So I am not sure how to continue here – revert all the ODE stuff and just keep the renaming of the embedded AD stuff? Omit this PR totally? The last part would be a little sad since I spent really a lot of time on this. Even the first is a little sad, but I fear I do not understand enough of the ODE and how to do it in a basis or the embedding to fix it here.

The renaming is certainly fine. The ODE part is very tricky and I think we should be careful here -- it'd very likely take quite a long time to polish it.

[kellertuer]

Thanks for the heads up. What do you think about maybe splitting this into two PRs then?

I could make a new branch, cherry pick the Renaming commits (or take the lines over), and we do a second one to clear the ODE stuff (maybe on a separate PR maybe on your PR for the ODE in charts)?

I would love to have the differentiation part (without ODE) finished to also finish my tutorial on the gradient conversion. I might take a break from the ODE then, though, to first look into the Hessian conversion. But that said, of course there is absolutely no problem if it takes a while :)

[mateuszbaran]

Sure, splitting it into two PRs is a good idea. One PR for renaming, then I can continue my PR to have the geodesics in charts, and then you can work on your approach in another PR.

I would love to have the differentiation part (without ODE) finished to also finish my tutorial on the gradient conversion. I might take a break from the ODE then, though, to first look into the Hessian conversion. But that said, of course there is absolutely no problem if it takes a while :)

Sure, the differentiation part should be simple to finish. Hessians (and, probably more importantly, Hessian-vector products) may be an interesting thing to do after that 🙂 .

[kellertuer]

The dimensions problems for the solve_ode_exp we will have to consider in another PR (will only do that after the Hessians).