projected SDE

[AshtonSBradley]

following up on #145 discussion on callbacks, I am trying to track down any issues that might arise in formulating the SDE in this way. There is a potential snag:

when doing the integration it is essential that the RHS of the SDE be projected. That is, I think there will be a fundamental problem in using a discrete callback to:

In the algorithm you want to do, you just need to project after doing the update. Because we know that this converges (under appropriate assumptions), I think it's fine to explicitly do that projection to the SDE. I think the best way to implement this algorithm as an SDE would be in two steps. The first step is to define f and g as though there is no projection. Then use a DiscreteCallback. DiscreteCallbacks can be set to happen after every update ...

because that would amount to solving a different equation of motion, i.e. the RHS has to be projected because it turns up in the sub-steps of the integration routine, even with a numerical tilmestep dt.

The problem is that, even in the drift terms, the notation (using the trapping potential and kinetic terms here for example)

d\psi(x,t) = P[-I(p^2/2m + V(x))\psi(x,t)dt]

is a shorthand for a set of equations for coefficients in a truncated basis expansion:

\psi(x)=\sum_n c_n(t) \phi_n(x)

of the form

dc_n(t) = \int \phi_m^*(x)[-i(p^2/2m + V(x))\psi(x,t)dt]

There is a projected functional calculus that goes along with this kind of equation, relating the general diffusive SDE evolution to a functional Fokker-Planck equation in this projected functional space.

So the equations are most certainly SDEs, but it is essential that the projector always appear. To first approximation we could think of them as a discretised field SDE. These objects are known to converge, in a certain sense, as the grids get fine enough, to a continuous field SDE. For the projected SDE, there is no need to make the grids finer: the projected SDE is a well defined SDE in a projected functional space. Anyway, the point is that the RHS must be projected at every evaluation. Even for the Hamiltonian case (turning off noise and damping), if the projector is there, then this is a purely Hamiltonian evolution (even with nonlinear terms, when evaluated correctly) that also conserves particle norm. If it is not there at every evaluation, then these properties are lost.

I assume callbacks do not give access to this sort of low-level meddling with the integrated variables inside dt?

A further point worth making - in the above equations, one would typically have an additional nonlinear term (s-wave scattering V(x) -> V(x) + g|\psi(x)|^2), but the single-particle basis is usually chosen to diagonalise the single-particle terms, so that this equation without nonlinearity just becomes diagonal in coefficient space:

dc_n(t) = -ie_n c_n dt

for single-particle eigenstate energies e_n.

I am a little worried about the direction of giving up this nice view of (a) getting single particle terms "for free" with the projection, and (b) having the nice properties of the projected equation.

So my question is this: say I have some generalised SDE that is expressed in a non-standard basis, but can be shown in another basis to be an SDE. How do I pass this information to the solver? Could I have a condition that the determinant of the diffusion matrix is evaluated to check that the process is a true diffusion, either in the coefficient basis, or in some basis of the user's choosing? On second thoughts, this would not work because the properties of the diffusion matrix are related to mapping to Fokker-Planck equation, not whether the SDE is in standard form.

This seems like a big ask, maybe you already know of a solution, or maybe there is a research question in there the needs to be answered first?

(btw - super cool to see you implemented the generalised noise process, I am looking forward to trying it out!)

[ChrisRackauckas]

I need to work on my advancement, so this is not going into DiffEq 2.0 which should start its release process this weekend. This is not going to make it, so it will take a bit.

This seems like a big ask, maybe you already know of a solution, or maybe there is a research question in there the needs to be answered first?

Can you clarify what exactly you need? Not the model, please abstract the computation away from the model and the physics. It'll help find out the true underlying form. If you were to write an Euler-Maruyama update on this equation, what would it look like? Would you

utilde = Pinv(uprev)
utmp = f(t,utilde)*dt + g(t,utilde)*dW
u = P(utmp)

? Or is

Anyway, the point is that the RHS must be projected at every evaluation.

Are you actually talking about a SDAE? A projection method like:

https://reference.wolfram.com/language/tutorial/NDSolveProjection.html

? I plan on doing this stuff with

http://docs.juliadiffeq.org/latest/types/refined_ode_types.html#Mathematical-Specification-of-a-Constrained-ODE-Problem-1

but it will take awhile to get it over to SDEs.

I assume callbacks do not give access to this sort of low-level meddling with the integrated variables inside dt?

It actually does... but only after the full step. It has access to all of the intermediate variables and the cache, but the full step completes before the callback is called.

maybe there is a research question in there the needs to be answered first?

Maybe, but I am unsure what the question is.

[AshtonSBradley]

good point, I will work on a more concrete description of the problem to be solved.

