Type error for primal feasibility report in non-Float64 precision

[mtanneau]

using JuMP

model = JuMP.GenericModel{Float32}()
@variable(model, x)
# linear constraints are OK
@constraint(model, x == 1f0)
# quadratic constraints are OK too (it seems)
@constraint(model, x*x == 0f0)

JuMP.primal_feasibility_report(model, Dict(x => 0f0))  # no error so far

# general nonlinear constraints throw an error
# note: a similar error occurs 
@constraint(model, cos(x) == 1f0)
@constraint(model, sqrt(x) == 1f0)

JuMP.primal_feasibility_report(model, Dict(x => 0f0))  # 💥

the stacktrace below suggests that something got converted

ERROR: distance_to_set using the distance metric MathOptInterface.Utilities.ProjectionUpperBoundDistance() for set type MathOptInterface.EqualTo{Float32} has not been implemented yet.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] distance_to_set(d::MathOptInterface.Utilities.ProjectionUpperBoundDistance, ::Float64, set::MathOptInterface.EqualTo{Float32})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/PUD05/src/Utilities/distance_to_set.jl:75
  [3] distance_to_set(point::Float64, set::MathOptInterface.EqualTo{Float32})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/PUD05/src/Utilities/distance_to_set.jl:71
  [4] _distance_to_set(point::Float64, set::MathOptInterface.EqualTo{Float32}, ::Type)
    @ JuMP ~/.julia/packages/JuMP/CU7H5/src/feasibility_checker.jl:172
  [5] _add_infeasible_constraints(model::GenericModel{Float32}, ::Type{GenericNonlinearExpr{GenericVariableRef{Float32}}}, ::Type{MathOptInterface.EqualTo{Float32}}, violated_constraints::Dict{Any, Float32}, point_f::Function, atol::Float32)
    @ JuMP ~/.julia/packages/JuMP/CU7H5/src/feasibility_checker.jl:143
  [6] primal_feasibility_report(point::Function, model::GenericModel{Float32}; atol::Float32, skip_missing::Bool)
    @ JuMP ~/.julia/packages/JuMP/CU7H5/src/feasibility_checker.jl:107
  [7] primal_feasibility_report
    @ ~/.julia/packages/JuMP/CU7H5/src/feasibility_checker.jl:99 [inlined]
  [8] #primal_feasibility_report#127
    @ ~/.julia/packages/JuMP/CU7H5/src/feasibility_checker.jl:58 [inlined]
  [9] primal_feasibility_report(model::GenericModel{Float32}, point::Dict{GenericVariableRef{Float32}, Float32})
    @ JuMP ~/.julia/packages/JuMP/CU7H5/src/feasibility_checker.jl:52
 [10] top-level scope
    @ REPL[152]:1

Notes:

• A similar issue occurs with other non-Float64 arithmetic
• No error is thrown when using Float64 arithmetic (see snippet below)

---

Same code as above, but Float64

using JuMP

model = JuMP.GenericModel{Float64}()
@variable(model, x)
# linear constraints are OK
@constraint(model, x == 1.0)
# quadratic constraints are OK too (it seems)
@constraint(model, x*x == 0.0)

JuMP.primal_feasibility_report(model, Dict(x => 0.0))  # no error so far

# general nonlinear constraints don't throw an error in Float64
@constraint(model, cos(x) == 1.0)
@constraint(model, sqrt(x) == 1.0)

JuMP.primal_feasibility_report(model, Dict(x => 0.0))

[odow]

Do you have a motivation for nonlinear models with non-Float64 precision?

[mtanneau]

I was using this to gauge how working in a lower precision (Float32) might affect the accuracy of constraint violations.
This is motivated by the fact that most machine learning tools use Float32 precision by default.

More details: consider a constraint h(x) = 0, and consider h64 and h32 the two methods using 64-bit and 32-bit precision.
Now consider a candidate solution x, typically the output of a machine learning model, hence in 32-bit precision.
The underlying question is: how does h32(x) compare with h64(Float64(x))?

This brought me to JuMP's generic model interface, which allows me to build the same model in different precision, then evaluate the same solution via JuMP.primal_feasibility_report

[odow]

Don't solver tolerances make much more of a difference?

[mtanneau]

Don't solver tolerances make much more of a difference?

Depends on the use case I guess.
Solvers' tolerances are usually in the range of sqrt(eps) where eps is the machine precision for the chosen arithmetic.

If we're trying to run Ipopt in Float32 precision, then I would expect the numerics to play a bigger role than tolerances (I've seen that with Tulip on LPs, where Float32 would often crash in later iterations).

The setting I'm considering is where one is training a machine learning model to predict good solutions. This typically involves evaluating the feasibility of said solutions, and everything is in 32-bit precision. So I was exploring whether I could trust the violation levels when evaluated in 32-bit precision... which led me here.

One possible example: consider the constraint x + y - z == 0, with x, y, z = 1f4, 1f-4, 1f4.
In Float32, this will evaluate to 0f0 and the constraint is satisfied. But if converted to Float64, It would evaluate to 1e-4, which would violate the usual 1e-6 or 1e-8 tolerances.

(obviously we're now past the original issue 🙃)

[odow]

One possible example

Isn't this example a counter example to

whether I could trust the violation levels when evaluated in 32-bit precision

If so, the answer is "no".