I’m experimenting with using Turing to handle some Gaussian random fields problems, and I’m struggling a bit with the interface. The following is supposed to do a Bayesian posterior sampling of the problem
y= u(x_i) + \eta_i
where u has a prior mean zero Gaussian with covariance (-d^2/dx^2)^{-1}. I’m trying to use some old code I wrote for FFT based sampling, which might be the problem. Here’s what I’m running:
using Turing
using Random
using Distributions
using FFTW
N = 32;
γ = 0.01;
x = LinRange(0,1,N+1)[2:end-1];
uᵗ(x) = x*(1-x);
Random.seed!(500); # set a seed for reproducibility
y_data = @. uᵗ(x) + γ * randn()
"""
`build_field` - Build a mean zero Gaussian random field with the (-Δ)^{-α} covariance operator in dimension one
### Fields
`ξ` - Vector of N(0,1) values
`α=1` - Smoothness parameter
"""
function build_field(ξ; α=1)
N = length(ξ)
uhat = zeros(ComplexF64,2*N); # preallocate space
# construct the eigenvalues
k = 1:N;
λ = @. 1/(π*k)^(2*α);
# fill in the nonzero entries
# NOTE we need to multiply by 2 *N for FFT scaling
@. uhat[2:N+1] = 2 * N * sqrt(λ) * sqrt(2) * ξ;
# invert and get the relevant imaginary part
u = imag.(ifft(uhat))[N+2:end];
return u
end
@model function mean_recovery(y_data)
ξ ~ MvNormal(zeros(N), 1.)
u = build_field(ξ);
for i in 1:length(y_data)
y_data[i]~Normal(u[i], γ^2)
end
end
model=mean_recovery(y_data)
chain = sample(model, HMC(0.1, 10), 10^4)
The code then seems to run for a few moments before spitting out:
TypeError: in new, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 11}
Stacktrace:... # omitted for brevity, but I can put in if requested.
You’ve constructed this to have a ComplexF64 eltype, but this is incompatible with operator-overloading automatic differentiation (AD) packages like ForwardDiff, ReverseDiff, and Tracker, which use a custom numeric type to compute gradients. Instead, try
uhat = zeros(complex(eltype(ξ)),2*N);
The next issue you’ll have is that not all AD packages support complex numbers. Zygote definitely does, but it doesn’t support mutation. But if you want to differentiate an FFT, I think Zygote is your only option. So here’s a non-mutating build_field:
function build_field(ξ; α=one(eltype(ξ)))
N = length(ξ)
# construct the eigenvalues
πk = π * (1:N);
# NOTE we need to multiply by 2 *N for FFT scaling
c = 2N * sqrt(2)
umid = @. c * ξ / πk^α;
uhat = [0; umid; zeros(N - 1)]
# invert and get the relevant imaginary part
u = @views imag.(ifft(uhat)[N+2:end]);
return u
end
See https://turing.ml/dev/docs/using-turing/autodiff for details.
Also, from your simulated data, it seems γ might be intended to be a standard deviation? Note that Normal takes a standard deviation as a parameter, so perhaps you mean to use Normal(u[i], γ) as your likelihood?
1 Like
I tried your suggested version of build_field, but now this generates the error:
type DataType has no field mutable
Stacktrace:
[1] getproperty
@ ./Base.jl:37 [inlined]
[2] adjoint
@ ~/.julia/packages/Zygote/zowrf/src/lib/lib.jl:281 [inlined]
[3] _pullback
@ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:65 [inlined]
# keeps going for a bit
But the switch to Zygote does appear to have addressed the other issue.
Also, are any of the interpolation packages compatible with Zygote for automatic differentiation?
Odd, this runs fine for me:
using Turing, Random, Distributions, FFTW, Zygote
Turing.setadbackend(:zygote)
N = 32;
γ = 0.01;
x = LinRange(0,1,N+1)[2:end-1];
uᵗ(x) = x*(1-x);
Random.seed!(500); # set a seed for reproducibility
y_data = @. uᵗ(x) + γ * randn()
"""
`build_field` - Build a mean zero Gaussian random field with the (-Δ)^{-α} covariance operator in dimension one
### Fields
`ξ` - Vector of N(0,1) values
`α=1` - Smoothness parameter
"""
function build_field(ξ; α=one(eltype(ξ)))
N = length(ξ)
# construct the eigenvalues
πk = π * (1:N);
# NOTE we need to multiply by 2 *N for FFT scaling
c = 2N * sqrt(2)
umid = @. c * ξ / πk^α;
uhat = [0; umid; zeros(N - 1)]
# invert and get the relevant imaginary part
u = @views imag.(ifft(uhat)[N+2:end]);
return u
end
@model function mean_recovery(y_data)
ξ ~ filldist(Normal(), N)
u = build_field(ξ);
y_data .~ Normal.(u, γ)
end
model=mean_recovery(y_data)
chain = sample(model, NUTS(), MCMCThreads(), 10^3, 4)
Note that I switched to NUTS and dropped the number of requested samples. This model samples so well with NUTS that the 4000 draws for most parameters have effective sample sizes of over 10,000.
I’m not familiar with the interpolation packages. I’d say just try and see.
Never mind. I think one of my packages was out of date. Interpolations appears to work.