SDE driven by a NoiseTransport noise

[Damiano]

Hi everybody,

first of all let me say that I still consider myself new to Julia (and even programming in general), so everything I say and write is likely to contain very naive mistakes. Also, if there is anything I should improve regarding how to post/ask questions, please let me know.

I will first write my final goal, just in case I am approaching it from a completely wrong/super inefficient angle and you know a better way to do it. After that I will write my actual question, which is (I think) much simpler and concerns basic usage of the DifferentialEquations.jl package.

Goal: simulate a system of coupled SDEs, let’s say 2N, driven by N complex Brownian motions. The first N equations are driven by the N complex Brownian motions, the second N equations are driven by the conjugate of the N Brownian motions.

Strategy: implement a custom noise using DiffEqNoiseProcess.jl. In particular, I want to construct a NoiseTransport process and then use it to drive my system of SDEs.

Problem: I don’t seem to be able to follow the instructions of the documentation. I link below the sections I read the most, but of course I tried to at least skim through all what seemed to be relevant.
DiffEqNoiseProcess.NoiseTransport
Stochastic Differential Equations
It seems that when it comes to using a NoiseTransport noise, I am not able to implement even a simple Brownian motion. Of course in this case I could just set up my problem as explained by the general tutorial of DifferentialEquations.jl (without the need of DiffEqNoiseProcess.jl), since Brownian motion is already the default noise of SDEProblem(…). The point here is that I am trying to implement the easiest nontrivial NoiseTranport noise I can think of.

Minimal (not working) example:

using DifferentialEquations

function f(u,p,t,rv)
    rv
end

t0 = 0.0
T = 10
u0 = 0.0 

W = NoiseTransport(t0, f, randn)

function drift(u,p,t)
    0
end

function diffusion(u,p,t)
    1
end

my_prob = SDEProblem(drift, diffusion, u0, (t0,T), noise = W)
my_sol = solve(my_prob)

but the values of the solution my_sol are simply zero for all time points.

I hope some of you can point out what I am making wrong.
I am sure I am missing something very obvious, so I apologize in advance.

Thanks,

Damiano

[Damiano]

I’m replying to my same question, with the idea that maybe there are people who know how to answer it, but did not have the time to read it before new questions made it go down the list. If this is inappropriate behaviour (or if it okay but I did not wait enough), please tell me: I will apologize and not do it again!

[Damiano]

Let me push it one more time. Just to be absolutely clear, I think my problem is either with the understanding of how to use NoiseTransport (and the other abstract noise processes) or with the goal that NoiseTransport tries to achieve. Over the last few days I tried to go through the documentation once more. I also opened two very small issues on GitHub regarding the documentation of DiffEqNoiseProcess (by the way it was my first real action ever on GitHub, so I hope I did everything right there). I can try to share my understanding (or lack thereof) of what I’m doing. I initially refrained to do so because very often such things can be distracting for people who answer: I don’t think one wants to go through reasonings and examples of somebody who is talking nonsense because they don’t understand anything of what they’re doing. But I don’t see what else I could do (a part from reading and trying to understand the source code, which frankly seems to be a bit beyond of my skills right now).

I will still write some thoughts, I guess with the only goal of showing that I did try to think about this on my own.
Maybe NoiseTransport was never supposed to be able to create a Brownian motion in the first place. Indeed, the expression W(t) = f(u, p, t, RV) written in DiffEqNoiseProcess.NoiseTransport seems to suggest that NoiseTransport can only create a noise made as following. First of all, I sample RV at time 0. Then RV will then be sampled again. I would then interpret the function f I defined above as a constant. The noise that multiplies the diffusion would then be the time derivative of a constant, which is 0. Thus it would be correct that the solution of the code I wrote in my first post is identically 0. This theory is confirmed by the fact that, if I replace f in the code of my first post by f(u,p,t,rv) = rv * t, then the solution I get is just a line, starting at 0, whose slope changes each time I run the whole code. However, if this is indeed the intended use, I don’t why NoiseTransport exists at all. Since if I only sample RV at time 0, the evolution of the dynamics is actually deterministic. Indeed, I could write something like

