Broadcasting variables constrained on creation

[odow]

We have a decision to make about what to do if a user writes the following:

model = Model()
@variable(model, x[1:2] in MOI.Integer())

Should we

• (a) throw a nice error message because MOI.Integer is a scalar-set and x is a vector; or
• (b) broadcast the scalar set over the vector.

Either case needs a method like:

function build_variable(
    _error::Function, 
    variables::Vector{<:ScalarVariable},
    set::MOI.AbstractScalarSet
)
    return VariablesConstrainedOnCreation.(variables, Ref(set))
end

along with some tests. (If (a), the return replaced by a nice error message.)

cc @Wikunia

[Wikunia]

I think if there are sets where broadcasting is not reasonable we should have at least a way to allow broadcasting for when it is useful i.e by having an extra set type like AbstractScalarSetAllowBroadcast <: AbstractScalarSet and then have this build_variable function for AbstractScalarSetAllowBroadcast . I would assume that in most cases broadcasting is natural and expected.

[blegat]

I would vote for (b) as it would also resolve a common feature requested from user coming from cvx. In cvx, you can create multiple symmetric matrices or multiples psd matrices at once (see http://cvxr.com/cvx/doc/basics.html):

Matrix-specific keywords can be applied to n-dimensional arrays as well: each 2-dimensional “slice” of the array is given the stated structure. So for instance, the declaration variable R(10,10,8) hermitian semidefinite constructs 8 10×10 complex Hermitian PSD matrices, stored in the 2-D slices of R.

where here it currently fails:

julia> @variable(model, x[1:10, 1:10, 1:8] in PSDCone())
ERROR: MethodError: no method matching build_variable(::JuMP.var"#_error#83"{Tuple{Symbol,Expr}}, ::Array{ScalarVariable{Float64,Float64,Float64,Float64},3}, ::PSDCone)
Closest candidates are:
  build_variable(::Function, ::Array{#s92,2} where #s92<:ScalarVariable, ::PSDCone) at /home/blegat/.julia/dev/JuMP/src/sd.jl:182
  build_variable(::Function, ::Array{#s92,1} where #s92<:ScalarVariable, ::AbstractVectorSet) at /home/blegat/.julia/dev/JuMP/src/sets.jl:6
  build_variable(::Function, ::Array{#s92,2} where #s92<:ScalarVariable, ::SymMatrixSpace) at /home/blegat/.julia/dev/JuMP/src/sd.jl:164
  ...
Stacktrace:
 [1] top-level scope at /home/blegat/.julia/dev/JuMP/src/macros.jl:45

The reasoning would be similar to broadcast in Base. Since PSD variables are matrices and the user want to create an Array{T,3}, we consider the PSD matrix dimensions to be along the 2 first of the array and broadcast along the third one. So here, x[:, :, 1] would be PSD, x[:, :, 2] would be PSD, ... This is consistent with what cvx does.

I would assume that in most cases broadcasting is natural and expected.

We could add a trait in case some sets want to throw a warning or error if it gets broadcasted but we can wait to see if there is a set that needs that.

[odow]

Urgh, the Array{3} business is really non-obvious. What happens in N-dimensions? Am I reading that comment right that

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

will create 9 2x2 PSD matrices? So x[:, :, i j] is PSD?

Does

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

create 4 SOC vectors?

I guess it is a natural extension from broadcasting scalar sets over a vector. We'd just have to make it consistent for all sets and dimensions.

[blegat]

@odow Yes, exactly, it's a bit odd but if we make it consistent with each set and dimension and with Julia Base broadcasting, it should not be too confusing once the user gets used to it.

[Wikunia]

I used the example of integer sets like @variable(m, x in IntegerSet([1,3,5])) which had this problem. How about using @variable(m, x, Integers([1,3,5]) instead? This works after adding:

function JuMP.build_variable(_error::Function, info::JuMP.VariableInfo, setType::Type{T}) where T<:MOI.AbstractScalarSet
    return JuMP.VariableConstrainedOnCreation(JuMP.ScalarVariable(info), setType())
end

function JuMP.build_variable(_error::Function, info::JuMP.VariableInfo, set::T) where T<:MOI.AbstractScalarSet
    return JuMP.VariableConstrainedOnCreation(JuMP.ScalarVariable(info), set)
end

The one above is only necessary if something like @variable(m, x, Odd) is used.

Here broadcasting is already implemented the same way as it is for Bin and Int.

[blegat]

Yes, this was kind of the old way to specify constrained variables. This is what is used in SumOfSquares and PolyJuMP where you can create containers of polynomials and broadcast works. The in syntax is more consistent with @constraint and more elegant but currently the use cases of these two syntax indeed overlap

[odow]

Copied from: https://github.com/jump-dev/JuMP.jl/pull/2416#discussion_r543767698, which has a suggestion that makes the syntax a bit clearer:

We don't support

@variable(model, x[2:3, 2:3] in SkewSymMatrixSpace())

because it isn't apparent you can have a symmetric DenseAxisArray.

But if users are forced to write

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

why not just

@variable(model, X in SkewSymMatrixSpace(3))

For anonymous variables, we could support

X = @variable(model, _ in SkewSymMatrixSpace(3))

Then

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

is a 3x3 matrix of 2x2 skew-symmetric matrices.

It should be possible to have a deprecated way of doing this, because we can distinguish between SkewSymMatrixSpace() and SkewSymMatrixSpace(2).

[odow]

Related issue: https://github.com/jump-dev/JuMP.jl/issues/2056

[odow]

#3184 resolves another array-in-scalar-set issue. But I'm not convinced about changing the vector-valued sets. It seems ripe to add more confusion, rather than less.

People can easily do things like

X = [
    @variable(model, [1:3, 1:3] in SkewSymMatrixSpace())
    for i in 1:2, j in 1:2
]

[odow]

@blegat are you okay with closing this?

it should not be too confusing once the user gets used to it.

I'm not convinced that the syntax:

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

is less confusing compared with the currently working:

x = [@variable(model, [1:2, 1:2] in PSDCone()) for i in 1:3, j in 1:3]
x = [@variable(model, [1:3] in SecondOrderCone()) for i in 1:4]

[blegat]

Yes, for vector set we could wait to see if there is a real need for it. But we could at least support it for abitrary scalar set. https://github.com/jump-dev/JuMP.jl/pull/3184/ is only for ComplexPlane so it does not work for abitrary scalar set IIUC

[odow]

But we could at least support it for abitrary scalar set.

We already supported scalar sets: https://github.com/jump-dev/JuMP.jl/pull/2657. A special method was needed for ComplexPlane because it isn't a MOI.AbstractSet.

Closing because this seems resolved. If any future reader has a strong wish for the vector syntax @variable(model, x[1:3, 1:4] in SecondOrderCone()) please comment and we can re-open this issue.

[blegat]

If any future reader has a strong wish for the vector syntax @variable(model, x[1:3, 1:4] in SecondOrderCone()) please comment and we can re-open this issue.

For the record, @araujoms just requested exactly that in https://discourse.julialang.org/t/error-when-constraining-variable-to-be-hermitian-psd/95223

there a more convenient way to declare a set of Hermitian variables? With YALMIP I can just do

E = sdpvar(d,d,n-1,'hermitian','complex');

which is much more concise and readable.