Complex ODE and Manifold projection

[benneti]

As it easier to see a simple example of the problem:

f(u, p, t) = - im * u; u0 = 1.0+0im
prob = ODEProblem(f,u0,(0.0,100.0))
function g(resid, u, p, t)
       resid = abs2(u) - 1
end
sol = solve(prob,Vern7(),callback=cb)

Rises the following error
ERROR: MethodError: no method matching nlsolve(::NLSolversBase.OnceDifferentiable{Array{Complex{Float64},1},Array{Complex{Float64},2},Array{Complex{Float64},1}}, ::Complex{Float64}).
Is there anyway to get around this, for this simple example I see, that one could equivalently do the example of the documentation and just use the real problem corresponding to it but for more complicated problems converting them to real problems is a bit cumbersome.

As an independent side note I want to thank you for this library and the videos you put on youtube, as there i stumbled over julia two weeks ago and already enjoy it very much!

[ChrisRackauckas]

@pkofod is NLsolve supposed to handle Complex as defined with these types? I thought it did?

[pkofod]

Don't think we ever have supported complex values only complex input. Are there no problems when calculating the Jacobian? Is the Complex-Complex derivative well-defined? @antoine-levitt

[antoine-levitt]

For methods with derivative no, in general there's no more structure to the jacobian than a real 2N by 2N matrix (you can represent it as Wirtinger derivatives, ie two complex NxN matrices, but I don't think that helps with matrices inversions). Broyden or Anderson do generalize to complex mappings, but are not equivalent to the flattened 2n real version. Not sure if there's much point in doing anything else than wrapping as real 2n arrays (although supporting the native Anderson and Broyden might be nice and would be quite easy)

[ChrisRackauckas]

The Jacobians work just fine for nonlinear solving. We use it internally in DiffEq complex solving.

[antoine-levitt]

Does it also work for non holomorphic functions? Try solving z->conj(z) - (1+i) or something similar.

edit: that might actually be too simple. Try A*conj(z) - b for a random complex matrix A and random complex vector b.

[benneti]

For this problem it would be sufficient to a real funciton with complex arguments.
Is there some way to tell the callback that I do not want resid to be like u but be a real number (Float64)?

[antoine-levitt]

OK, I did not look at the OP carefully, this is a totally different problem: you want to solve a differential equation (in this case, a very simplified Schrödinger equation) and ensure that you stay on the sphere. In that case, what you want is not the generic nlsolve, but a custom callback that renormalizes the solution (u /= norm(u))

@ChrisRackauckas :

1 is the approach above sufficient to preserve the order of the integrator?

2 that might be a good time to revisit our earlier discussion on integrating the manifolds from Optim. Can you try to get the OP's example working using

using Optim
M = Optim.Sphere()
Optim.retract!(M, u)

to perform the normalization? Maybe as a separate method ManifoldProjection(::Optim.Manifold)? Then if @pkofod agrees I can split Manifolds to a separate package on which Optim depends. Ideally, we could even have an ImplicitManifold <: Manifold that uses NLSolve to perform the retraction, simplifying the interface and allowing Optim users to get that also, but I'm not sure how to go about that (we probably don't want Optim to depend on NLSolve)

In general, the user could specify an ODE that does not stay explicitly on the manifold (ie a vector field f that is not in the tangent space). In that case one would have to use Optim.project_tangent as well.

[ChrisRackauckas]

1 is the approach above sufficient to preserve the order of the integrator?

Yes. If it's applied every once in awhile to the solution then it's an order-preserving operation. It breaks the interpolant but that's it.

2 that might be a good time to revisit our earlier discussion on integrating the manifolds from Optim.

Yes, this is definitely a case where a non-general method should be used, dispatching on the type of manifold to do a much simpler operation than a nonlinear solve.

[antoine-levitt]

Yes. If it's applied every once in awhile to the solution then it's an order-preserving operation

Yes, I guess even if it's applied every time-step it works too, because the deviation from the manifold is of high order anyway. In practice do you apply it every step or do you have a tolerance on the residual?

[ChrisRackauckas]

Yes, I guess even if it's applied every time-step it works too, because the deviation from the manifold is of high order anyway. In practice do you apply it every step or do you have a tolerance on the residual?

Yes, that's exactly the case. It's a triangle inequality argument to show that the perturbation is less than the order of the algorithm (if done every step) and therefore the order is kept. In practice it's only done every once in awhile when above some tolerance. We just call NLsolve each time though, and it auto-exits if above tolerance.

[pkofod]

OK, I did not look at the OP carefully, this is a totally different problem: you want to solve a differential equation (in this case, a very simplified Schrödinger equation) and ensure that you stay on the sphere. In that case, what you want is not the generic nlsolve, but a custom callback that renormalizes the solution (u /= norm(u))

wow I didn't either! It's even in the title :)

Yes, that's exactly the case. It's a triangle inequality argument to show that the perturbation is less than the order of the algorithm (if done every step) and therefore the order is kept. In practice it's only done every once in awhile when above some tolerance. We just call NLsolve each time though, and it auto-exits if above tolerance.

Yes, I remember this discussion. This is why you were interested in having cache variable structs, right ? So these calls are essentially free to the extent that you're only really calculating the norm, not taking any steps.

Maybe as a separate method ManifoldProjection(::Optim.Manifold)? Then if @pkofod agrees I can split Manifolds to a separate package on which Optim depends. Ideally, we could even have an ImplicitManifold <: Manifold that uses NLSolve to perform the retraction, simplifying the interface and allowing Optim users to get that also, but I'm not sure how to go about that (we probably don't want Optim to depend on NLSolve)

