SDE with matrix values, question on DifferentialEquations.jl



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?


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


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.


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.


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.