Constrained variables and variable bridges

[blegat]

Currently, if a variable belongs to a set or cone, it should be entered as:

x = MOI.add_variable(model)
MOI.add_constraint(model, x in set)

or

x = MOI.add_variables(model, MOI.dimension(cone))
MOI.add_constraint(model, x in cone)

This means that the model first need to create free variables and is only informed afterwards that the variables are in fact not free and belong to a set or cone.

This approach has two major drawbacks:

1. For many solvers, this prevents supporting default_copy_to hence it should implement Allocate-Load or a custom copy_to and it prevents the solver from being available in direct mode without a cache. This is because variables are handled differently depending on the cone in which they belong. For instance:
   • In Mosek, SDP variables are treated differently, they are created with a different functions than scalar variables and their coefficients are stored in a different matrix for constraints. Previously, scalar variables were created and then when an SDP constraint was added on them, an SDP variable was created and set equal to the scalar variable but that gave a model very slow to solve (as evidenced by Avoiding creation of scalar variables in semidefinite programs JuMP.jl#1882) and it is currently resolved by a hack that was introduced in Use Matrix Variables MosekTools.jl#22 which is to not create the variable before it is used and if it was used before being set to SDP, we delete the created scalar variables to replace them by SDP variables and move their coefficients form the matrix of scalar variables to the matrices for psd variables. I discussed about this with @erling-d-andersen and eliminating the scalar variables in presolve that would be created in this manner is unlikely to happen soon. IMO relying on presolve instead of generating the correct problem isn't the cleanest solution anyway.
   • In SDPT3, when creating constraints, the coefficients are stored in different matrices depending on their order.
   • In SeDuMi, CDCS, the coefficients in the matrix are stored in a different order depending on the cone in which the variables belong.
2. When creating a VectorOfVariables-in-Cone constraint for which the Cone is not supported, the constraint may be bridged to a VectorAffineFunction-in-SupportedCone which is then bridged to a VectorOfVariables-in-SupportedCone by adding slacks (consider for instance Cone=RotatedSecondOrderCone and SupportedCone=SecondOrderCone or Cone=Nonpositives and SupportedCone=Nonnegatives). The issue here is that twice the number of variables needed have been created. What should have been done is creating variables in SupportedCone and then create dummy MOI variables that do not correspond to variables of the solver (e.g. by giving them negative index as suggested by @odow and @joaquimg during JuMP-dev) and return these dummy variables and every time they are used in a constraint they are replaced by their corresponding expressions in the variables created (this would be done transparently by the LazyBridgeOptimizer). This would be called variable bridges but it is currently impossible to do because the variable creation and the variable constraint creation are done in two distinct MOI calls.
   • In SumOfSquares, a prototype of variable bridges is already implemented https://github.com/JuliaOpt/SumOfSquares.jl/blob/ce6edd78cd4a4ca7f89571f06c8f88964738525c/src/SumOfSquares.jl#L36-L40 to avoid creating unnecessary variables as described above.
   • For CSDP, SDPA, DSDP, SDPLR, only nonnegative and PSD variables are supported. For this reason, many variable bridges are needed. The current solution is to create a layer called SemidefiniteOptInterface for which most of the work is hardcoding these variables bridges (e.g. for a nonpositive variable, we create a nonnegative and we remember to multiply by -1 every time it is used, for a variable GreaterThan(1.0) we create a nonnegative one and remember to shift it every time it is used, for a free variable, we create two nonnegative variables and use the difference, ...). However, it is only hardcoded for a certain number of sets and for other ones, we will have the issue of adding to twice the number variables needed.
   • For SCS, ECOS, SeDuMi, CDCS, the solution is to encode the problem the primal form (for SCS, ECOS) and dual form (for SeDuMi, CDCS) depending on which form has only free variables.
   • For Mosek and SDPT3, the solution used by SCS, ECOS, SeDuMI, CDCS is not an option since Mosek supports MOI.Integer and MOI.Interval and SDPT3 supports log-det in the objective for SDP variables (and log in the objective for nonnegative variables). That means that they suffer from this issue of the creating of unnecessary variables for VectorOfVariables-in-set if set is not natively supported !

The solution we discussed at JuMP-dev is to add a new function

function add_constrained_variables(model, set)
    variables = MOI.add_variables(model, MOI.dimension(set))
    ci = MOI.add_constraint(model, MOI.VectorOfVariables(variables), set)
    return variables, ci
end

and

MOI.supports_variable(::MOI.ModelLike, ::MOI.AbstractSet)

For instance, SemidefiniteOptInterface would implement

MOI.supports_variable(::Optimizer, ::MOI.Reals) = false

which means that MOI.add_variable/MOI.add_variables is not supported; they are equivalent to MOI.add_constrained_variables(model, MOI.Reals(...)) but MOI.Reals is only used in MOI.supports_variable.
Then it would implement

MOI.supports_variable(::Optimizer, ::MOI.Nonnegatives) = true
MOI.supports_variable(::Optimizer, ::MOI.PositiveSemidefinieConeTriangle) = true

and implement add_constrained_variables for these two sets, the rest would be handled by variable bridges!

In JuMP, we could create the new syntax

@variable(model, x[1:3, 1:3] in PSDCone())

instead of

@variable(model, x[1:3, 1:3], PSD)

and it would generalize for any cone.
For instance

@variable(model, x[1:3] in SecondOrderCone())

would translate to a MOI.add_constraint_variables(backend(model), MOI.SecondOrderCone(3)).
If the user does instead

@variable(model, x[1:3])
@constraint(model, x in SecondOrderCone())

then it would be equivalent if there is a CachingOptimizer as the CachingOptimizer will call add_constraint_variables in copy_to but it won't work in direct mode. So in direct mode, only @variable(model, x[1:3], SecondOrderCone()) will work if the optimizer implements add_constraint_variables and not the VectorOfVariables-in-MOI.SecondOrderCone constraint, but that seems faithful with the philosophy of direct mode.

Note that creating constrained variables is similar to the feature of SumofSquares that allows to directly create an SOS polynomial variables (instead of creating a free polynomial and then constrain it to be SOS which would create unnecessary variables). That difference with other SOS modeling framework such as YALMIP allows to create smaller problem that can be solved faster as @mforets showed in JuMP-dev (see https://www.youtube.com/watch?v=4mf2W4y7PJo).