I am doing high dimensional Monte Carlo integrations in a maximum likelihood estimation using Optim.jl with the autodiff = :forward
option (based on ForwardDiff.jl). I use up to millions of quasirandom draws (e.g., Halton sequence and the like) for the numerical integration and thus require high precision summations for adding up the simulated numbers. However, current implementations of most of the high precision summation functions do not work with ForwardDiff, though I suppose they are fixable. I have limited knowledge of ForwardDiff and the related type constranits, thus I am asking for advice here.
I am considering the following functions (listed from slow to fast):

KahanSummation.sum_kbn()
: generic Julia code. 
psum_kbn()
: a multithreading extension ofKahanSummation.sum_kbn()
proposed by @goerch from this thread. 
AccurateArithmetic.sum_kbn()
andAccurateArithmetic.sum_oro()
: an implementation usingllvmcall
, which is close to the speed of the regularsum
.
Currently, only KahanSummation.sum_kbn()
, the slowest, is compatible with ForwardDiff.jl. According to the doc, ForwardDiff requires that the target function be composed of generic Julia functions. I am sure psum_kbn()
is generic Julia, and the problem is in type conversions. I am less certain of whether the other two functions from AccurateArithmetic are considered āgeneric Juliaā by ForwardDiff.
Here is a MWE.
using StatsFuns, Distributions, HaltonSequences, Optim
using KahanSummation, AccurateArithmetic
data = [0.755, 1.710, 0.891, 2.889, 0.881, 0.763, 0.365]
draws = 200 # potentially in millions for high dimension problems
halton = Halton(length=draws)
function LLexample(algo, log_Ļįµ¤Ā², log_Ļįµ„Ā², e, LDSdraw) # algo: summation algorithm
Ļįµ„ = exp(0.5*log_Ļįµ„Ā²)
Ļįµ¤ = exp(0.5*log_Ļįµ¤Ā²)
u = quantile(Normal(0, Ļįµ¤), 0.5 * LDSdraw .+ 0.5)
loglike = 0.0
for i = 1:length(e)
temp = normpdf.(0, Ļįµ„, e[i] .+ u)
loglike += log(algo(temp)/length(LDSdraw)) # algo: summation algorithm
end
return loglike
end
Using AccurateArithmetic.sum_kbn
(or similarly, Accurate
Arithmetic.sum_oro`), I got the following errors.
julia> func = TwiceDifferentiable(vars > LLexample(AccurateArithmetic.sum_kbn, vars[1], vars[2], data, halton), ones(2); autodiff = :forward);
julia> res = optimize(func, [0.1, 0.1], Newton())
ERROR: MethodError: no method matching accumulate(::Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#17#18", Float64}, Float64, 2}}}, ::AccurateArithmetic.Summation.var"#1#2"{typeof(AccurateArithmetic.EFT.fast_two_sum)}, ::Val{:scalar}, ::Val{2}, ::Val{0})
Closest candidates are:
accumulate(::Tuple{Vararg{AbstractArray{T}, A}}, ::F, ::Any, ::Val{Ushift}, ::Val{Prefetch}) where {F, A, T<:Union{Float32, Float64}, Ushift, Prefetch} at C:\Users\King\.julia\packages\AccurateArithmetic\449tA\src\Summation.jl:25
accumulate(::Tuple{Vararg{AbstractArray{T}, A}}, ::F, ::Any, ::Val{Ushift}) where {F, A, T<:Union{Float32, Float64}, Ushift} at C:\Users\King\.julia\packages\AccurateArithmetic\449tA\src\Summation.jl:25
accumulate(::Tuple{Vararg{AbstractArray{T}, A}}, ::F, ::Any) where {F, A, T<:Union{Float32, Float64}} at C:\Users\King\.julia\packages\AccurateArithmetic\449tA\src\Summation.jl:25
To test psum_kbn()
, which is a multithreading extension of KahanSummationās sum_kbn()
, the patch KahanSummation_patch.jl
may be needed, which can be downloaded here; see also the discussion here.
using InitialValues, Folds
include("KahanSummation_patch.jl") # for psum_kbn()
psum_kbn(f, X) = singleprec(Folds.mapreduce(f, InitialValues.asmonoid(plus_kbn), X)) # credit to @goerch
psum_kbn(X) = psum_kbn(identity, X)
Then,
julia> func = TwiceDifferentiable(vars > LLexample(psum_kbn, vars[1], vars[2], data, halton), ones(2); autodiff = :forward);
julia> res = optimize(func, [0.1, 0.1], Newton())
ERROR: MethodError: convert(::Type{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#16", Float64}, Float64, 2}}, ::TwicePrecisionN{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#16", Float64}, Float64, 2}}) is ambiguous. Candidates:
convert(::Type{T}, x::TwicePrecisionN) where T in Main at e:\temp\KahanSummation_patch.jl:74
convert(::Type{ForwardDiff.Dual{T, V, N}}, x) where {T, V, N} in ForwardDiff at C:\Users\King\.julia\packages\ForwardDiff\PBzup\src\dual.jl:384
Possible fix, define
convert(::Type{ForwardDiff.Dual{T, V, N}}, ::TwicePrecisionN) where
{T, V, N}
Suggestions to fix the problems will be appreciated.