具有指定噪声的随机微分方程灵敏度分析

Question

我正在尝试计算给定噪声的特定实现的随机微分方程 (SDE) 解的泛函的梯度。如果不指定噪声，我可以成功计算出这些梯度，如 DiffEqFlux.jl: Using Other Differential Equations. I can also successfully obtain the solution to my SDE for a specific noise realization, like shown in DifferentialEquations.jl: NoiseWrapper Example 所示。但是，当我尝试将两者放在一起时，代码 returns 出错了。

这是根据上面引用的两个独立示例改编的最小工作示例：

using StochasticDiffEq, DiffEqBase, DiffEqNoiseProcess, DiffEqSensitivity, Zygote

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
function lotka_volterra_noise(du,u,p,t)
  du[1] = 0.1u[1]
  du[2] = 0.1u[2]
end
dt = 1//2^(4)
u0 = [1.0,1.0]
p = [2.2, 1.0, 2.0, 0.4]
prob1 = SDEProblem(lotka_volterra,lotka_volterra_noise,u0,(0.0,10.0),p)
sol1 = solve(prob1,EM(),dt=dt,save_noise=true)

W2 = NoiseWrapper(sol1.W)
prob2 = SDEProblem(lotka_volterra,lotka_volterra_noise,u0,(0.0,10.0),p,noise=W2)
sol2 = solve(prob2,EM(),dt=dt)

function predict_sde1(p)
  Array(concrete_solve(remake(prob1,p=p),EM(),dt=dt,sensealg=ForwardDiffSensitivity(),saveat=0.1))
end
loss_sde1(p) = sum(abs2,x-1 for x in predict_sde1(p))

loss_sde1(p)

# This gradient is successfully calculated
Zygote.gradient(loss_sde1,p)

function predict_sde2(p)
  W2 = NoiseWrapper(sol1.W)
  Array(concrete_solve(remake(prob2,p=p,noise=W2),EM(),dt=dt,sensealg=ForwardDiffSensitivity(),saveat=0.1))
end
loss_sde2(p) = sum(abs2,x-1 for x in predict_sde2(p))

# This loss is successfully calculated
loss_sde2(p)

# This gradient calculation raises and error
Zygote.gradient(loss_sde2,p)

我在运行这段代码结束时得到的错误是

TypeError: in setfield!, expected Float64, got ForwardDiff.Dual{Nothing,Float64,4}

Stacktrace:
 [1] setproperty! at ./Base.jl:21 [inlined]
...

接着是对堆栈跟踪的没完没了的结论（如果您认为它有帮助，我可以 post 它，但由于它比这个问题的其余部分更长，我不想把事情搞得一团糟蝙蝠）。

当前不支持使用指定噪声实现为 SDE 问题计算梯度，还是我只是没有进行适当的函数调用？我很容易相信后者，因为要达到上述代码的工作部分工作的程度有点困难，但是我找不到任何关于我在逐步完成后提供的错误信息的线索使用 Juno 调试器编写代码。

Answer 1

作为 Whosebug 解决方案，您可以使用 ForwardDiffSensitivity(convert_tspan=false) 解决此问题。工作代码：

using StochasticDiffEq, DiffEqBase, DiffEqNoiseProcess, DiffEqSensitivity, Zygote

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
function lotka_volterra_noise(du,u,p,t)
  du[1] = 0.1u[1]
  du[2] = 0.1u[2]
end
dt = 1//2^(4)
u0 = [1.0,1.0]
p = [2.2, 1.0, 2.0, 0.4]
prob1 = SDEProblem(lotka_volterra,lotka_volterra_noise,u0,(0.0,10.0),p)
sol1 = solve(prob1,EM(),dt=dt,save_noise=true)

W2 = NoiseWrapper(sol1.W)
prob2 = SDEProblem(lotka_volterra,lotka_volterra_noise,u0,(0.0,10.0),p,noise=W2)
sol2 = solve(prob2,EM(),dt=dt)

function predict_sde1(p)
  Array(concrete_solve(remake(prob1,p=p),EM(),dt=dt,sensealg=ForwardDiffSensitivity(convert_tspan=false),saveat=0.1))
end
loss_sde1(p) = sum(abs2,x-1 for x in predict_sde1(p))

loss_sde1(p)

# This gradient is successfully calculated
Zygote.gradient(loss_sde1,p)

function predict_sde2(p)
  Array(concrete_solve(prob2,EM(),prob2.u0,p,dt=dt,sensealg=ForwardDiffSensitivity(convert_tspan=false),saveat=0.1))
end
loss_sde2(p) = sum(abs2,x-1 for x in predict_sde2(p))

# This loss is successfully calculated
loss_sde2(p)

# This gradient calculation raises and error
Zygote.gradient(loss_sde2,p)

作为开发人员...这不是一个很好的解决方案，我们的默认设置应该更好。我会努力的。您可以在此处跟踪开发 https://github.com/JuliaDiffEq/DiffEqSensitivity.jl/issues/204。大概一个小时左右就能搞定。

编辑：修复已发布，您的原始代码有效。

具有指定噪声的随机微分方程灵敏度分析

Stochastic differential equation sensitivity analysis with specified noise

stochastic

julia

differentialequations.jl