Can't figure out why log-likelihood function doesn't work with autodiff?

I’m trying to practice Julia by coding an ordinal regression model. I have defined a log-likelihood function and it works as expected when I used Optim.jl for optimization. To gain additional speed, I want to try using the autodiff features in Julia. However, they don’t seem to work with my log-likelihood function and I can’t figure out why. Any help is appreciated.

Log-Likelihood code

function neg_loglik(params::Vector, X::Matrix, Y_ord::Vector)

    # Define the link function
    F(x) = cdf(Normal(0,1), x)

    # Number of categories
    K = maximum(Y_ord)

    # Define the subset of the params vector for the cutpoints and the betas
    c = params[1:(K-1)]
    β = params[K:end]

    # Number of observations
    n = length(Y_ord)

    # Initialize a vector to store likelihood values
    L = zeros(n)

    # Return Inf negative log-likelihood if the cutpoints are out of order
    if !issorted(c)
        return(Inf)
    end

    for i in 1:n
        # Current observation
        y = Y_ord[i]
        # Linear predictor
        η = dot(X[i,:], β)

        if y == 1
            L[i] = F(c[1] - η)
        elseif y == K
            L[i] = 1 - F(c[K - 1] - η)
        else
            L[i] = F(c[y] - η) - F(c[y-1] - η)
        end
    end

    ℓ = log.(L)

    return(-sum(ℓ))
end

neg_loglik_c(params::Vector) = neg_loglik(params, X, Y_ord)

Trying autodiff

ForwardDiff.gradient(neg_loglik_c, vcat(sort!(rand(K)), zeros(p)))

Error message

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
 [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7})
   @ Base ./number.jl:7
 [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}, i1::Int64)
   @ Base ./array.jl:966
 [3] neg_loglik(params::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}}, X::Matrix{Float64}, Y_ord::Vector{Int64})
   @ Main ./code:36
 [4] neg_loglik_c(params::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}})
   @ Main ./code:88
 [5] chunk_mode_gradient(f::typeof(neg_loglik_c), x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}}})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/gradient.jl:123
 [6] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}}}, ::Val{true})
   @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/gradient.jl:21
 [7] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}}}) (repeats 2 times)
   @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/gradient.jl:17
 [8] top-level scope
   @ REPL[7]:1

Usually that error comes up because you’re trying to convert something to Float64 somewhere in your code that needs to accept Dual numbers to work with ForwardDiff. Not obvious to me in your code where that’s happening, but I would look for that.

3 Likes

It’s the zeros(n). Try instead zeros(eltype(params), n) and see how it goes.

3 Likes