Inequality constraints over symmetric matrices are cumbersome

[LebedevRI]

Unless i'm missing something very obvious, there is no nice way to exploit the symmetry of matrices (to avoid adding duplicate constraints) when using <= and >= constraints. Here's what i've tried:

using JuMP
import LinearAlgebra
import MathOptInterface as MOI

function main()
    model = Model() 

    @variable(model, v[1:3,1:3], Symmetric)

    # Results in nice triangular matrix in Zeros()
    @constraint(model, c0, LinearAlgebra.Symmetric(v.-1)  == LinearAlgebra.Symmetric(fill(41, 3,3))) # Yay!

    # Unsupported matrix in vector-valued set.
    # @constraint(model, c1, LinearAlgebra.Symmetric(v.-1)  >= LinearAlgebra.Symmetric(fill(41, 3,3))) # No.

    # Works, but results in 9 constraints, and all 9 elements are used, instead of always using upper-triangular elements.
    @constraint(model, c2, LinearAlgebra.Symmetric(v.-1) .>= LinearAlgebra.Symmetric(fill(41, 3,3))) # Yes, but no symmetry.

    # Unsupported matrix in vector-valued set.
    # @constraint(model, c3, LinearAlgebra.Symmetric(v.-42) >= LinearAlgebra.Symmetric(fill(0, 3,3))) # No.

    # Unsupported matrix in vector-valued set.
    #@constraint(model, c4, LinearAlgebra.Symmetric(v.-42) >= 0) # No

    # Works, but still no notion of symmetry.
    @constraint(model, c5, vec(LinearAlgebra.Symmetric(v.-42)) >= 0)

    # MethodError: no method matching triangular_form
    # @constraint(model, c6, vec(MOI.triangular_form(LinearAlgebra.Symmetric(v.-42))) in MOI.Nonnegatives(3*3)) # No

    # RHS is not lower-triangular, you get `-42 in Nonnegatives` elts. Oops.
    @constraint(model, c7, vec(LinearAlgebra.LowerTriangular(v.-1)) >= vec(fill(41, 3,3))) # Yes, kind of.

    # Finally, this works, no duplicate constraints, but the vector is still 9-element (padded with zeros)
    @constraint(model, c8, vec(LinearAlgebra.LowerTriangular(v.-1)) >= vec(LinearAlgebra.LowerTriangular(fill(41, 3,3)))) # FINALLY!

    print(model)
end

main()

It seems to be it's a bit too brittle, no? I would have guessed there should be Nonnegatives/Nonpositives sets, similar to the existing Zeros argument-less set, and such comparisons between symmetric matrices should work similarly with equality-comparison.

[blegat]

We can probably make this work

@constraint(model, Symmetric(A) == Symmetric(B))

by adding something like

function build_constraint(error_fn, Q::Symmetric, set::Union{Zeros,Nonnegatives,Nonpositives})
    n = LinearAlgebra.checksquare(Q)
    shape = SymmetricMatrixShape(n)
    return VectorConstraint(
        vectorize(Q, shape),
        set,
        shape,
    )
end

For the inequality, I'd prefer an explicit:

@constraint(model, Symmetric(A) >= Symmetric(B), Nonnegatives())

because most user would actually expect a PSD constraint, not elementwise inequality when you do an inequality between matrices so it would be too confusing to make it work without the extra argument. Currently, Nonnegatives will be added by default and you get an error: https://github.com/jump-dev/JuMP.jl/blob/56b186fc18d854be33a7f558bc2ba1918dd7c86d/src/sd.jl#L575-L582 Adding the build_constraint I suggest above wouldn't distinguish between the case where Nonnegatives is added explicitly. We would need a way to distinguish it before allow Nonnegatives in build_constraint.

[LebedevRI]

We can probably make this work

@constraint(model, Symmetric(A) == Symmetric(B))

Err, note that for equality, it already works.

For the inequality, I'd prefer an explicit:

@constraint(model, Symmetric(A) >= Symmetric(B), Nonnegatives())

Without ever having used the PSD side of things, this syntax looks kinda weird to me. What happens if the wrong set is specified?

[odow]

because most user would actually expect a PSD constraint

Agreed. We discussed this, see https://github.com/jump-dev/JuMP.jl/pull/3281#discussion_r1133340687

@constraint(model, Symmetric(A) >= Symmetric(B), Nonnegatives())

:+1: this solves a lot of issues.

What happens if the wrong set is specified?

There is no such thing as the "wrong" set. The constraint is rewritten to Symmetric(A) - Symmetric(B) in Set()

[LebedevRI]

What happens if the wrong set is specified?

There is no such thing as the "wrong" set. The constraint is rewritten to Symmetric(A) - Symmetric(B) in Set()

What i mean is, there's both the inequality comparison between A and B, and a statement about the set to which the A-B belong to. One of these is redundant.

A >= B, is A - B >= 0 aka A - B in Nonnegatives, sure, absolutely.

But what would A - B >= 0, Nonpositives aka A - B in Nonnegatives in Nonpositives mean? That's a bogus constraint (it should be an equality-with-zero constraint), that is not exactly easy to spot, so i would very much hope that it would be invalid and rejected with a diagnostic. And by that point, what's the point of accepting both the set and the inequality operator? It's a footgun. (At most, perhaps it should be some syntactic sugar placeholder set that is otherwise unused.)

[odow]

See https://jump.dev/JuMP.jl/stable/manual/constraints/#Set-inequality-syntax

The set defines the partial ordering of >=.

@constraint(model, X >= Y, Nonpositives() means X - Y in Nonpositives(), or X - Y <= 0.

[LebedevRI]

Oh wow, that is the most exquisite design decision (strictly from the user-facing view, i'm sure it makes sense from the internal implementation detail) i've seen in last few years...

[odow]

I might be misinterpreting your comment, but this is not a design decision that is unique to JuMP.

It follows from the syntax for the generalized inequality of a cone:

$x \succeq_K y \iff x - y \in K$

See, e.g., Section 2.4.1 of Boyd and Vandenberghe.

[LebedevRI]

I might be misinterpreting your comment,

It is more that i'm failing to fully convey my point :)

It follows from the syntax for the generalized inequality of a cone:

x⪰Ky⟺x−y∈K

See, e.g., Section 2.4.1 of Boyd and Vandenberghe.

Yes, but that's ⪰, which is not the same symbol as ≥, though they may look similar in small font size. I just feel like that ≥ have a very well defined meaning, so reusing them in this context is confusing.

Anyways, that would be better than nothing i suppose :)

[odow]

but that's ⪰, which is not the same symbol as ≥,

Sure. We intentionally conflate the two because ⪰ is hard to type in ASCII.