This seems like a natural/Julian thing to do. We just need to make sure that the projections can still be applied where needed. We might have to change some things in Optim. But if it's just a matter of moving things out, then it shouldn't be too much of a problem.

[pkofod]

@antoine-levitt I've created https://github.com/JuliaNLSolvers/ManifoldProjections.jl/

[antoine-levitt]

Great! Can you give me push rights?

[pkofod]

Done

[benneti]

How does this projection work currently, as my example in the OP was very simplified would it be possible to have a state vector in C^n and only project a part of it like x[k:n] on a sphere of dimension n-k?

And is there any way, I as a julia newb can help?

[antoine-levitt]

Sure, you just have to define a new manifold type that does this. Look at how Sphere is defined in the above repository to see how this is done.

If you want to contribute, a great way to start is to find out how to solve your problem using a callback and ManifoldProjections, and to make a pull request to add that as an example to the DiffEq docs (maybe just for the sphere as it's more widely applicable).

[pkofod]

Would it be too specific to allow this encoding in ManifoldProjections.jl? So you could define your sphere with an index specification (defaulting to Colon [or just a fallback that doesn't index])? Say an "Indexed Sphere Projection" would just read

using Parameters, LinearAlgebra
@with_kw struct Sphere{T} <: Manifold
    indices::T = nothing
end

retract!(S::Sphere{<:Nothing}, x) = normalize!(x)
project_tangent!(S::Sphere{<:Nothing},g,x) = (g .-= real(dot(x,g)).*x)
@views retract!(S::Sphere, x) = normalize!(x[S.indices])
@views project_tangent!(S::Sphere,g,x) = (g[S.indices] .-= real(dot(x[S.indices],g[S.indices])).*x[S.indices])

[antoine-levitt]

It really takes three lines to define such a custom manifold so I'd rather we just document how it's done rather than add complexity to poor old Sphere. Also note this can be achieved with the Product manifold (although its interface is currently pretty bad)

[benneti]

Do you think it would be helpful to have a sphere with an radius, as I think the current implementation could lead to unnecessary numerical problem.
For example if one has a sphere in R^n with a big radius r and a point close to (r, 0, 0, 0, ...) using the current implementation one would need to normalize it and then scale it up again.
I think the numerically more stable approach would be to implement the sphere like this:

struct Sphere <: Manifold
    r
end
retract!(S::Sphere, x) = (x .= x .* (S.r / norm(x)))

where the numbers we divide should have similar magnitude.

[antoine-levitt]

Due to the way floating point numbers are stored, multiplication is not numerically unstable, ie very little precision is lost by dividing and multiplying by a given number (unless you get to extremely large or small numbers, which you probably shouldn't do anyway). It's addition that's problematic. That said you can always define a new Sphere type, but (as with the other type) I think it's simple enough to do that it doesn't need to be put in the library.

[benneti]

Ok, so the best solution to project on a sphere with a radius would be the code above (but replacing Sphere with SphereR)?

[antoine-levitt]

Yes! You also typically want r to be statically typed:

struct Sphere{T} <: Manifold where {T <: Real}
    r::T
end

[benneti]

Thank you very much, then I'll try to create a Callback using this.

[benneti]

This function creates a callback as does ManifoldProjection but using a https://github.com/JuliaNLSolvers/ManifoldProjections.jl/ .
@ChrisRackauckas do you think this should be in some way added to manifolds.jl?

function ManifoldRetraction(M::Mfd; save=true) where {Mfd <: Manifold}
    affect!(integrator) = retract!(M, integrator.u)
    condition = (u, t, integrator) -> true
    save_positions = (false,save)
    DiscreteCallback(condition, affect!,
                     save_positions=save_positions)
end

@antoine-levitt To the sphere projection without considering the numerics, wouldn't it still be useful to change the sphere manifold to also be able to change the radius. I think the complexity added to the library is marginally is big -- still you are right that writing a own SphereR class is not really hard.
Should we move this discussion to the issue tracker of ManifoldProjections.jl or are you absolutely sure that it is unnecessary clutter?

[pkofod]

Should we move this discussion to the issue tracker of ManifoldProjections.jl or are you absolutely sure that it is unnecessary clutter?

Yes I think this is a discussion worth having over there.

[ChrisRackauckas]

@ChrisRackauckas do you think this should be in some way added to manifolds.jl?

Yes, I think it would be nice to have a pre-built callback where you just add a manifold via ManifoldProjections.jl and it projects to that manifold. I would accept that in DiffEqCallbacks and the DiffEqDocs since that seems like a pretty good addition to have.

[benneti]

Any prefereces on function naming?

[ChrisRackauckas]

I think it would be nice if ManifoldProjection worked by multiple dispatch and giving it this kind of manifold is a valid alternative input. That might be wonky though, so I'm not sure...

[paolo-mgi]

Please bear with me I am not an advaced programmer. It's not clear to me if there is a solution for implementing a manifoldprojection callback when the differential equation involves complex numbers as in the original OP.
From the discussion it seems to me that it is feasible but I do not understand exactly how.
Many thanks!

[ChrisRackauckas]

This was solved with the new infrastructure. ManifoldProjection now uses NonlinearSolve.jl which handles complex. If you have any issues, please open an MWE in DiffEqCallbacks.jl