Why does `solve` for `AbstractNoiseProblem` copies the noise?

[aleCombi]

Question❓

I looked at the implementation of the following function:

function DiffEqBase.__solve(prob::AbstractNoiseProblem,
        args::Union{Nothing, SciMLBase.DEAlgorithm}...; dt = 0.0,
        kwargs...)

The noise is copied at every call.
This seems to be the main bottleneck in performance.
I ran the following snippet using Geometric Brownian motion as an example:

using StochasticDiffEq, DiffEqNoiseProcess, BenchmarkTools, Profile

function simple_solve(r, σ, S₀, T)
    drift_part = (r - 0.5 * σ^2) * T
    diffusion_part = σ * sqrt(T)
    S = S₀ .* exp.(drift_part .+ diffusion_part .* randn())
    return S
end

r, σ, t₀, S₀ = 0.02, 0.2, 0.0, 1.0
T = 1.0
tspan = (t₀, T)
noise = GeometricBrownianMotionProcess(r, σ, t₀, S₀)
noise_problem = NoiseProblem(noise, tspan)

@btime StochasticDiffEq.solve(noise_problem; dt=T)
@btime simple_solve(r, σ, S₀, T)
@profile for i in 1:100000 StochasticDiffEq.solve(noise_problem; dt=T) end
Profile.print()

This is my output, suggesting that copying is the main performance reducer:

julia> include("noise_solve_speed.jl")
10.600 μs (47 allocations: 2.72 KiB)
55.172 ns (1 allocation: 16 bytes)
Overhead ╎ [+additional indent] Count File:Line; Function

