AssertionError: k == rinfo.nnz + N

[basberr]

Hi all,

After the update from v1.46 to v.1.47, some of my codes now trigger this error.

This might be related to this recent change: #2882
(which btw has sped up the resolution of some of my problems by a factor ~2, so that's very nice!)

I did not manage to isolate a MWE yet but I verified that my code works with v.1.46.

Here is the stacktrace in case it evokes something to you (it's with Ipopt but I get the same with MadNLP)

ERROR: AssertionError: k == rinfo.nnz + N Stacktrace: [1] _indirect_recover_structure(rinfo::MathOptInterface.Nonlinear.ReverseAD.Coloring.RecoveryInfo) @ MathOptInterface.Nonlinear.ReverseAD.Coloring ~/.julia/packages/MathOptInterface/9lCV5/src/Nonlinear/ReverseAD/Coloring/Coloring.jl:447 [2] hessian_color_preprocess(edgelist::Set{…}, num_total_var::Int64, seen_idx::MathOptInterface.Nonlinear.ReverseAD.Coloring.IndexedSet) @ MathOptInterface.Nonlinear.ReverseAD.Coloring ~/.julia/packages/MathOptInterface/9lCV5/src/Nonlinear/ReverseAD/Coloring/Coloring.jl:489 [3] MathOptInterface.Nonlinear.ReverseAD._FunctionStorage(nodes::Vector{…}, const_values::Vector{…}, num_variables::Int64, coloring_storage::MathOptInterface.Nonlinear.ReverseAD.Coloring.IndexedSet, want_hess::Bool, subexpressions::Vector{…}, dependent_subexpressions::Vector{…}, subexpression_linearity::Vector{…}, subexpression_edgelist::Vector{…}, subexpression_variables::Vector{…}, moi_index_to_consecutive_index::Dict{…}) @ MathOptInterface.Nonlinear.ReverseAD ~/.julia/packages/MathOptInterface/9lCV5/src/Nonlinear/ReverseAD/types.jl:97 [4] initialize(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, requested_features::Vector{Symbol}) @ MathOptInterface.Nonlinear.ReverseAD ~/.julia/packages/MathOptInterface/9lCV5/src/Nonlinear/ReverseAD/mathoptinterface_api.jl:123 [5] initialize(evaluator::MathOptInterface.Nonlinear.Evaluator{…}, features::Vector{…}) @ MathOptInterface.Nonlinear ~/.julia/packages/MathOptInterface/9lCV5/src/Nonlinear/evaluator.jl:105 [6] _setup_model(model::IpoptMathOptInterfaceExt.Optimizer) @ IpoptMathOptInterfaceExt ~/.julia/packages/Ipopt/qJgV6/ext/IpoptMathOptInterfaceExt/MOI_wrapper.jl:1341 [7] optimize!(model::IpoptMathOptInterfaceExt.Optimizer) @ IpoptMathOptInterfaceExt ~/.julia/packages/Ipopt/qJgV6/ext/IpoptMathOptInterfaceExt/MOI_wrapper.jl:1356 [8] optimize! @ ~/.julia/packages/MathOptInterface/9lCV5/src/Bridges/bridge_optimizer.jl:367 [inlined] [9] optimize! @ ~/.julia/packages/MathOptInterface/9lCV5/src/MathOptInterface.jl:122 [inlined] [10] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{…}) @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/9lCV5/src/Utilities/cachingoptimizer.jl:370 [11] optimize!(model::JuMP.Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{}) @ JuMP ~/.julia/packages/JuMP/e83v9/src/optimizer_interface.jl:609 [12] optimize! @ ~/.julia/packages/JuMP/e83v9/src/optimizer_interface.jl:560 [inlined]

[mlubin]

MWE, or just any standalone example, would be greatly appreciated.

[odow]

Hi @basberr I'll take a look.

Can you share the model?

using JuMP
write_to_file(model, "model.mof.json")

[odow]

I've been mucking around with this, and I can't get the assert to trigger. Nothing seems obviously wrong reading through the code, so I guess I need an example.

[basberr]

Here is the model.
On my side, this works on v1.46 but not v1.47 where I get the assertion error.

model.mof.json

[odow]

Thanks! I can reproduce:

julia> using Revise

julia> using JuMP, Ipopt

julia> model = read_from_file("/Users/odow/Downloads/model.mof.json"; use_nlp_block = false)
A JuMP Model
├ solver: none
├ objective_sense: MIN_SENSE
│ └ objective_function_type: NonlinearExpr
├ num_variables: 195
├ num_constraints: 312
│ ├ NonlinearExpr in MOI.EqualTo{Float64}: 100
│ ├ AffExpr in MOI.EqualTo{Float64}: 11
│ ├ QuadExpr in MOI.EqualTo{Float64}: 50
│ ├ VariableRef in MOI.GreaterThan{Float64}: 88
│ ├ VariableRef in MOI.LessThan{Float64}: 55
│ └ VariableRef in MOI.Parameter{Float64}: 8
└ Names registered in the model: none

julia> set_optimizer(model, Ipopt.Optimizer)

julia> optimize!(model)
ERROR: AssertionError: k == rinfo.nnz + N
Stacktrace:
  [1] _indirect_recover_structure(rinfo::MathOptInterface.Nonlinear.ReverseAD.Coloring.RecoveryInfo)
    @ MathOptInterface.Nonlinear.ReverseAD.Coloring ~/git/jump-dev/MathOptInterface/src/Nonlinear/ReverseAD/Coloring/Coloring.jl:447

Working on this now.

[odow]

julia> using JuMP, Ipopt

julia> begin
           model = Model(Ipopt.Optimizer)
           @variable(model, x[1:7])
           @objective(
               model,
               Min,
               x[1]*x[5] + x[2]*x[4] + x[2]*x[6] + x[3]*x[5] + (x[1]*x[4]+x[3]*x[6])*x[7]
           )
           optimize!(model)
       end
ERROR: AssertionError: k == rinfo.nnz + N

[mlubin]

Is the sparsity detection correct? Are the parentheses in the last term necessary?

[odow]

Simpler and pure MOI. All bilinear terms appear necessary.

julia> import MathOptInterface as MOI

julia> begin
           nlp = MOI.Nonlinear.Model()
           x = MOI.VariableIndex.(1:6)
           MOI.Nonlinear.set_objective(
               nlp,
               Expr(
                   :call,
                   :+,
                   :($(x[1]) * $(x[4])),
                   :($(x[1]) * $(x[5])),
                   :($(x[2]) * $(x[4])),
                   :($(x[2]) * $(x[6])),
                   :($(x[3]) * $(x[5])),
                   :($(x[3]) * $(x[6])),
               )
           )
           backend = MOI.Nonlinear.SparseReverseMode()
           evaluator = MOI.Nonlinear.Evaluator(nlp, backend, x)
           MOI.initialize(evaluator, [:Hess])
       end
ERROR: AssertionError: k == rinfo.nnz + N
Stacktrace:
 [1] _indirect_recover_structure(rinfo::MathOptInterface.Nonlinear.ReverseAD.Coloring.RecoveryInfo)
   @ MathOptInterface.Nonlinear.ReverseAD.Coloring ~/git/jump-dev/MathOptInterface/src/Nonlinear/ReverseAD/Coloring/Coloring.jl:447
 [2] hessian_color_preprocess(edgelist::Set{…}, num_total_var::Int64, seen_idx::MathOptInterface.Nonlinear.ReverseAD.Coloring.IndexedSet)
   @ MathOptInterface.Nonlinear.ReverseAD.Coloring ~/git/jump-dev/MathOptInterface/src/Nonlinear/ReverseAD/Coloring/Coloring.jl:489
 [3] MathOptInterface.Nonlinear.ReverseAD._FunctionStorage(nodes::Vector{…}, const_values::Vector{…}, num_variables::Int64, coloring_storage::MathOptInterface.Nonlinear.ReverseAD.Coloring.IndexedSet, want_hess::Bool, subexpressions::Vector{…}, dependent_subexpressions::Vector{…}, subexpression_linearity::Vector{…}, subexpression_edgelist::Vector{…}, subexpression_variables::Vector{…}, moi_index_to_consecutive_index::Dict{…})
   @ MathOptInterface.Nonlinear.ReverseAD ~/git/jump-dev/MathOptInterface/src/Nonlinear/ReverseAD/types.jl:97
 [4] initialize(d::MathOptInterface.Nonlinear.ReverseAD.NLPEvaluator, requested_features::Vector{Symbol})
   @ MathOptInterface.Nonlinear.ReverseAD ~/git/jump-dev/MathOptInterface/src/Nonlinear/ReverseAD/mathoptinterface_api.jl:104
 [5] initialize(evaluator::MathOptInterface.Nonlinear.Evaluator{…}, features::Vector{…})
   @ MathOptInterface.Nonlinear ~/git/jump-dev/MathOptInterface/src/Nonlinear/evaluator.jl:105
 [6] top-level scope
   @ REPL[283]:19
Some type information was truncated. Use `show(err)` to see complete types.

[mlubin]

Does the error go away if you add diagonal x[i]^2 terms?

[odow]

Nope. Sparsity detection is correct. It's something in the colouring.

julia> begin
           nlp = MOI.Nonlinear.Model()
           x = MOI.VariableIndex.(1:6)
           MOI.Nonlinear.set_objective(
               nlp,
               Expr(
                   :call,
                   :+,
                   :($(x[1]) * $(x[4])),
                   :($(x[1]) * $(x[5])),
                   :($(x[2]) * $(x[4])),
                   :($(x[2]) * $(x[6])),
                   :($(x[3]) * $(x[5])),
                   :($(x[3]) * $(x[6])),
               )
           )
           backend = MOI.Nonlinear.SparseReverseMode()
           evaluator = MOI.Nonlinear.Evaluator(nlp, backend, x)
           MOI.initialize(evaluator, [:Hess])
       end
edgelist = Set([(4, 2), (5, 1), (4, 1), (6, 2), (5, 3), (6, 3)])
ERROR: AssertionError: k == rinfo.nnz + N
Stacktrace:
 [1] _indirect_recover_structure(rinfo::MathOptInterface.Nonlinear.ReverseAD.Coloring.RecoveryInfo)
   @ MathOptInterface.Nonlinear.ReverseAD.Coloring ~/git/jump-dev/MathOptInterface/src/Nonlinear/ReverseAD/Coloring/Coloring.jl:447

[mlubin]

Try force-adding the diagonal entries to the sparsity pattern.

[odow]

julia> using MathOptInterface.Nonlinear.ReverseAD: Coloring

julia> begin
           I = [4, 5, 4, 6, 5, 6]
           J = [2, 1, 1, 2, 3, 3]
           g = Coloring.UndirectedGraph(I, J, length(I))
           color, num_colors = Coloring.acyclic_coloring(g)
           rinfo = Coloring.recovery_preprocess(g, color, num_colors, 1:6)
           I, J = Coloring._indirect_recover_structure(rinfo)
       end
ERROR: AssertionError: k == rinfo.nnz + N
Stacktrace:
 [1] _indirect_recover_structure(rinfo::MathOptInterface.Nonlinear.ReverseAD.Coloring.RecoveryInfo)
   @ MathOptInterface.Nonlinear.ReverseAD.Coloring ~/git/jump-dev/MathOptInterface/src/Nonlinear/ReverseAD/Coloring/Coloring.jl:447
 [2] top-level scope
   @ REPL[33]:7

Diagonals don't help:

julia> begin
           I = [1, 2, 3, 4, 5, 6, 4, 5, 4, 6, 5, 6]
           J = [1, 2, 3, 4, 5, 6, 2, 1, 1, 2, 3, 3]
           g = Coloring.UndirectedGraph(I, J, length(I))
           color, num_colors = Coloring.acyclic_coloring(g)
           rinfo = Coloring.recovery_preprocess(g, color, num_colors, 1:6)
           I, J = Coloring._indirect_recover_structure(rinfo)
       end
ERROR: AssertionError: k == rinfo.nnz + N
Stacktrace:
 [1] _indirect_recover_structure(rinfo::MathOptInterface.Nonlinear.ReverseAD.Coloring.RecoveryInfo)
   @ MathOptInterface.Nonlinear.ReverseAD.Coloring ~/git/jump-dev/MathOptInterface/src/Nonlinear/ReverseAD/Coloring/Coloring.jl:447
 [2] top-level scope
   @ REPL[34]:7

[odow]

I need to have some lunch. Will be back in a bit

[odow]

So the colouring changed.

julia> using MathOptInterface.Nonlinear.ReverseAD: Coloring

julia> I = [4, 5, 4, 6, 5, 6]
6-element Vector{Int64}:
 4
 5
 4
 6
 5
 6

julia> J = [2, 1, 1, 2, 3, 3]
6-element Vector{Int64}:
 2
 1
 1
 2
 3
 3

julia> g = Coloring.UndirectedGraph(I, J, length(I))
MathOptInterface.Nonlinear.ReverseAD.Coloring.UndirectedGraph([5, 4, 4, 6, 5, 6, 2, 1, 1, 3, 2, 3], [2, 3, 1, 4, 5, 6, 1, 3, 2, 5, 4, 6], [1, 3, 5, 7, 9, 11, 13], [(4, 2), (5, 1), (4, 1), (6, 2), (5, 3), (6, 3)])

julia> color, num_colors = Coloring.acyclic_coloring(g)
([1, 1, 1, 2, 2, 2], 2)

julia> rinfo = Coloring.recovery_preprocess(g, color, num_colors, 1:6)
MathOptInterface.Nonlinear.ReverseAD.Coloring.RecoveryInfo([[4, 2, 5, 1, 6, 3]], [[4, 3, 6, 5, 2, 1]], [[0, 1, 6, 3, 2, 5]], [1, 1, 1, 2, 2, 2], 2, 6, [1, 2, 3, 4, 5, 6])

julia> I, J = Coloring._indirect_recover_structure(rinfo)
ERROR: AssertionError: k == rinfo.nnz + N

where previously

(moi) pkg> st
Status `/private/tmp/moi/Project.toml`
⌃ [b8f27783] MathOptInterface v1.46.0
Info Packages marked with ⌃ have new versions available and may be upgradable.

julia> using MathOptInterface.Nonlinear.ReverseAD: Coloring

julia> I = [4, 5, 4, 6, 5, 6]
6-element Vector{Int64}:
 4
 5
 4
 6
 5
 6

julia> J = [2, 1, 1, 2, 3, 3]
6-element Vector{Int64}:
 2
 1
 1
 2
 3
 3

julia> g = Coloring.UndirectedGraph(I, J, length(I))
MathOptInterface.Nonlinear.ReverseAD.Coloring.UndirectedGraph([5, 4, 4, 6, 5, 6, 2, 1, 1, 3, 2, 3], [2, 3, 1, 4, 5, 6, 1, 3, 2, 5, 4, 6], [1, 3, 5, 7, 9, 11, 13], [(4, 2), (5, 1), (4, 1), (6, 2), (5, 3), (6, 3)])

julia> color, num_colors = Coloring.acyclic_coloring(g)
([1, 1, 1, 2, 2, 3], 3)

julia> rinfo = Coloring.recovery_preprocess(g, color, num_colors, 1:6)
MathOptInterface.Nonlinear.ReverseAD.Coloring.RecoveryInfo([[4, 2, 5, 1, 3], [6, 2, 3]], [[2, 5, 3, 4, 1], [2, 3, 1]], [[0, 1, 4, 1, 3], [0, 1, 1]], [1, 1, 1, 2, 2, 3], 3, 6, [1, 2, 3, 4, 5, 6])

julia> I, J = Coloring._indirect_recover_structure(rinfo)
([1, 2, 3, 4, 5, 6, 4, 5, 5, 4, 6, 6], [1, 2, 3, 4, 5, 6, 2, 3, 1, 1, 2, 3])

[odow]

The culprit might actually be #2885