SDE with matrix values, question on DifferentialEquations.jl

General Usage
1

I have the following SDE
du = f(u,p,t)dt + sum_i g_i(u,p,t) dW_i
where du, f, and g_i are all matrices, while W_i are all independent Wiener processes. I have no problem in setting up the ODE when g_i is zero, but I can’t figure out how to write the above problem in DifferentialEquations.jl in the general case, where all functions g_i have to be packed in a single function g.
Is g going to be a vector of matrices (each element being g_i), or a 3-index tensor such that g[:,:,i] = g_i[:,:] ? In any case it is not clear how can I set the noise_rate_prototype, or whether I should use a completely different approach.
Any ideas?
Thanks

2

Did you see this part of the tutorial which “translates” between the matrix and summation forms?

http://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Example-4:-Systems-of-SDEs-with-Non-Diagonal-Noise-1

3

Yes, indeed my “poor man” solution is based on vectorising everything: matrix u is transformed into a vector, so is f and g_i. In that case g becomes a matrix (whose columns are g_i), so I can use the example that you mentioned. However, in f and g I have to “devectorise” u to transform it back into a matrix, are my f and g are doing some non-trivial matrix operations that are not straightforward in the vectorised notation. So I think that my poor-man solution is not optimal, and I don’t see any reason why in julia it shouldn’t be possible to work directly with matrix-valued functions by playing with data types, rather than using my ugly trick of vectorising and devectorising the function each time.

4

Multiplication between a 3-tensor and 2-tensor isn’t well-defined though. Generalizing this definition in a way that doesn’t require a bunch of g functions passed in is hard, but having a bunch of functions rather than one in-place function would be very very inefficient. That’s the main issue that keeps it a noise matrix for now. But you can define it as any type you want and then define * against dW to have the correct interpretation. That would be a useful AbstractArray to work out and add to the library if you’re up for it.

Using reshape is cheap/free though since it’s just a view. There shouldn’t be a cost to this.

5

Thanks, I like your idea of defining * against a type and dW. I will see if I manage do to something, but otherwise I’ll keep my poor man solution.
Cheers