╎357 @base\client.jl:557; _start()
╎ 357 @base\client.jl:338; exec_options(opts::Base.JLOptions)
╎ 357 @base\client.jl:421; run_main_repl(interactive::Bool, quiet::Bool, banner::Bool, history_file::Bool, color_set::Bool)
╎ 357 @base\essentials.jl:889; invokelatest
╎ 357 @base\essentials.jl:892; #invokelatest#2
╎ 357 @base\client.jl:437; (::Base.var"#1014#1016"{Bool, Bool, Bool})(REPL::Module)
╎ ╎ 357 @repl\src\REPL.jl:375; run_repl(repl::REPL.AbstractREPL, consumer::Any)
╎ ╎ 357 @repl\src\REPL.jl:389; run_repl(repl::REPL.AbstractREPL, consumer::Any; backend_on_current_task::Bool, backend::Any)
╎ ╎ 357 @repl\src\REPL.jl:228; kwcall(::NamedTuple, ::typeof(REPL.start_repl_backend), backend::REPL.REPLBackend, consumer::Any)
╎ ╎ 357 @repl\src\REPL.jl:231; start_repl_backend(backend::REPL.REPLBackend, consumer::Any; get_module::Function)
╎ ╎ 357 @repl\src\REPL.jl:246; repl_backend_loop(backend::REPL.REPLBackend, get_module::Function)
╎ ╎ ╎ 357 @repl\src\REPL.jl:150; eval_user_input(ast::Any, backend::REPL.REPLBackend, mod::Module)
╎ ╎ ╎ 357 @base\boot.jl:385; eval
╎ ╎ ╎ 357 @base\client.jl:494; include(fname::String)
╎ ╎ ╎ 357 @base\loading.jl:2206; _include(mapexpr::Function, mod::Module, _path::String)
╎ ╎ ╎ 357 @base\loading.jl:2146; include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
╎ ╎ ╎ ╎ 357 @base\boot.jl:385; eval
╎ ╎ ╎ ╎ 357 C:\repos\Hedgehog.jl\examples\noise_solve_speed.jl:42; top-level scope
╎ ╎ ╎ ╎ 357 @Profile\src\Profile.jl:44; macro expansion
3╎ ╎ ╎ ╎ 357 C:\repos\Hedgehog.jl\examples\noise_solve_speed.jl:18; macro expansion
9╎ ╎ ╎ ╎ 9 @base\boot.jl:622; NamedTuple
╎ ╎ ╎ ╎ 345 @DiffEqBase\src\solve.jl:1223; kwcall(::@NamedTuple{dt::Float64}, ::typeof(solve), ::NoiseProblem{NoiseProcess{Float64, 1, Float64, Flo…
╎ ╎ ╎ ╎ ╎ 345 @DiffEqBase\src\solve.jl:1224; #solve#49
╎ ╎ ╎ ╎ ╎ 345 @DiffEqNoiseProcess\src\solve.jl:1; __solve
1╎ ╎ ╎ ╎ ╎ 1 @DiffEqNoiseProcess\src\solve.jl:0; __solve(::NoiseProblem{NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoiseP…
╎ ╎ ╎ ╎ ╎ 339 @DiffEqNoiseProcess\src\solve.jl:7; __solve(::NoiseProblem{NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoiseP…
╎ ╎ ╎ ╎ ╎ 3 @DiffEqNoiseProcess\src\copy_noise_types.jl:39; copy
╎ ╎ ╎ ╎ ╎ 3 @DiffEqNoiseProcess\src\types.jl:189; NoiseProcess
╎ ╎ ╎ ╎ ╎ ╎ 1 @DiffEqNoiseProcess\src\types.jl:196; (NoiseProcess{false})(t0::Float64, W0::Float64, Z0::Nothing, dist::DiffEqNoiseProcess.Geome…
╎ ╎ ╎ ╎ ╎ ╎ 1 @ResettableStacks\src\core.jl:6; ResettableStack
╎ ╎ ╎ ╎ ╎ ╎ 1 @base\boot.jl:496; Array
1╎ ╎ ╎ ╎ ╎ ╎ 1 @base\boot.jl:477; Array
╎ ╎ ╎ ╎ ╎ ╎ 1 @DiffEqNoiseProcess\src\types.jl:198; (NoiseProcess{false})(t0::Float64, W0::Float64, Z0::Nothing, dist::DiffEqNoiseProcess.Geome…
╎ ╎ ╎ ╎ ╎ ╎ 1 @ResettableStacks\src\core.jl:6; ResettableStack
╎ ╎ ╎ ╎ ╎ ╎ 1 @base\boot.jl:496; Array
1╎ ╎ ╎ ╎ ╎ ╎ 1 @base\boot.jl:477; Array
╎ ╎ ╎ ╎ ╎ ╎ 1 @DiffEqNoiseProcess\src\types.jl:214; (NoiseProcess{false})(t0::Float64, W0::Float64, Z0::Nothing, dist::DiffEqNoiseProcess.Geome…
╎ ╎ ╎ ╎ ╎ ╎ 1 @base\array.jl:163; vect
1╎ ╎ ╎ ╎ ╎ ╎ 1 @base\boot.jl:477; Array
╎ ╎ ╎ ╎ ╎ 336 @DiffEqNoiseProcess\src\copy_noise_types.jl:44; copy
24╎ ╎ ╎ ╎ ╎ 24 @DiffEqNoiseProcess\src\copy_noise_types.jl:0; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNois…
31╎ ╎ ╎ ╎ ╎ 31 @DiffEqNoiseProcess\src\copy_noise_types.jl:3; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNois…
41╎ ╎ ╎ ╎ ╎ 41 @DiffEqNoiseProcess\src\copy_noise_types.jl:4; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNois…
13╎ ╎ ╎ ╎ ╎ 13 @DiffEqNoiseProcess\src\copy_noise_types.jl:5; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNois…
55╎ ╎ ╎ ╎ ╎ 55 @DiffEqNoiseProcess\src\copy_noise_types.jl:7; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNois…
40╎ ╎ ╎ ╎ ╎ 40 @DiffEqNoiseProcess\src\copy_noise_types.jl:9; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNois…
7╎ ╎ ╎ ╎ ╎ 10 @DiffEqNoiseProcess\src\copy_noise_types.jl:10; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
╎ ╎ ╎ ╎ ╎ ╎ 3 @RecursiveArrayTools\src\utils.jl:20; recursivecopy(a::Vector{Float64})
3╎ ╎ ╎ ╎ ╎ ╎ 3 @base\array.jl:411; copy
3╎ ╎ ╎ ╎ ╎ 3 @DiffEqNoiseProcess\src\copy_noise_types.jl:11; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
8╎ ╎ ╎ ╎ ╎ 12 @DiffEqNoiseProcess\src\copy_noise_types.jl:12; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
4╎ ╎ ╎ ╎ ╎ ╎ 4 @base\Base.jl:37; getproperty(x::ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, f::Symbol)
17╎ ╎ ╎ ╎ ╎ 18 @DiffEqNoiseProcess\src\copy_noise_types.jl:13; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
1╎ ╎ ╎ ╎ ╎ ╎ 1 @base\Base.jl:37; getproperty(x::ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, f::Symbol)
18╎ ╎ ╎ ╎ ╎ 52 @DiffEqNoiseProcess\src\copy_noise_types.jl:14; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
1╎ ╎ ╎ ╎ ╎ ╎ 1 @base\Base.jl:37; getproperty(x::ResettableStacks.ResettableStack{Tuple{Float64, Float64, Nothing}, false}, f::Symbol)
╎ ╎ ╎ ╎ ╎ ╎ 33 @RecursiveArrayTools\src\utils.jl:13; recursivecopy(a::Vector{Tuple{Float64, Float64, Nothing}})
╎ ╎ ╎ ╎ ╎ ╎ 33 @base\deepcopy.jl:26; deepcopy
2╎ ╎ ╎ ╎ ╎ ╎ 2 @base\deepcopy.jl:0; deepcopy_internal(x::Vector{Tuple{Float64, Float64, Nothing}}, stackdict::IdDict{Any, Any})
╎ ╎ ╎ ╎ ╎ ╎ 3 @base\deepcopy.jl:89; deepcopy_internal(x::Vector{Tuple{Float64, Float64, Nothing}}, stackdict::IdDict{Any, Any})
╎ ╎ ╎ ╎ ╎ ╎ 3 @base\abstractdict.jl:17; haskey
╎ ╎ ╎ ╎ ╎ ╎ 3 @base\iddict.jl:192; in
2╎ ╎ ╎ ╎ ╎ ╎ ╎ 3 @base\iddict.jl:102; get
╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @base\Base.jl:37; getproperty
╎ ╎ ╎ ╎ ╎ ╎ 12 @base\deepcopy.jl:92; deepcopy_internal(x::Vector{Tuple{Float64, Float64, Nothing}}, stackdict::IdDict{Any, Any})
╎ ╎ ╎ ╎ ╎ ╎ 12 @base\deepcopy.jl:97; _deepcopy_array_t
8╎ ╎ ╎ ╎ ╎ ╎ 8 @base\array.jl:411; copy
2╎ ╎ ╎ ╎ ╎ ╎ 2 @base\iddict.jl:86; setindex!(d::IdDict{Any, Any}, val::Any, key::Any)
╎ ╎ ╎ ╎ ╎ ╎ 1 @base\iddict.jl:91; setindex!(d::IdDict{Any, Any}, val::Any, key::Any)
1╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @base\Base.jl:37; getproperty
1╎ ╎ ╎ ╎ ╎ ╎ 1 @base\iddict.jl:96; setindex!(d::IdDict{Any, Any}, val::Any, key::Any)
╎ ╎ ╎ ╎ ╎ ╎ 16 @base\iddict.jl:48; IdDict
9╎ ╎ ╎ ╎ ╎ ╎ 16 @base\iddict.jl:30; IdDict
7╎ ╎ ╎ ╎ ╎ ╎ 7 @base\boot.jl:477; Array
1╎ ╎ ╎ ╎ ╎ 1 @DiffEqNoiseProcess\src\copy_noise_types.jl:15; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
4╎ ╎ ╎ ╎ ╎ 4 @DiffEqNoiseProcess\src\copy_noise_types.jl:16; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
3╎ ╎ ╎ ╎ ╎ 3 @DiffEqNoiseProcess\src\copy_noise_types.jl:17; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
5╎ ╎ ╎ ╎ ╎ 5 @DiffEqNoiseProcess\src\copy_noise_types.jl:18; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
1╎ ╎ ╎ ╎ ╎ 1 @DiffEqNoiseProcess\src\copy_noise_types.jl:24; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
1╎ ╎ ╎ ╎ ╎ 23 @DiffEqNoiseProcess\src\copy_noise_types.jl:25; copy!(Wnew::NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoi…
╎ ╎ ╎ ╎ ╎ ╎ 22 @randomnumbers\src\Xorshifts\xoroshiro128.jl:55; copy(src::RandomNumbers.Xorshifts.Xoroshiro128Plus)
╎ ╎ ╎ ╎ ╎ ╎ 21 @randomnumbers\src\Xorshifts\xoroshiro128.jl:28; RandomNumbers.Xorshifts.Xoroshiro128Plus()
╎ ╎ ╎ ╎ ╎ ╎ 21 @randomnumbers\src\utils.jl:24; gen_seed
╎ ╎ ╎ ╎ ╎ ╎ 9 @base\tuple.jl:391; Tuple
9╎ ╎ ╎ ╎ ╎ ╎ 9 @base\tuple.jl:428; _totuple
╎ ╎ ╎ ╎ ╎ ╎ 12 @random\src\Random.jl:293; rand
╎ ╎ ╎ ╎ ╎ ╎ 12 @random\src\Random.jl:290; rand
╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @base\boot.jl:494; Array
╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @base\boot.jl:486; Array
1╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @base\boot.jl:477; Array
╎ ╎ ╎ ╎ ╎ ╎ ╎ 11 @random\src\Random.jl:269; rand!
╎ ╎ ╎ ╎ ╎ ╎ ╎ 11 @random\src\RNGs.jl:32; rand!
11╎ ╎ ╎ ╎ ╎ ╎ ╎ 11 @base\libc.jl:384; getrandom!
╎ ╎ ╎ ╎ ╎ ╎ 1 @randomnumbers\src\Xorshifts\xoroshiro128.jl:51; copyto!
1╎ ╎ ╎ ╎ ╎ ╎ 1 @base\Base.jl:41; setproperty!
╎ ╎ ╎ ╎ ╎ 4 @DiffEqNoiseProcess\src\solve.jl:19; __solve(::NoiseProblem{NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoise…
1╎ ╎ ╎ ╎ ╎ 1 @base\Base.jl:0; setup_next_step!
╎ ╎ ╎ ╎ ╎ 1 …rc\noise_interfaces\noise_process_interface.jl:149; setup_next_step!
╎ ╎ ╎ ╎ ╎ 1 @base\operators.jl:378; >
1╎ ╎ ╎ ╎ ╎ ╎ 1 @base\float.jl:536; <
╎ ╎ ╎ ╎ ╎ 2 …rc\noise_interfaces\noise_process_interface.jl:174; setup_next_step!
╎ ╎ ╎ ╎ ╎ 2 @ResettableStacks\src\core.jl:35; copyat_or_push!
╎ ╎ ╎ ╎ ╎ ╎ 2 @base\array.jl:1119; push!
2╎ ╎ ╎ ╎ ╎ ╎ 2 @base\array.jl:1072; _growend!
╎ ╎ ╎ ╎ ╎ 1 @DiffEqNoiseProcess\src\solve.jl:24; __solve(::NoiseProblem{NoiseProcess{Float64, 1, Float64, Float64, Nothing, Nothing, DiffEqNoise…
╎ ╎ ╎ ╎ ╎ 1 …src\noise_interfaces\noise_process_interface.jl:13; accept_step!
╎ ╎ ╎ ╎ ╎ 1 …rc\noise_interfaces\noise_process_interface.jl:30; accept_step!
╎ ╎ ╎ ╎ ╎ ╎ 1 @base\array.jl:1119; push!
1╎ ╎ ╎ ╎ ╎ ╎ 1 @base\array.jl:1072; _growend!
Total snapshots: 359. Utilization: 100% across all threads and tasks. Use the groupby kwarg to break down by thread and/or task.

Cleary my case is not general as only the final point in time is taken.
Is there a possibility to avoid this copying? This creates a lot of allocation especially with Ensemble solvers.

[ChrisRackauckas]

If you don't copy then solving twice would give you the same result. There's a way to turn it off but we found that's the behavior most users expected.

[aleCombi]

Is that the case also for EnsembleProblem defined on a NoiseProblem? Or is there a best practice to make it faster?

[ChrisRackauckas]

If each of the noise problem solves are writing to the same noiseprocess array, then they will collide. At some point you need to copy it to new arrays if you're going to get many solutions of it.

[aleCombi]

That makes sense, I am just wondering why it is more expensive than the "manual implementation".
Perhaps my example is too elementary...

[ChrisRackauckas]

S = S₀ .* exp.(drift_part .+ diffusion_part .* randn()) you're missing a cumulative sum? It's not valid to generate all of the brownian at once, that would be a jagged path.

[aleCombi]

I wanted to test just the terminal value generation rather than the whole path. That's also why I took dt=T is the NoiseProblem solver.
Perhaps that's not a good use case of a NoiseProblem?

I noticed that solve(NoiseProblem,...) scales very well with N if dt=T/N, which makes sense a copying is done once per trajectory.

[ChrisRackauckas]

NoiseProblem is mostly about generating the trajectories. If you just want to get to the end, you can skip taking steps and directly use those statistics instead of time stepping by making dt=T yes. But if you can do that, then you also don't need to generate new arrays because you just have the scalar end value. So that is something we don't have a good interface to be fully performant... we might want to make an interface for that.

[aleCombi]

That makes sense, thanks for the clarification! I'd love to contribute if that is deemed interesting at some point.

[ChrisRackauckas]

It's definitely interesting! There's just only so many hours in a day, so a contribution would be amazing.