Add set inequality syntax for matrices (#3766)

docs/src/manual/constraints.md

@@ -164,27 +164,46 @@ julia> @constraint(model, A * x .>= b)
  3 x[1] + 4 x[2] ≥ 6
 ```
 
-### Vectorized matrix constraints
+## Matrix inequalities
 
-In most cases, you cannot use the non-broadcasting syntax for general matrices.
-For example:
+Inequalities between matrices are not supported, due to the common ambiguity
+between elementwise inequalities and a [`PSDCone`](@ref) constraint.
 
-```jldoctest
+```jldoctest symmetric_matrix
 julia> model = Model();
 
-julia> @variable(model, X[1:2, 1:2])
-2×2 Matrix{VariableRef}:
- X[1,1]  X[1,2]
- X[2,1]  X[2,2]
+julia> @variable(model, x[1:2, 1:2], Symmetric);
+
+julia> @variable(model, y[1:2, 1:2], Symmetric);
 
-julia> @constraint(model, X >= 0)
-ERROR: At none:1: `@constraint(model, X >= 0)`: Unsupported matrix in vector-valued set. Did you mean to use the broadcasting syntax `.>=` instead? Alternatively, perhaps you are missing a set argument like `@constraint(model, X >= 0, PSDCone())` or `@constraint(model, X >= 0, HermitianPSDCone())`.
+julia> @constraint(model, x >= y)
+ERROR: At none:1: `@constraint(model, x >= y)`: Unsupported matrix in vector-valued set. Did you mean to use the broadcasting syntax `.>=` instead? Alternatively, perhaps you are missing a set argument like `@constraint(model, X >= 0, PSDCone())` or `@constraint(model, X >= 0, HermitianPSDCone())`.
 Stacktrace:
 [...]
 ```
 
-Instead, to represent matrix inequalities you must always use the element-wise
-broadcasting `.==`, `.>=`, or `.<=`, or use the [Set inequality syntax](@ref).
+Instead, use the [Set inequality syntax](@ref) to specify a set like
+[`PSDCone`](@ref) or [`Nonnegatives`](@ref):
+
+```jldoctest symmetric_matrix
+julia> @constraint(model, x >= y, PSDCone())
+[x[1,1] - y[1,1]  x[1,2] - y[1,2]
+ ⋯                x[2,2] - y[2,2]] ∈ PSDCone()
+
+julia> @constraint(model, x >= y, Nonnegatives())
+[x[1,1] - y[1,1]  x[1,2] - y[1,2]
+ ⋯                x[2,2] - y[2,2]] ∈ Nonnegatives()
+
+julia> @constraint(model, x >= y, Nonpositives())
+[x[1,1] - y[1,1]  x[1,2] - y[1,2]
+ ⋯                x[2,2] - y[2,2]] ∈ Nonpositives()
+
+julia> @constraint(model, x >= y, Zeros())
+[x[1,1] - y[1,1]  x[1,2] - y[1,2]
+ ⋯                x[2,2] - y[2,2]] ∈ Zeros()
+```
+
+### Special cases
 
 There are two exceptions: if the result of the left-hand side minus the
 right-hand side is a `LinearAlgebra.Symmetric` matrix or a `LinearAlgebra.Hermitian`

docs/styles/config/vocabularies/JuMP/accept.txt

@@ -19,6 +19,7 @@ discretize(d|s)
 docstring(?s)
 doctest(?s)
 embeddable
+elementwise
 [Ee]num(?s)
 errored
 flamegraph

src/sd.jl

@@ -570,8 +570,6 @@ function build_constraint(
     return VectorConstraint(x, MOI.Zeros(length(x)), shape)
 end
 
-reshape_set(::MOI.Zeros, ::SymmetricMatrixShape) = Zeros()
-
 function build_constraint(error_fn::Function, ::AbstractMatrix, ::Nonnegatives)
     return error_fn(
         "Unsupported matrix in vector-valued set. Did you mean to use the " *
@@ -657,3 +655,41 @@ function build_variable(
     x = _vectorize_variables(error_fn, variables)
     return VariablesConstrainedOnCreation(x, set, SymmetricMatrixShape(n))
 end
+
+moi_set(::Nonnegatives, dim::Int) = MOI.Nonnegatives(dim)
+moi_set(::Nonpositives, dim::Int) = MOI.Nonpositives(dim)
+moi_set(::Zeros, dim::Int) = MOI.Zeros(dim)
+
+shape(f::LinearAlgebra.Symmetric) = SymmetricMatrixShape(size(f, 1))
+
+reshape_set(::MOI.Nonnegatives, ::SymmetricMatrixShape) = Nonnegatives()
+reshape_set(::MOI.Nonpositives, ::SymmetricMatrixShape) = Nonpositives()
+reshape_set(::MOI.Zeros, ::SymmetricMatrixShape) = Zeros()
+
+shape(f::Array) = ArrayShape(size(f))
+
+reshape_set(::MOI.Nonnegatives, ::ArrayShape) = Nonnegatives()
+reshape_set(::MOI.Nonpositives, ::ArrayShape) = Nonpositives()
+reshape_set(::MOI.Zeros, ::ArrayShape) = Zeros()
+
+function build_constraint(
+    error_fn::Function,
+    f::Union{Array,LinearAlgebra.Symmetric},
+    ::Nonnegatives,
+    set::Union{Nonnegatives,Nonpositives,Zeros},
+)
+    s = shape(f)
+    x = vectorize(f, s)
+    return VectorConstraint(x, moi_set(set, length(x)), s)
+end
+
+function build_constraint(
+    error_fn::Function,
+    ::Union{Array,LinearAlgebra.Symmetric},
+    ::Nonpositives,
+    set::Union{Nonnegatives,Nonpositives,Zeros},
+)
+    return error_fn(
+        "The syntax `x <= y, $set` not supported. Use `y >= x, $set` instead.",
+    )
+end

src/shapes.jl

@@ -168,3 +168,31 @@ struct VectorShape <: AbstractShape end
 reshape_vector(vectorized_form, ::VectorShape) = vectorized_form
 
 vectorize(x, ::VectorShape) = x
+
+"""
+    ArrayShape{N}(dims::NTuple{N,Int}) where {N}
+
+An [`AbstractShape`](@ref) that represents array-valued constraints.
+
+## Example
+
+```jldoctest
+julia> model = Model();
+
+julia> @variable(model, x[1:2, 1:3]);
+
+julia> c = @constraint(model, x >= 0, Nonnegatives())
+[x[1,1]  x[1,2]  x[1,3]
+ x[2,1]  x[2,2]  x[2,3]] ∈ Nonnegatives()
+
+julia> shape(constraint_object(c))
+ArrayShape{2}((2, 3))
+```
+"""
+struct ArrayShape{N} <: AbstractShape
+    dims::NTuple{N,Int}
+end
+
+reshape_vector(x, shape::ArrayShape) = reshape(x, shape.dims)
+
+vectorize(x, ::ArrayShape) = vec(x)

test/test_constraint.jl

@@ -1887,4 +1887,58 @@ function test_set_normalized_coefficient_quadratic_batch()
     return
 end
 
+function test_symmetric_matrix_inequality()
+    model = Model()
+    @variable(model, x[1:2, 1:2], Symmetric)
+    @variable(model, y[1:2, 1:2], Symmetric)
+    set_start_value.(x, [1 2; 2 3])
+    set_start_value.(y, [6 4; 4 7])
+    g = [x[1, 1] - y[1, 1], x[1, 2] - y[1, 2], x[2, 2] - y[2, 2]]
+    for set in (Nonnegatives(), Nonpositives(), Zeros())
+        c = @constraint(model, x >= y, set)
+        o = constraint_object(c)
+        @test isequal_canonical(o.func, g)
+        @test o.set == moi_set(set, 3)
+        @test o.shape == SymmetricMatrixShape(2)
+        @test reshape_set(o.set, o.shape) == set
+        primal = value(start_value, c)
+        @test primal isa LinearAlgebra.Symmetric
+        @test primal == LinearAlgebra.Symmetric([-5.0 -2.0; -2.0 -4.0])
+        @test_throws_runtime(
+            ErrorException(
+                "In `@constraint(model, x <= y, set)`: The syntax `x <= y, $set` not supported. Use `y >= x, $set` instead.",
+            ),
+            @constraint(model, x <= y, set),
+        )
+    end
+    return
+end
+
+function test_matrix_inequality()
+    model = Model()
+    @variable(model, x[1:2, 1:3])
+    @variable(model, y[1:2, 1:3])
+    set_start_value.(x, [1 2 3; 4 5 6])
+    set_start_value.(y, [7 9 11; 8 12 13])
+    g = vec(x .- y)
+    for set in (Nonnegatives(), Nonpositives(), Zeros())
+        c = @constraint(model, x >= y, set)
+        o = constraint_object(c)
+        @test isequal_canonical(o.func, g)
+        @test o.set == moi_set(set, 6)
+        @test o.shape == ArrayShape((2, 3))
+        @test reshape_set(o.set, o.shape) == set
+        primal = value(start_value, c)
+        @test primal isa Matrix{Float64}
+        @test primal == [-6.0 -7.0 -8.0; -4.0 -7.0 -7.0]
+        @test_throws_runtime(
+            ErrorException(
+                "In `@constraint(model, x <= y, set)`: The syntax `x <= y, $set` not supported. Use `y >= x, $set` instead.",
+            ),
+            @constraint(model, x <= y, set),
+        )
+    end
+    return
+end
+
 end  # module