[ChrisRackauckas]

If it is not there at every evaluation, then these properties are lost.

Why does it have to be at every evaluation instead of at the end of the timestep? The latter seems more numerically stable (and can be shown to be for ODEs), and there's no loss of order by only projecting once. I will be modifying all of the standard problems to allow for a function manifold(u) which is a definition of a manifold that the solution should live on. From what I've read simply projecting to the manifold at the end of the step is sufficient for structure-preserving properties and conservation of order, and it's actually better for a few reasons too.

Would this kind of projection solve your problem? Then you can define the SDE, and the manifold it should live on. The manifold can just be some constraint equations that must be satisfied, and it will project to a nearby solution on the manifold each step. I don't know of theory for this on SDEs, but for ODEs this is sound.

[ChrisRackauckas]

End of step projections seem justified for SDEs:

https://arxiv.org/pdf/1601.04157.pdf

[AshtonSBradley]

Interesting! Our current approach is to work in a basis that allows for very accurate projection (i.e. at machine precision) to be built into the ODE evolution (probably excessive accuracy in the SDE context). Maybe it would be ok to keep using the basis approach but implement projection as a callback after all. This would give a large speed improvement.

How would manifold(u) implement the projection numerically?
Does the need to truncate noises (page 6, https://arxiv.org/pdf/1601.04157.pdf) have any implications for the adaptive routines?

[ChrisRackauckas]

Does the need to truncate noises (page 6, https://arxiv.org/pdf/1601.04157.pdf) have any implications for the adaptive routines?

If you interpret the truncation as part of the method and not of the underlying Brownian path, I believe the convergence results due to Gaines and Lyons still holds. This means the NoiseProcess should be un-truncated, as the truncation there may bias the Brownian Bridge statistics. However, all uses inside of the actual methods could just check if a manifold is given and if so apply truncation.

Milstein actually has another paper which shows that this kind of truncation also allows for fully implicit methods (not just implicit in the drift term, but also the diffusion term). But in general it seems to be more of a theoretical tool. Truncation methods were developed to get around the non-boundedness of the moments for the inverted Normal distribution by truncating and thus forcing it to be bounded. However, there's a lot of very recent work to show that convergence still occurs properly without the truncation:

http://www.sciencedirect.com/science/article/pii/S0377042715003210
http://www.sciencedirect.com/science/article/pii/S0096300313008722
https://arxiv.org/pdf/1203.5809.pdf
https://arxiv.org/pdf/1609.08101.pdf
https://arxiv.org/pdf/1010.3756.pdf

there's instead a tamed Euler scheme which does a different scaling. Generally, what all of these are showing is that you need a certain relation between the truncation scheme and the Lipschitz coefficients (or whatever other bound in the assumptions) for the convergence result to hold s.t. they go to infinity in the correct way. This means that the truncation bound can be as high as you want, as long as you take it to infinity as you decrease the stepsize to zero in a specific way... which is to say it's a proof technique moreso than an actual numerical technique. But I can easily introduce an option to allow for truncating.

How would manifold(u) implement the projection numerically?

Rootfinding. Since it's O(h^p) away, it should find a piece of the manifold that is very close. The paper just uses Newton iterations, but starting close enough Newton = Trust-Region, so it seems that it truly doesn't matter what rootfinding technique you use here as long as its some local form. So I'll start by making it a call to NLsolve, and generalize it so the user can pass in whatever nonlinear solver method they want.

[ChrisRackauckas]

I created a callback for DiffEqCallbacks.jl that implements the manifold projection:

using OrdinaryDiffEq, Base.Test, DiffEqBase, DiffEqCallbacks

u0 = ones(2)
f = function (t,u,du)
  du[1] = u[2]
  du[2] = -u[1]
end
prob = ODEProblem(f,u0,(0.0,100.0))

function g(u,resid)
  resid[1] = u[2]^2 + u[1]^2 - 2
  resid[2] = 0
end

cb = ManifoldProjection(g)
sol = solve(prob,Vern7())
@test !(sol[end][1]^2 + sol[end][2]^2 ≈ 2) # Does not conserve energy
sol = solve(prob,Vern7(),callback=cb)
@test sol[end][1]^2 + sol[end][2]^2 ≈ 2 # Conserves energy

using Plots; plot(sol,vars=(1,2))

Documenting and releasing soon. Even if this doesn't solve your problem, it'll be useful elsewhere!

[ChrisRackauckas]

Projection docs: http://docs.juliadiffeq.org/latest/features/callback_library.html#Manifold-Conservation-and-Projection-1 . For implicitly defined g(u)=0 manifolds and conservation laws.

[ChrisRackauckas]

As far as I can tell, the callbacks have solved this.