Implicit Manifolds

[ChrisRackauckas]

In DiffEq we define implicit manifolds and project using a Newton method here:

https://github.com/JuliaDiffEq/DiffEqCallbacks.jl/blob/master/src/manifold.jl

We set this up with the generic nlsolve interface:

http://docs.juliadiffeq.org/latest/features/linear_nonlinear.html#Nonlinear-Solvers:-nlsolve-Specification-1

which defaults to using NLsolve.jl, this extra part can easily be dropped. Basically it's just that g is a residual function and you project to this manifold by calling a non-linear solver, not much more to it.

[antoine-levitt]

I see. The fact that the jacobian is singular and the solution non-unique is not a problem for NLSolve? (it does a SVD of the jacobian in the Newton method?)

Not sure if it's good to have this library depend on NLSolve though (because it would mean Optim takes on a dependency on NLSolve). Will there be a good solution to this problem post Pkg3?

[ChrisRackauckas]

I see. The fact that the jacobian is singular and the solution non-unique is not a problem for NLSolve? (it does a SVD of the jacobian in the Newton method?)

It seems to work fine.

Not sure if it's good to have this library depend on NLSolve though (because it would mean Optim takes on a dependency on NLSolve). Will there be a good solution to this problem post Pkg3?

Not necessarily. Requires.jl. What you can do is make it allow a generic nonlinear solver like I have it, and then just don't give a default nonlinear solver. That would be dependency free but require some docs for how to set it up with a nonlinear solver.

[antoine-levitt]

It seems to work fine.

Interesting. If the jacobian is not computed exactly (for instance, with finite differences) then that could create instabilities (small but nonzero eigenvalues in the jacobian, that then get inverted). I guess in practice the nullspace of the jacobian is aligned with the axes, so FD doesn't screw up and everything works fine. It could see it becoming more troublesome for more complicated manifolds though (e.g. the set of projectors, ie matrices satisfying X^2 = X).

Not necessarily. Requires.jl. What you can do is make it allow a generic nonlinear solver like I have it, and then just don't give a default nonlinear solver. That would be dependency free but require some docs for how to set it up with a nonlinear solver.

Then there's not much point in doing that in this library, it just adds a layer without actually saving any code: since the caller has to plug its solver into ImplicitManifold, it's simpler to just define a subtype of Manifold there and implement retract! on it.

One more question: in DiffEq you assume that the trajectory stays on the manifold (ie you just use retract and not project_tangent)? Do you have any interest in supporting "forcing" the ODE to stay in the manifold (and therefore changing the ODE x' = f(x) to x' = P_x f(x) where P_x is the projection on the tangent space at x)?

[ChrisRackauckas]

One more question: in DiffEq you assume that the trajectory stays on the manifold (ie you just use retract and not project_tangent)? Do you have any interest in supporting "forcing" the ODE to stay in the manifold (and therefore changing the ODE x' = f(x) to x' = P_x f(x) where P_x is the projection on the tangent space at x)?

The SSPRK methods specifically have step limiters for that:

http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws)-1

So if Manifolds.jl has projection operators that would make things easy for users.