Setting heuristic solution involving symmetric matrix in JuMP callback

Specific Domains Optimization (Mathematical)
1

Dear All,

I am experimenting with JuMP callback with Gurobi’s MIQCP feature, where I am storing the heuristic values based on tips from https://discourse.julialang.org/t/jump-cplex-how-to-set-heuristic-solution-using-sparseaxisarrays-in-callback/56938. However, I see a Dimension mismatch while setting my variables to heuristic values.

using Gurobi, JuMP

const GRB_ENV = Gurobi.Env() 

model = direct_model(Gurobi.Optimizer(GRB_ENV))

set_optimizer_attribute(model, "NonConvex", 2)
# "NonConvex" => 2 tells Gurobi to use its nonconvex algorithm

# declare the variables
# ----------------------------
@variable(model, x[1:2] >= 0)
@variable(model, y >= 0)
@variable(model, z >=0)
@variable(model, M[1:3,1:3], Symmetric)

# add constraints
# --------------------
@constraint(model, x[1] + x[2] + y + z <= 10)
@constraint(model, x[1]*y + x[2]*z <= 1)
@constraint(model, x[1]*z + x[2]*y == 1)

# add objective
# -----------------
@objective(model, Min, x[1] + x[2])

# this is the callback function to implement the heuristic solution
# -----------------------------------------------------------------------------

function my_heuristic_callback_function(cb_data)
    # load the current values of the variables at a node of the branch-and-bound tree
    x_val = callback_value.(cb_data, x)
    y_val = callback_value.(cb_data, y)
    z_val = callback_value.(cb_data, z)
    heuristic_flag = 0
    # println(α_val)
    if abs(x_val[2]) <= 1e-4 # this is the condition for applying the heuristic
        @info "[🐼 ] Heuristic condition satisfied: Status of the callback node is: $(callback_node_status(cb_data, model))"  # show the status of the callback node
        @info "[🎤 ] The values of the variables at the callback node $([x_val; y_val; z_val])"
        x_heuristic = [0.10102051067754236; 0.0]
        y_heuristic = 0.0
        z_heuristic = 9.898979489322459
        M_heuristic = [0 for i in 1:3, j in 1:3]
        # Load the JuMP variables in a vertical vector pointwise
        JuMP_variables = vcat(
           [model[:x][i] for i in eachindex(model[:x])],
           y,
           z,
           [model[:M][i_j] for i_j in eachindex(model[:M])] 
           # 💀: the line above causes the problem, wrapping it into vec(⋅) seems to fix the issue
           )
        # Load the heuristic solution values in a vertical vector pointwise
        heuristic_values = vcat(
           [x_heuristic[i] for i in eachindex(model[:x])],
           y_heuristic,
           z_heuristic,
           [M_heuristic[i_j] for i_j in eachindex(model[:M])] 
           # 💀: the line above causes the problem, wrapping it into vec(⋅) seems to fix the issue
           )
        # Submit the heuristic solution for potentially improving the current solution
       status = MOI.submit(
             model, MOI.HeuristicSolution(cb_data), JuMP_variables, heuristic_values
             )
      println("[🙀 ] Status of the submitted heuristic solution is: ", status) # The status shows if the submitted heuristic solution is accepted or not
  end
end

# IMPORTANT: This enables the heuristic
MOI.set(model, MOI.HeuristicCallback(), my_heuristic_callback_function)

# time to optimize
# --------------------
optimize!(model)

# store the solution
# ----------------------
x_opt = value.(x)
y_opt = value(y)
z_opt = value(z)
M_opt = value.(M)
p_star = objective_value(model)

Running the code above gives the error: DimensionMismatch("tried to assign 2-element array to 2×3 destination"). The error goes away if I put the problematic lines

[model[:M][i_j] for i_j in eachindex(model[:M])]

and

 [M_heuristic[i_j] for i_j in eachindex(model[:M])]

in vec(.) function . That being said, I am not sure if this vec(.) operation is the right way to store the heuristic values in general, as the JuMP documentation provides somewhat simple examples, so any tips regarding the right practice will be appreciated.

2

Not sure if this solves the problem, but I did not use eachindex() nor i_j. You could try,

[model[:M][i,j] for i in 1:3 for j in 1:3]

and

 [M_heuristic[i,j] for i in 1:3 for j in 1:3]
1 Like
3

Thanks @mike_k ! Yes, this would work for this simple example, but I am not sure if it will extend to other container types in https://jump.dev/JuMP.jl/stable/manual/containers/#Containers, as eachindex(.) seems like the unified way to loop over a generic container.

For example, consider the more complicated case (this particular code snippet runs fine, because I have used the vec operation for Z), where each variable is of a different type:

JuMP_variables = vcat(
     # λ is a dense axis array over a dictionary of custom types
     [BnB_PEP_model[:λ][i_j_λ] for i_j_λ in eachindex(BnB_PEP_model[:λ])],
     # ν is a scalar
     BnB_PEP_model[:ν],
     # Z is a symmetric matrix
     vec([BnB_PEP_model[:Z][i_j] for i_j in eachindex(BnB_PEP_model[:Z])]),
     # α is sparse axis array that ranges over (i = 1:5, j in 1:i-1)
     [BnB_PEP_model[:α][i_j] for i_j in eachindex(BnB_PEP_model[:α]) ]
)

In this particular case, your suggestion will not extend to \lambda and \alpha. I was wondering if there is a unified way to construct these vcat objects over different types of containers?

4

I get the point, however, I do not now how to deal with eachindex(), sorry. In my application I did some index book-keeping (because I needed that anyway quite often). Something like

[z_val[i,j,w] for w in W for i in I for j in Rout[i]]

where W and I are index sets, e.g., W=1:100, and Rout::Vector{Vector{Int}} stores the needed indices. E.g, Rout[1] = [2,5,7]. This way it worked for me.

1 Like
5

Your problem is that

julia> model = Model();

julia> @variable(model, x[1:3, 1:3], Symmetric)
3×3 LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}:
x[1,1] x[1,2] x[1,3]
x[1,2] x[2,2] x[2,3]
x[1,3] x[2,3] x[3,3]

julia> [x[i] for i in eachindex(x)]
3×3 Matrix{VariableRef}:
x[1,1] x[1,2] x[1,3]
x[1,2] x[2,2] x[2,3]
x[1,3] x[2,3] x[3,3]

You could use reshape:

julia> reshape([x[i] for i in eachindex(x)], length(x),)
9-element Vector{VariableRef}:
x[1,1]
x[1,2]
x[1,3]
x[1,2]
x[2,2]
x[2,3]
x[1,3]
x[2,3]
x[3,3]
6

Thanks @odow! Would not the reshape operation be the same as vec that I mentioned? Or would the latter mess up the indices of x or something?

1 Like
7

Actually I think vec(x) = reshape(x, length(x)), so they’re exactly equivalent.

8

Thanks @odow, makes sense !