noise_rate_prototype for OOP problems

[rvignolo]

Hi,

for me, it is not clear what is the best way to define a SDEProblem when there is non-diagonal noise.

First, a common header for all:

using BenchmarkTools
using StochasticDiffEq
using StaticArrays
using SparseArrays

f(u, p, t) = u
function g(u, p, t)
    return @SMatrix [0.30u[1] 0 0 0.12u[2]
                     0        0 0 1.80u[2]]
end

f(du, u, p, t) = du .= u # copyto!(du, u)
function g(du, u, p, t) 
    du[1,1] = 0.30u[1]
    du[1,4] = 0.12u[2]
    du[2,4] = 1.80u[2]
end

oop_u0 = @SVector ones(2)
iip_u0 = ones(2)
dt = 1//2^(4)
tspan = (0., 1.)

Out of Place version

My first idea was using a noise_rate_prototype as a StaticArray as well:

oop_gprototype = SMatrix{2,4,eltype(oop_u0)}([1 0 0 1
                                              0 0 0 1])

oop_prob = SDEProblem{false}(f, g, oop_u0, tspan, noise_rate_prototype = oop_gprototype)
oop_sol = solve(oop_prob, EM(), dt = dt)
@btime solve($oop_prob, EM(), dt = $dt) # 2.668 μs (31 allocations: 3.05 KiB)

Adding sparse capabilities seems to worsen the situation:

oop_gprototype = sparse(oop_gprototype)

oop_prob = SDEProblem{false}(f, g, oop_u0, tspan, noise_rate_prototype = oop_gprototype)
oop_sol = solve(oop_prob, EM(), dt = dt)
@btime solve($oop_prob, EM(), dt = $dt) # 13.295 μs (262 allocations: 25.86 KiB)

This could be related to the fact that the underlying data for the SparseMatrixCSC1 seems to be an Array, not a StaticArray:

julia> oop_gprototype.nzval
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

Do you think this an issue that we should file in SparseArrays.jl?

Finally, just for completion, I decided to test the noise_rate_prototype as an Array with and without sparse capabilities:

oop_gprototype = eltype(oop_u0)[1 0 0 1
                                0 0 0 1]
oop_prob = SDEProblem{false}(f, g, oop_u0, tspan, noise_rate_prototype = oop_gprototype)
oop_sol = solve(oop_prob, EM(), dt = dt)
@btime solve($oop_prob, EM(), dt = $dt) # 7.620 μs (153 allocations: 13.27 KiB)

oop_gprototype = sparse(oop_gprototype)

oop_prob = SDEProblem{false}(f, g, oop_u0, tspan, noise_rate_prototype = oop_gprototype)
oop_sol = solve(oop_prob, EM(), dt = dt)
@btime solve($oop_prob, EM(), dt = $dt) # 13.332 μs (262 allocations: 25.86 KiB)

showing that the generated sparse object with oop_gprototype as a StaticArray or Array is the same object. Also, working with Array is better than with a SparseArray but worse than using a StaticArray.

The conclusion seems to be that the noise_rate_prototype that performs best for out of place problems is a StaticArray object.

In Place version

Using an Array as noise rate prototype yields to:

iip_gprototype = zeros(2, 4)

iip_prob = SDEProblem{true}(f, g, iip_u0, tspan, noise_rate_prototype = iip_gprototype)
iip_sol = solve(iip_prob, EM(), dt = dt)
@btime solve($iip_prob, EM(), dt = $dt) # 6.978 μs (101 allocations: 7.00 KiB)

Adding sparse capabilities:

iip_gprototype[1,1] = 1
iip_gprototype[1,4] = 1
iip_gprototype[2,4] = 1
iip_gprototype = sparse(iip_gprototype)

iip_prob = SDEProblem{true}(f, g, iip_u0, tspan, noise_rate_prototype = iip_gprototype)
iip_sol = solve(iip_prob, EM(), dt = dt)
@btime solve($iip_prob, EM(), dt = $dt) # 14.001 μs (222 allocations: 19.73 KiB)

worsen the situation, which seems odd. Is there an issue here?

The conclusion, for this particular case, seems to be that working with Arrays instead of SparseArrays is better...

These are the first MWE that I could think of. Please, let me know if you want me to elaborate more tests. For me, it is not clear if there are issues involved or everything is working correctly. Maybe that would requiere inspecting the operation g(u, p, t) * dW that are happening under the hood for all the cases involved.

Thanks!

[ChrisRackauckas]

worsen the situation, which seems odd. Is there an issue here?

The array isn't sparse, so that just adds a bunch of overhead.

It looks like the only issue is OOP?

[rvignolo]

@ChrisRackauckas, sorry for the delay regarding this issue!

Let me enumerate the issues that I see and might be relevant:

1. OOP: applying sparse to a StaticArray yields to an object with an underlying Array. In this sense, the noise_rate_prototype performs better if we keep it as an StaticArray instead of converting it to a SparseMatrixCSC1. This might be an issue for SparseArrays.jl.

2. IIP: I have shown that the solver performs faster if the noise_rate_prototype is kept as an Array instead of a SparseMatrixCSC1, which seems to be a problem as well. It should perform better with sparse matrix multiplication, right?

Thank you!

[ChrisRackauckas]

How sparse? You need like, really sparse, 0.01% filled, for sparse matrices to be "worth it". And they have to be big, like 100x100 or larger. Small sparse matrices are just slower because the indexing scheme is bad.

[rvignolo]

great, thanks!

In the OOP case, it is not clear to me if there is anything wrong besides the fact that applying sparse to a StaticArray returns an object with an Array. I mean, it is not clear that why do we need a noise_rate_prototype when using OOP since there no need to create a cache object. I suppose that the matrix-vector product g(u,p,t) * dW is Static product. Does this makes sense to you?

Finally... should I raise an issue in SparseArrays.jl?

Thanks!

[ChrisRackauckas]

I mean, it is not clear that why do we need a noise_rate_prototype when using OOP since there no need to create a cache object. I suppose that the matrix-vector product g(u,p,t) * dW is Static product. Does this makes sense to you?

We need to know the number of Brownian motions.

Finally... should I raise an issue in SparseArrays.jl?

I don't think so. Mixing sparse with static just doesn't make all that much sense, at least without a sparse static type (which would specialize on the zero locations)

[rvignolo]

We need to know the number of Brownian motions.

yes, that is correct, but I was thinking that since it is an OOP problem, we could just evaluate g(u, p, t) directly, i.e. we don't need du and get the number of brownian motions from there: size(g(u0, p, tspan[1]), 2).

Asking for the noise rate prototype is better in the sense that it will fail if g is not defined correctly.

I don't think so. Mixing sparse with static just doesn't make all that much sense, at least without a sparse static type (which would specialize on the zero locations)

This is what I suspected! That is why I wanted to ask you first.

Thank you!

[ChrisRackauckas]

we could just evaluate g(u, p, t) directly, i.e. we don't need du and get the number of brownian motions from there: size(g(u0, p, tspan[1]), 2).

Yes we could. That's an interesting proposal.

Something else I'd like to try is g(u,p,t,dW), which would allow users to write more efficient code for sparsity handling. However... if the relationship to dW isn't linear then this would allow users to write incorrect code that would not be an SDE... and violate the assumptions of the solver... very easily...