具有指定噪声的随机微分方程灵敏度分析
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了具有指定噪声的随机微分方程灵敏度分析相关的知识,希望对你有一定的参考价值。
鉴于噪声的特定实现,我正在尝试计算随机微分方程(SDE)解决方案的函数的梯度。如果未指定噪声,则可以成功计算这些梯度,如DiffEqFlux.jl: Using Other Differential Equations所示。我也可以成功获得针对特定噪声实现的SDE解决方案,如DifferentialEquations.jl: NoiseWrapper Example中所示。但是,当我尝试将两者放在一起时,代码将返回错误。
这是从上面引用的两个单独示例改编而成的最小工作示例:
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]
...
紧随其后的是栈跟踪的一个无休止的结论(如果您认为这会有所帮助,我可以将其发布,但是由于它比该问题的其余部分还要长,所以我不想把事情弄得一团糟。
是否正在为当前不支持指定噪声实现的SDE问题计算梯度,或者我只是不进行适当的函数调用?我很容易相信后者,因为要努力达到上述代码的工作部分的工作有点费劲,但是我找不到任何有关逐步完成此操作后错误提供的线索。 Juno调试器编写代码。
作为StackOverflow解决方案,您可以使用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跟踪开发情况。它可能会在一小时左右的时间内解决。
编辑:已发布此修复程序,并且您的原始代码有效。
以上是关于具有指定噪声的随机微分方程灵敏度分析的主要内容,如果未能解决你的问题,请参考以下文章