Fix Explicit approach (#45)

* Fix Explicit approach

* Fix format

Author: blegat

Project.toml

@@ -25,6 +25,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 
 [compat]

src/MacaulayMatrix.jl

@@ -8,6 +8,8 @@ import JuMP
 import MultivariateMoments as MM
 import Printf
 
+import SumOfSquares
+
 import CommonSolve: solve, solve!, init, step!
 
 abstract type AbstractIteration end

src/hankel.jl

@@ -121,9 +121,19 @@ function MM.moment_matrix(
 )
     model = JuMP.Model(solver)
     JuMP.@variable(model, γ)
-    JuMP.@objective(model, Min, γ)
-    JuMP.@variable(model, λ[i in eachindex(polynomials)])
-    con_ref = JuMP.@constraint(model, γ + dot(polynomials, λ) in SOSCone())
+    JuMP.@objective(model, Max, γ)
+    vars = MP.variables(polynomials)
+    JuMP.@variable(
+        model,
+        λ[i in eachindex(polynomials)],
+        SumOfSquares.PolyJuMP.Poly(
+            MP.monomials(vars, 0:(maxdegree-MP.maxdegree(polynomials[i]))),
+        ),
+    )
+    con_ref = JuMP.@constraint(
+        model,
+        LinearAlgebra.dot(polynomials, λ) - γ in SumOfSquares.SOSCone()
+    )
     JuMP.optimize!(model)
     if JuMP.termination_status(model) == JuMP.MOI.DUAL_INFEASIBLE
         return

test/Project.toml

@@ -1,5 +1,5 @@
 [deps]
-CSDP = "0a46da34-8e4b-519e-b418-48813639ff34"
+Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacaulayMatrix = "bcb046b7-f005-4e9a-a986-fb3acd014398"

test/macaulay.jl

@@ -6,7 +6,7 @@ import MultivariateBases as MB
 import MultivariateMoments as MM
 using MacaulayMatrix
 using JuMP
-import CSDP
+import Clarabel
 
 function __test_monomial_ideal_generators(standard, generators)
     basis = MB.MonomialBasis(standard)
@@ -48,14 +48,16 @@ function test_monomial_ideal_generators()
 end
 
 function _test_moments(polys, solver, d, expected; tol = 1e-4)
-    ν = MM.moment_matrix(polys, solver, d)
-    @test MM.value_matrix(ν)[1, 1] ≈ 1 rtol = tol
-    η = MM.atomic_measure(ν, tol)
-    if isnothing(expected)
-        @test isnothing(η)
-    else
-        sols = [η.atoms[i].center for i in eachindex(η.atoms)]
-        _test_sols(sols, expected)
+    @testset "$approach" for approach in [Explicit(), Hankel()]
+        ν = MM.moment_matrix(polys, solver, d, approach)
+        @test MM.value_matrix(ν)[1, 1] ≈ 1 rtol = tol
+        η = MM.atomic_measure(ν, tol)
+        if isnothing(expected)
+            @test isnothing(η)
+        else
+            sols = [η.atoms[i].center for i in eachindex(η.atoms)]
+            _test_sols(sols, expected)
+        end
     end
 end
 
@@ -114,7 +116,10 @@ function test_dreesen1()
         @testset "psd_hankel" begin
             _test_moments(
                 ps,
-                optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true),
+                optimizer_with_attributes(
+                    Clarabel.Optimizer,
+                    MOI.Silent() => true,
+                ),
                 d,
                 d == 3 ? nothing : expected,
             )
@@ -134,11 +139,9 @@ function test_univariate()
         sols = solve_system([q], column_maxdegree = d)
         _test_sols(sols, [[exp]])
     end
-    solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
-    @test psd_hankel([q], solver, 3) === nothing
-    @testset "d=$d" for d in 4:8
-        sols = psd_hankel([q], solver, d)
-        _test_sols(sols, [[exp]])
+    solver = optimizer_with_attributes(Clarabel.Optimizer, MOI.Silent() => true)
+    @testset "d=$d" for d in 3:8
+        _test_moments([q], solver, d, d == 3 ? nothing : [[exp]])
     end
 end