Add support for SumOfSquares cone (#844)

* Add support for SumOfSquares cone

* Fix format

* Update to MOI interface

* Use LowRankOpt

* Fixes

* Fix format

* Drop Julia v1.6

* Add bound for LowRankOpt

* Add tests

* Bump doc

* Fixes

* Fix format

Author: blegat

Project.toml

@@ -10,6 +10,7 @@ GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
+LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 PolynomialRoots = "3a141323-8675-5d76-9d11-e1df1406c778"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -25,6 +26,7 @@ DocStringExtensions = "0.9"
 GenericLinearAlgebra = "0.3"
 IterativeSolvers = "0.9"
 LinearMaps = "3"
+LowRankOpt = "0.2"
 MathOptInterface = "1"
 PolynomialRoots = "1"
 Requires = "1"

src/Hypatia.jl

@@ -44,6 +44,7 @@ const VI = MOI.VariableIndex
 const SAF = MOI.ScalarAffineFunction
 const VV = MOI.VectorOfVariables
 const VAF = MOI.VectorAffineFunction
+import LowRankOpt as LRO
 include("MathOptInterface/cones.jl")
 include("MathOptInterface/transform.jl")
 include("MathOptInterface/wrapper.jl")

src/MathOptInterface/cones.jl

@@ -635,6 +635,32 @@ function cone_from_moi(
     return Cones.WSOSInterpNonnegative{T, R}(cone.U, cone.Ps, use_dual = cone.use_dual)
 end
 
+const _PrimalRankOnePSD{T <: Real, F <: AbstractVector{T}} = LRO.SetDotProducts{
+    LRO.WITHOUT_SET,
+    MOI.PositiveSemidefiniteConeTriangle,
+    LRO.TriangleVectorization{T, LRO.Factorization{T, F, LRO.One{T}}},
+}
+
+const _DualRankOnePSD{T <: Real, F <: AbstractVector{T}} = LRO.LinearCombinationInSet{
+    LRO.WITHOUT_SET,
+    MOI.PositiveSemidefiniteConeTriangle,
+    LRO.TriangleVectorization{T, LRO.Factorization{T, F, LRO.One{T}}},
+}
+
+function cone_from_moi(
+    ::Type{T},
+    cone::Union{_PrimalRankOnePSD, _DualRankOnePSD},
+) where {T <: Real}
+    return cone_from_moi(
+        T,
+        WSOSInterpNonnegativeCone{T, T}(
+            length(cone.vectors),
+            [reduce(vcat, [v.matrix.factor' for v in cone.vectors])],
+            cone isa _DualRankOnePSD,
+        ),
+    )
+end
+
 """
 $(TYPEDEF)
 
@@ -785,6 +811,8 @@ const SupportedCone{T <: Real} = Union{
     MOI.DualExponentialCone,
     MOI.LogDetConeTriangle,
     MOI.RelativeEntropyCone,
+    _PrimalRankOnePSD{T},
+    _DualRankOnePSD{T},
 }
 
 Base.copy(cone::HypatiaCones) = cone # maybe should deep copy the cone struct, but this is expensive

test/moicones.jl

@@ -13,8 +13,8 @@ using Test
 using LinearAlgebra
 import SparseArrays
 import Random
-import MathOptInterface
-const MOI = MathOptInterface
+import MathOptInterface as MOI
+import LowRankOpt as LRO
 import Hypatia
 import Hypatia.Cones
 
@@ -383,6 +383,37 @@ function test_moi_cones(T::Type{<:Real})
         @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
         @test hyp_cone.Ps == Ps
 
+        @testset "LowRankOpt" begin
+            P = Ps[1]
+            moi_cone = LRO.SetDotProducts{LRO.WITHOUT_SET}(
+                MOI.PositiveSemidefiniteConeTriangle(2),
+                LRO.TriangleVectorization.([
+                    LRO.positive_semidefinite_factorization(P[1, :]),
+                    LRO.positive_semidefinite_factorization(P[2, :]),
+                    LRO.positive_semidefinite_factorization(P[3, :]),
+                ]),
+            )
+            hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
+            @test hyp_cone isa Cones.WSOSInterpNonnegative{T, T}
+            @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
+            @test only(hyp_cone.Ps) == P
+            @test hyp_cone.use_dual_barrier
+
+            moi_cone = LRO.LinearCombinationInSet{LRO.WITHOUT_SET}(
+                MOI.PositiveSemidefiniteConeTriangle(2),
+                LRO.TriangleVectorization.([
+                    LRO.positive_semidefinite_factorization(P[1, :]),
+                    LRO.positive_semidefinite_factorization(P[2, :]),
+                    LRO.positive_semidefinite_factorization(P[3, :]),
+                ]),
+            )
+            hyp_cone = Hypatia.cone_from_moi(T, moi_cone)
+            @test hyp_cone isa Cones.WSOSInterpNonnegative{T, T}
+            @test MOI.dimension(moi_cone) == Cones.dimension(hyp_cone) == 3
+            @test only(hyp_cone.Ps) == P
+            @test !hyp_cone.use_dual_barrier
+        end
+
         Ps = [rand(Complex{T}, 4, 3), rand(Complex{T}, 4, 2)]
         moi_cone = Hypatia.WSOSInterpNonnegativeCone{T, Complex{T}}(4, Ps)
         hyp_cone = Hypatia.cone_from_moi(T, moi_cone)

test/runmoitests.jl

@@ -10,8 +10,8 @@ run MOI tests
 =#
 
 using Test
-import MathOptInterface
-const MOI = MathOptInterface
+import MathOptInterface as MOI
+import LowRankOpt as LRO
 import Hypatia
 include(joinpath(@__DIR__, "moicones.jl"))
 
@@ -27,6 +27,17 @@ include(joinpath(@__DIR__, "moicones.jl"))
         end
     end
 
+    @testset "Supports $T" for T in [Float64, BigFloat]
+        model = Hypatia.Optimizer{T}()
+        for S in [
+            Hypatia._PrimalRankOnePSD{T, Vector{T}}
+            Hypatia._DualRankOnePSD{T, Vector{T}}
+        ]
+            @test MOI.supports_add_constrained_variables(model, S)
+            @test MOI.supports_constraint(model, MOI.VectorAffineFunction{T}, S)
+        end
+    end
+
     # real types, tolerances, and tests to include for MOI.Test tests
     test_T = [
         (Float64, 2 * sqrt(sqrt(eps())), 4, String[], String[]),
@@ -49,19 +60,20 @@ include(joinpath(@__DIR__, "moicones.jl"))
             T,
         )
         MOI.set(model, MOI.Silent(), true)
+        config = MOI.Test.Config(
+            T,
+            atol = tol_test,
+            rtol = tol_test,
+            exclude = Any[
+                MOI.ConstraintBasisStatus,
+                MOI.VariableBasisStatus,
+                MOI.ObjectiveBound,
+                MOI.SolverVersion,
+            ],
+        )
         MOI.Test.runtests(
             model,
-            MOI.Test.Config(
-                T,
-                atol = tol_test,
-                rtol = tol_test,
-                exclude = Any[
-                    MOI.ConstraintBasisStatus,
-                    MOI.VariableBasisStatus,
-                    MOI.ObjectiveBound,
-                    MOI.SolverVersion,
-                ],
-            ),
+            config,
             include = includes,
             exclude = vcat(
                 excludes,
@@ -72,5 +84,6 @@ include(joinpath(@__DIR__, "moicones.jl"))
                 ],
             ),
         )
+        MOI.Test.runtests(model, config, test_module = LRO.Test)
     end
 end;