Possible bug in automatic differentiation

[jywu20]

Consider the following demo, where we keep a matrix to be semidefinite during the optimization:

using NonconvexSemidefinite, NonconvexIpopt, LinearAlgebra, Test
using Distributions, DistributionsAD, ChainRulesTestUtils, Random

K = 3
g = 1

Mij(x⁴, x², n) = begin
    if n == 0
        return 1.0
    end
    if isodd(n)
        return 0.0
    end
    if n == 2
        return x²
    end
    if n == 4
        return x⁴
    end
    
    Mij(x⁴, x², n - 2) * x² + Mij(x⁴, x², n - 4)
end

sd_constraint((x⁴, x²)) = [Mij(x⁴, x², i + j) for i in 0 : K, j in 0 : K]

# Objective function
f((x⁴, x²)) = 2x² + 3x⁴

model = Model(f)

x0 = [1.0, 0.8]
lbs = [0.0, 0.0]
ubs = [Inf, Inf]
addvar!(model, lbs, ubs)

add_sd_constraint!(model, sd_constraint)

alg = SDPBarrierAlg(sub_alg=IpoptAlg())
options = SDPBarrierOptions(sub_options=IpoptOptions(max_iter=400))
result = optimize(model, alg, x0, options = options)

minimum, minimizer, optimal_ind = result.minimum, result.minimizer, result.optimal_ind
m_mat_minimizer = sd_constraint(minimizer)

@show minimum
@show eigen(m_mat_minimizer).values

The demo works well and finds the expected solution (just letting x⁴ == 0.0 and x² == 0.0). However, the following code

using NonconvexSemidefinite, NonconvexIpopt, LinearAlgebra, Test
using Distributions, DistributionsAD, ChainRulesTestUtils, Random

K = 3
g = 1

Mij(x⁴, x², n) = begin
    if n == 0
        return 1.0
    end
    if isodd(n)
        return 0.0
    end
    if n == 2
        return x²
    end
    if n == 4
        return x⁴
    end
    
    (n - 4) * Mij(x⁴, x², n - 2) * x² + Mij(x⁴, x², n - 4)
end

sd_constraint((x⁴, x²)) = [Mij(x⁴, x², i + j) for i in 0 : K, j in 0 : K]

# Objective function
f((x⁴, x²)) = 2x² + 3x⁴

model = Model(f)

x0 = [1.0, 0.8]
lbs = [0.0, 0.0]
ubs = [Inf, Inf]
addvar!(model, lbs, ubs)

add_sd_constraint!(model, sd_constraint)

alg = SDPBarrierAlg(sub_alg=IpoptAlg())
options = SDPBarrierOptions(sub_options=IpoptOptions(max_iter=400))
result = optimize(model, alg, x0, options = options)

minimum, minimizer, optimal_ind = result.minimum, result.minimizer, result.optimal_ind
m_mat_minimizer = sd_constraint(minimizer)

@show minimum
@show eigen(m_mat_minimizer).values

throws a strange error, saying

LoadError: MethodError: no method matching getindex(::Nothing, ::Int64)

The only difference between the two versions is that the second version has an additional (n - 4) factor at the end of the definition of Mij. It can be verified by Mathematica that the second problem also has a solution at x⁴ == 0.0 and x² == 0.0.

Considering that the error message contains lines like

 (::Zygote.var"#609#612"{Array{Union{Nothing, Tuple{Float64,Float64}},2},Tuple{UnitRange{Int64},UnitRange{Int64}}})(::Int64) 

it seems there is some hidden bug in Zygote.

[mohdibntarek]

I will look into it but opening an issue in Zygote can also get some people's attention.

[jywu20]

I want to, but since I don't know much about how NonconvexSemidefinite.jl is implemented and how Zygote is called, I'm not sure what to ask... Just copying this issue, or trying to identify the part concerning Zygote?

[mohdibntarek]

Just copying this issue, or trying to identify the part concerning Zygote?

Preferably the latter but if you can't, just copy the issue.

[mohdibntarek]

as an alternative for now, considering using ForwardDiff for Mij.

[mohdibntarek]

Something like this will do the trick.

using AbstractDifferentiation, ForwardDiff, ChainRulesCore

function ChainRulesCore.rrule(::typeof(Mij), x⁴, x², n)
   val, ds = AD.value_and_derivative(AD.ForwardDiffBackend(), Mij, x⁴, x², n)
   return val, Δ -> (NoTangent(), (Δ .* ds)...)
end