# SDE with matrix values, question on DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/sde-with-matrix-values-question-on-differentialequations-jl/10412>\
**Category:** General Usage\
**Tags:** sde\
**Created:** [April 18, 2018, 1:03pm UTC](https://discourse.julialang.org/t/sde-with-matrix-values-question-on-differentialequations-jl/10412 "2018-04-18T13:03:56Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![benkj](https://avatars.discourse-cdn.com/v4/letter/b/848f3c/32.png) [@benkj](https://discourse.julialang.org/u/benkj)\
**Post date:** [April 18, 2018, 1:03pm UTC](https://discourse.julialang.org/t/sde-with-matrix-values-question-on-differentialequations-jl/10412/1 "2018-04-18T13:03:56Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 18, 2018, 4:01pm UTC](https://discourse.julialang.org/t/sde-with-matrix-values-question-on-differentialequations-jl/10412/2 "2018-04-18T16:01:42Z")

</div>

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](http://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Example-4:-Systems-of-SDEs-with-Non-Diagonal-Noise-1)

---

<div class="post-metadata">

**Author:** ![benkj](https://avatars.discourse-cdn.com/v4/letter/b/848f3c/32.png) [@benkj](https://discourse.julialang.org/u/benkj)\
**Post date:** [April 18, 2018, 5:06pm UTC](https://discourse.julialang.org/t/sde-with-matrix-values-question-on-differentialequations-jl/10412/3 "2018-04-18T17:06:22Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 18, 2018, 5:17pm UTC](https://discourse.julialang.org/t/sde-with-matrix-values-question-on-differentialequations-jl/10412/4 "2018-04-18T17:17:02Z")

</div>

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.

> [@benkj](#):
>
> 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

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

---

<div class="post-metadata">

**Author:** ![benkj](https://avatars.discourse-cdn.com/v4/letter/b/848f3c/32.png) [@benkj](https://discourse.julialang.org/u/benkj)\
**Post date:** [April 18, 2018, 5:55pm UTC](https://discourse.julialang.org/t/sde-with-matrix-values-question-on-differentialequations-jl/10412/5 "2018-04-18T17:55:13Z")

</div>

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