SDEs

[mateuszbaran]

See: https://arxiv.org/abs/2009.10113

[mateuszbaran]

Also efficient Brownian motion: https://arxiv.org/abs/2202.00959

[kellertuer]

Very interesting! But we would need probably (yet again) distributions on manifolds more thoroughly.

[AdrienAngeAndre]

I also recommend the geodesic Euler method of the paper https://arxiv.org/pdf/2312.14882, which is easy to implement in the package. I am developing higher order methods, but I need a few more months to write the paper and implement the methods in ManifoldDiffEq.

[mateuszbaran]

Very interesting! But we would need probably (yet again) distributions on manifolds more thoroughly.

Probably, I will have to investigate it. For now I wanted to keep an issue with interesting papers about SDEs on manifolds.

I also recommend the geodesic Euler method of the paper https://arxiv.org/pdf/2312.14882, which is easy to implement in the package. I am developing higher order methods, but I need a few more months to write the paper and implement the methods in ManifoldDiffEq.

Thanks! BTW, I am planning to work on getting some basic support for SDEs in ManifoldDiffEq later this year.

[mateuszbaran]

@AdrienAngeAndre I'm reading these papers about SDEs on manifolds and one thing I can't find is a generic and intrinsic definition of what an SDE on a manifold is. The paper you linked has a nice intrinsic definition but it seems to be limited to Langevin diffusion (we don't have an arbitrary drift term). Could you recommend a paper that deals with general SDEs? I mean, I can just replace $\nabla \phi$ by anything, especially in the algorithm, but I'm not sure if that's mathematically fine.

[AdrienAngeAndre]

@mateuszbaran You are indeed right. One of the main focus with SDEs on manifolds lies in Riemannian Langevin dynamics, with "additive noise" and a drift taking the form of a potential. The stochastic numerics litterature for general SDEs on manifolds is quite limited. I am writing a paper on the subject with new geometric methods and (hopefully) a general analysis, but I need a bit more time to finish it. I can send you a draft when the preprint is in a better shape and we can plan a meeting to discuss it if you want.

[mateuszbaran]

Thanks! I'm currently working on an idea of providing a mathematically consistent description of concepts appearing in state estimation in robotics with concrete numerical examples using JuliaManifolds. I'm not going to delve deep into SDEs but some people use them to describe system dynamics. I'd like to show how to solve it using geometric Euler–Maruyama, with an implication that it can be easily swapped for a higher order integrator when it is available.

I'm looking forward to reading your preprint when it's ready.

[mateuszbaran]

I've made the following script:

using Manifolds
using Markdown

using Distributions

using Makie
using GLMakie
using LinearAlgebra

@doc raw"""
    simple_stochastic_integrator(M::AbstractManifold, F, p0, h::Real, N:Int)

Simulate the stochastic process on manifold described by

``math
d p(t) = F(p(t)) dt + dB^M(t) 
``
where ``p(t) \in \mathcal{M}``, ``F\colon \mathcal{M} \to T\mathcal{M}`` is the drift term
and ``B^M`` is the noise term.

Discretization step `h`, `N` time steps

See: [Bharath:2023](@cite).
"""
function simple_stochastic_integrator(
    M::AbstractManifold,
    F,
    p0,
    h::Real,
    N::Int;
    noise=MvNormal(zeros(manifold_dimension(M)), 1.0),
    retr::AbstractRetractionMethod=ExponentialRetraction(),
    B::AbstractBasis{ℝ}=DefaultOrthonormalBasis(),
)
    ps = [p0]
    p_i = p0
    for i in 1:N
        ξ = rand(noise)
        # in normal coordinates metric corrections vanish
        X_i = h * F(M, p_i) + sqrt(h) * get_vector(M, p_i, ξ, B)
        p_i = retract(M, p_i, X_i, retr)
        push!(ps, p_i)
    end
    return ps
end

function test()
    M = Manifolds.Sphere(2)
    h = 0.01
    N = 1000
    p0 = [0.0, 0.0, 1.0]
    no_drift(M, p_i) = zero_vector(M, p_i)
    whirly_drift(M, p_i) = cross(p_i, [1.0, 0.0, 0.0])
    noise = MvNormal(zeros(manifold_dimension(M)), 0.3)
    ps = simple_stochastic_integrator(M, whirly_drift, p0, h, N; noise=noise)
    lines([p[1] for p in ps], [p[2] for p in ps], [p[3] for p in ps])
end

which seems to implement the manifold SDE integration:

Does that looks OK?

I'm going to generalize it a bit but I'm not entirely sure how to handle changes to noise term with time. I'd like to handle at least as general noise terms as considered here: https://arxiv.org/abs/2108.09387 (though I know they do things a bit differently).

[AdrienAngeAndre]

The code is good to my opinion, though I would not say that this is a "simple_stochastic_integrator". ;)
Just to check: the line "get_vector(M, p_i, ξ, B)" returns an element of the tangent space at p_i whose coordinates in an orthonormal basis each are Gaussians, right? If yes, the code is good I think.
The extension for multiplicative noise is not completely clear to me yet.

[mateuszbaran]

Thanks! I didn't have a better naming idea.

Just to check: the line "get_vector(M, p_i, ξ, B)" returns an element of the tangent space at p_i whose coordinates in an orthonormal basis each are Gaussians, right? If yes, the code is good I think.

Yes, that's what it does.

The extension for multiplicative noise is not completely clear to me yet.

I think additive noise is enough for me but I'm currently considering having an arbitrary covariance matrix for the Gaussian distribution and transporting it along at each step. Lie group Kalman filter people seem to just implicitly use a flat Cartan-Schouten connection for transport but that choice seems somewhat arbitrary to me.