Add set for low-rank constrained SDP

[blegat]

Closes #2197

Note that once this is done, you can write an SDP in SDPA format directly in JuMP with something like

C = rand(d, d)
A = [rand(d, d) for i in 1:m]
model = Model()
@variable(model, x[1:m+1] in MOI.SetWithDotProducts(
    MOI.PositiveSemidefiniteConeTriangle(d),
    MOI.TriangleVectorization.([A; [C]]),
))
@constraint(model, x[1:m] == b)
@objective(model, Min, x[end])

or

model = Model()
@variable(model, y[1:m])
@constraint(model, [y; -1] in MOI.LinearCombinationInSet([A; [C]]))
@objective(model, Max, dot(y, b))

And Dualization.dualize should take you from one to the other one!

Solvers implementing this:

• Hypatia.LinMatrixIneqCone : A1 PSD and A2, ..., A[m] arbitrary
• Hypatia.WSOSInterpNonnegativeCone : A = u u' (Rank-1 PSD)
• DSDP : A = a * u * u' (Rank-1) Add support for rank-1 constraints DSDP.jl#37
• SDPLR : A = U * Diagonal(d) * U' (Low-Rank) Add support for low rank constraint SDPLR.jl#26

Basic

• Add a new AbstractScalarSet or AbstractVectorSet to src/sets.jl
• If isbitstype(S) == false, implement Base.copy(set::S)
• If isbitstype(S) == false, implement Base.:(==)(x::S, y::S)
• If an AbstractVectorSet, implement dimension(set::S), unless the
  dimension is given by set.dimension.

Utilities

• If an AbstractVectorSet, implement Utilities.set_dot,
  unless the dot product between two vectors in the set is equivalent to
  LinearAlgebra.dot
• If an AbstractVectorSet, implement Utilities.set_with_dimension in
  src/Utilities/matrix_of_constraints.jl
• Add the set to the @model macro at the bottom of src/Utilities.model.jl

Documentation

• Add a docstring, which gives the mathematical definition of the set,
  along with an ## Example block containing a jldoctest
• Add the docstring to docs/src/reference/standard_form.md
• Add the set to the relevant table in docs/src/manual/standard_form.md

Tests

• Define a new _set(::Type{S}) method in src/Test/test_basic_constraint.jl
  and add the name of the set to the list at the bottom of that files
• If the set has any checks in its constructor, add tests to test/sets.jl

MathOptFormat

• Open an issue at https://github.com/jump-dev/MathOptFormat to add
  support for the new set {{ replace with link to the issue }} I don't think we should add any set specialization of that set or at least not yet

Optional

• Implement dual_set(::S) and dual_set_type(::Type{S})
• Add new tests to the Test submodule exercising your new set
• Add new bridges to convert your set into more commonly used sets

[odow]

This is getting a little out of hand. PositiveSemidefiniteConeSquare, PositiveSemidefiniteConeTriangle, ScaledPositiveSemidefiniteConeTriangle, HermitianPositiveSemidefiniteConeTriangle...

Do we really need yet another PSD cone? Why isn't this just a solver enhancement where they detect and exploit low-rank structures in the general matrix? Because it's a constraint?

[blegat]

You're right, maybe this isn't really a PSD think, we could generalize it for any cone.

model = Model()
@variable(model, (AX[1:m], X[1:n, 1:n]) in ConeWithInnerProducts(A, PSDCone()))
@constraint(model, AX == b)
@objective(model, Min, dot(X, C))

and

model = Model()
@variable(model, y[1:m])
@constraint(model, (y, C) in AffineSpanInCone(A, PSDCone()))
@objective(model, Max, dot(y, b))

Why isn't this just a solver enhancement where they detect and exploit low-rank structures in the general matrix? Because it's a constraint?

That can be very expensive so it wouldn't be so useful in practice. It's also defeating the whole purpose of MOI of being able to transmit custom structure about the problem through JuMP. Even if we have an MOI meta solver that does that, it would need to define a new set to communicate the low-rank structure anyway.

[odow]

Do you have a numerical example of what the A matrix would be?

Also, (y, C) in isn't going to work because we don't support tuples as arguments.

[blegat]

Also, (y, C) in isn't going to work because we don't support tuples as arguments.

No, the JuMP syntax won't work, it's just to give an idea.

[odow]

Other questions:

• Clarify low-rank property of the input data vs low-rank constraints on the output variables
• How can the set be serialized to MOF?
• What solvers support this?

[blegat]

What solvers support this?

I try to keep the list in #2197

How can the set be serialized to MOF?

I think some special case might not be serializable but that's already the case for other parametrized sets like Scaled and Indicator. Same with MOI.Utililities.Model, we cannot aim to every every set there, other will be handled a bit less efficiently by UniversalFallback anyway, what's important there is to have the LP sets since LP solvers can handle millions of them

[odow]

@mlubin's concern around the serialization wasn't parameterized sets, but if we allow

struct ConeWithInnerProducts{A<:AbstractMatrix,S<:AbstractVectorSet}
    matrix::A
    set::S
end

How would we serialize matrix::A to JSON if we don't know the matrix type?

We need to think about doing something like:

struct ConeWithInnerProducts{T,S<:AbstractVectorSet}
    matrix::LinearAlgebra.LowerTriangular{Matrix{T}}
    set::S
end

[blegat]

If we have your former example, what prevents serializing ConeWithInnerProducts{LinearAlgebra.LowerTriangular{Matrix{T}}} ?

[odow]

Nothing. We could do that. But the point is that we'd need to pick the matrix type. We couldn't support arbitrary user-defined matrices.

[blegat]

But the point is that we'd need to pick the matrix type. We couldn't support arbitrary user-defined matrices.

Maybe we could only pick the matrix type for the purpose of MOF and MOI.Utilities.Model and leave it parametrized so as to allow additional flexibility for use cases that don't need writing to file and for which the performance of UniversalFallback is just fine. Of course, it will be easier to decide this once I have a few examples working :)

[blegat]

An additional point made by @mlubin that is worth writing down is that, for solvers supporting both low-rank solutions and low-rank constraints like SDPLR (and soon https://github.com/JuliaAlgebra/BMSOS.jl), you can start targeting SDP constraints ⟨A, X⟩ of size 1M of the rank of A and X is O(1). In that case, if A is low-rank but not sparse, communicating it in full and letting the solver do an SVD to recover the low-rank structure is not an option since the full A would have 10^12 entries.