du = ff(u,p,t)dt + g(u,p,t)dW 
     = ff(u,p,t)dt+ g(u,p,t) df(u,p,t,RV) 
     = ff(u,p,t)dt+ g(u,p,t) f'(u,p,t,RV)dt
     = (ff(u,p,t) + g(u,p,t) f'(u,p,t,RV))dt 
     = fff(u,[p,RV],t)dt

so that after an initial sample of RV I could just pass it as a parameter in an usual ODE solver. Indeed, the following code (added at the end of the code of my previous post, for the sake of having a minimal working example)

function fff(u,p,t)
    p
end
my_prob_2 = ODEProblem(fff,u0,(t0,T),randn())
my_sol_2 = solve(my_prob_2)

seems to produce exactly the same output as the code of my original post after replacing f(u,p,t,rv)=rv by f(u,p,t,rv)=rv*t.

I think I will stop it here because I feel like I already went quite long. I am sure I am missing something obvious, also in view of the sentence (italics added by me) “Lastly, the NoiseFunction allows you to use any function of time as the noise process, while NoiseTransport lets you define a random process as the transport of a random variable or a random vector by a time-dependent function. Together, these functionalities allow you to define any colored noise process and use it efficiently and accurately in your simulations.” that can found on the documentation pages I linked above. Of course the kind of noise I described in the last few lines if very far from “any colored noise process”.

As an extreme measure, I will allow myself to invoke @ChrisRackauckas (this is the only name I can link with this package that also actively answers here, but I am sure there are many more: sorry for calling out you exactly just out of my ignorance).

Thank you a lot to everybody that read through this and/or can offer any light on the matter.

[ChrisRackauckas]

@Damiano

Goal: simulate a system of coupled SDEs, let’s say 2N, driven by N complex Brownian motions. The first N equations are driven by the N complex Brownian motions, the second N equations are driven by the conjugate of the N Brownian motions.

This should just happen if your u0 is a complex number. Did you try just using a complex number? I don’t know what the rest of your discussion is about but if this is what you’re trying to do then:

using DifferentialEquations
function lorenz(du, u, p, t)
    du[1] = 10.0(u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end

function σ_lorenz(du, u, p, t)
    du[1] = 3.0
    du[2] = 3.0
    du[3] = 3.0
end

prob_sde_lorenz = SDEProblem(lorenz, σ_lorenz, [1.0 + 0im, 0.0 + 0im, 0.0 + 0im], (0.0, 10.0))
sol = solve(prob_sde_lorenz)

[Damiano]

First of all thank you a lot for answering! Also thanks for the documentation fixes, which are already in place.

I had tried to use a complex number as initial condition, as implied by the documentation of this funciton DiffEqNoiseProcess.RealWienerProcess. However this does not seem to do what I want. Indeed, consider the following simplification of the code you suggested:

using DifferentialEquations
function lorenz(du, u, p, t)
    du[1] = 0
    du[2] = 0
    du[3] = 0
end

function σ_lorenz(du, u, p, t)
    du[1] = 1
    du[2] = 1
    du[3] = 1
end

prob_sde_lorenz = SDEProblem(lorenz, σ_lorenz, [0.0 + 0im, 0.0 + 0im, 0.0 + 0im], (0.0, 10.0))
sol = solve(prob_sde_lorenz)

Based on what I can deduce from the piece of documentation I linked above, this code would simulate a system of three independent complex brownian motions. In other words, it would simulate the process

[W_1 + im * W_2, W_3 + im * W_4, W_5 + im * W_6]

where W_i, for each i = 1, ... , 6 is an independent real brownian motion. What I am trying to do is to implement some combination of diffusion/noise such that the solution of SDEProblem, in the case of drift constantly equal to 0, simulates the process

[W_1 + im * W_2, W_3 + im * W_4, W_1 - im * W_2, W_3 - im * W_4]

(as I wrote in piece of text you quoted, I always have an even number of equations and the noise for the second half is the conjugate of the noise of the first half).
As a confirmation that the real and imaginary parts of the brownian motions of the solution you proposed are not related in any way, the execution in the REPL of

sol.u

gives

109-element Vector{Vector{ComplexF64}}:
 [0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im]
 [0.0017131513590476093 - 0.0006250927629785429im, 0.0010197924502198723 + 0.0012471454943345324im, 0.000653257322638947 - 0.00017127908183037224im]
 [0.002207526263179617 + 0.0036781088418133966im, 0.002768519434565341 + 0.0026836406975075384im, -0.0017787045115214978 - 0.0020589970592874867im]
 [0.0027493475256336816 + 0.006009030839977256im, 0.002342665904362015 + 0.0024599074181595386im, -0.0005205704035833108 - 0.0003792961799878405im]
 [0.002647937671740059 + 0.0029110705490700555im, 0.0022402047691080984 + 0.003099232176846538im, -0.0026815161100179866 - 0.0028733868444415565im]
 [0.0024633868435726338 + 0.0033989325723413825im, 0.00025733359934868895 + 0.0019914774608010854im, -0.00019074332595502615 - 0.0023377832144141518im]
 ⋮
 [-1.8339887287606143 - 0.007526533831593737im, 0.2314740926962144 + 2.9764883882672994im, -0.47139960064430053 + 1.5955407929057597im]
 [-1.2844429937970574 + 0.9974373672210978im, 1.3281472728191668 + 1.9679572471149926im, 0.14055302484188031 + 1.7411722109444436im]
 [-0.5167516435167138 + 1.3958015464686786im, 0.8937907782727026 + 1.4142167739454101im, 0.04944291205192954 + 2.3818698278487895im]
 [-0.5647013934554475 + 0.31442037382691046im, 0.9941871310384436 + 1.9051200769480703im, 0.5958414463749292 + 1.6565415356954278im]
 [-0.405325999187691 - 0.7213116491948023im, 0.8294968381151699 + 2.99538005174791im, 0.16415631227578742 + 1.600492506435119im]
 [0.2620229251008528 - 0.22678144794864222im, 1.327186728685124 + 3.228625345114374im, 0.04304087987500906 + 1.6367654340488544im]

and there are no equal numbers.

I actually managed to get the desired result by implementing a non diagonal noise, thanks to Example 4: Systems of SDEs with Non-Diagonal Noise. However the diffusion function I wrote seems (to me) unnecessarily complicated. It does quite a number of allocations, and as a result it is much slower than the implementation you suggested (which does not seem to simulate my model), while being, in principle, simpler (my model actually has less randomness than the one of the implementation you suggested). This is what brought me to the DiffEqNoiseProcess.jl package.

The idea of being missing something very simple kept me from posting my working, slow and complicated implementation: I did not think it would have encouraged people to answer. But if you want to have a look at it, I am happy to post it.

[ChrisRackauckas]

@Damiano

I actually managed to get the desired result by implementing a non diagonal noise, thanks to Example 4: Systems of SDEs with Non-Diagonal Noise.

Yes, I was going to say that this would be the right approach. You can write the diffusion to be non-allocating this way and it’ll be the simplest representation that the solver can specialize on. If you start going down the DiffEqNoiseProcess.jl path then you might have issues with scaling the solver, I’d stick to the simple non-diagonal noise representation.

[Damiano]

I see. My implementation using non diagonal noise, at the moment, is everything but non allocating (I mean, it has quite a number of allocations according to @btime). Let me work on it a bit more and/or read again some of the resources on memory allocation you already mentioned multiple times in other posts of yours. If I manage to improve it to the point where it is good enough for me, I’ll mark your reply as answer. Otherwise I will post my best implementation (cleared of the other details of the model I’m trying to implement) and ask again for your suggestion.
Thanks again.