Adapted from: Example 6.14 p. 102 of [D13]
[D13] Dreesen, Philippe. Back to the Roots: Polynomial System Solving Using Linear Algebra Ph.D. thesis (2013)
using LinearAlgebra
using TypedPolynomials
using MacaulayMatrix
using JuMP
using MultivariateMoments
Consider the system given in [D13, Example 6.14] which corresponds to Systems/dreesen10 of macaulaylab:
@polyvar x[1:4]
system = [
x[1] + x[2] - 1,
x[1] * x[3] + x[2] * x[4],
x[1] * x[3]^2 + x[2] * x[4]^2 - 2/3,
x[1] * x[3]^3 + x[2] * x[4]^3,
]
4-element Vector{MultivariatePolynomials.VectorPolynomial{Float64, MultivariatePolynomials.Term{Float64, TypedPolynomials.Monomial{(x₁, x₂, x₃, x₄), 4}}}}:
-1.0 + x₂ + x₁
x₂x₄ + x₁x₃
-0.6666666666666666 + x₂x₄² + x₁x₃²
x₂x₄³ + x₁x₃³
With the classical MacaulayMatrix approach, the nullity increases by 2 at every degree because of the positive dimensional solution set at infinity.
solve_system(system, column_maxdegree = 6)
[ Info: Added 56 rows to complete columns up to degree 4 [ Info: Nullspace of dimensions (65, 15) computed from Macaulay matrix of dimension (56, 65) in 0.000884172 seconds. [ Info: Added 69 rows to complete columns up to degree 5 [ Info: Nullspace of dimensions (120, 17) computed from Macaulay matrix of dimension (125, 120) in 0.001944609 seconds. [ Info: Added 121 rows to complete columns up to degree 6 [ Info: Nullspace of dimensions (203, 18) computed from Macaulay matrix of dimension (246, 203) in 0.007104496 seconds. [ Info: Found 2 real solutions
2-element Vector{Vector{Float64}}:
[0.4999999999999999, 0.5000000000000001, 0.8164965809277263, -0.8164965809277261]
[0.5000000000000001, 0.5, -0.816496580927726, 0.8164965809277256]
Let's compute the moment matrix of degree 4 from the nullspace of the Macaulay matrix of degree 5 using the Clarabel SDP solver:
import Clarabel
ν = moment_matrix(system, Clarabel.Optimizer, 5)
[ Info: Nullspace of dimensions (120, 17) computed from Macaulay matrix of dimension (125, 120) in 0.001968263 seconds.
-------------------------------------------------------------
Clarabel.jl v0.9.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 17
constraints = 121
nnz(P) = 0
nnz(A) = 2025
cones (total) = 2
: Zero = 1, numel = 1
: PSDTriangle = 1, numel = 120
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 0.0000e+00 -0.0000e+00 0.00e+00 4.01e-01 3.42e+00 1.00e+00 3.22e+00 ------
1 0.0000e+00 -1.1323e-03 1.13e-03 5.81e-02 4.66e-01 1.19e-01 4.64e-01 8.61e-01
2 0.0000e+00 2.8687e-04 2.87e-04 1.56e-03 1.23e-02 3.48e-03 1.27e-02 9.74e-01
3 0.0000e+00 2.8697e-06 2.87e-06 1.56e-05 1.23e-04 3.48e-05 1.28e-04 9.90e-01
4 0.0000e+00 2.8697e-08 2.87e-08 1.56e-07 1.23e-06 3.48e-07 1.28e-06 9.90e-01
5 0.0000e+00 2.8697e-10 2.87e-10 1.56e-09 1.23e-08 3.48e-09 1.28e-08 9.90e-01
6 0.0000e+00 2.8697e-10 2.87e-10 1.56e-09 1.23e-08 3.48e-09 1.28e-08 0.00e+00
---------------------------------------------------------------------------------------------
Terminated with status = solved (reduced accuracy)
solve time = 12.2ms
[ Info: Terminated with ALMOST_OPTIMAL (ALMOST_SOLVED) in 0.012238786000000001 seconds.
┌ Warning: * Solver : Clarabel
│
│ * Status
│ Result count : 1
│ Termination status : ALMOST_OPTIMAL
│ Message from the solver:
│ "ALMOST_SOLVED"
│
│ * Candidate solution (result #1)
│ Primal status : NEARLY_FEASIBLE_POINT
│ Dual status : NEARLY_FEASIBLE_POINT
│ Objective value : -0.00000e+00
│ Dual objective value : -2.86973e-10
│
│ * Work counters
│ Solve time (sec) : 1.22388e-02
│ Barrier iterations : 6
└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:25
MomentMatrix with row/column basis:
MonomialBasis([1, x[4], x[3], x[2], x[1], x[4]^2, x[3]*x[4], x[3]^2, x[2]*x[4], x[2]*x[3], x[2]^2, x[1]*x[4], x[1]*x[3], x[1]*x[2], x[1]^2])
And entries in a 15×15 MultivariateMoments.SymMatrix{Float64}:
1.0 3.725744452537372e-10 … 0.2500000002434462
3.725744452537372e-10 0.6666666666666679 9.315270776966145e-11
-3.7257428739390086e-10 -0.6666666666666662 -9.313425725077096e-11
0.5000000000000245 1.8628763028688544e-10 0.12499999987823852
0.499999999999976 1.8628704828715925e-10 0.1250000003652074
0.6666666666666679 2.4838786583103456e-10 … 0.16666666650438633
-0.6666666666666662 -2.4838329223259015e-10 -0.16666666682896436
0.6666666666666669 2.4838344055144734e-10 0.1666666665043526
1.8628763028688544e-10 0.333333333333317 4.6547795999607455e-11
-1.8628716624835562e-10 -0.3333333333333493 -4.6587067537018356e-11
0.25000000024349434 9.315274246413097e-11 … 0.2555649744273872
1.8628704828715925e-10 0.33333333333335063 4.660461599970134e-11
-1.8628748023330477e-10 -0.333333333333317 -4.654744515178444e-11
0.2499999997565302 9.313424337298315e-11 -0.13056497454914884
0.2500000002434462 9.315270776966145e-11 0.2555649749143563
The nullspace of this moment matrix is a Macaulay matrix
M = nullspace(ν, ShiftNullspace())
0×0 Macaulay matrix for polynomials: 5.897711921664038e-14 + 0.9999999999999994*x[4] + x[3] - 1.4682525325425826e-14*x[4]^2 -0.5000000000001729 + 3.2596772937398524e-16*x[4] + x[2] + 1.0484472896961149e-13*x[4]^2 -0.49999999999983163 + 6.156031452731937e-17*x[4] + x[1] - 1.0394986062844766e-13*x[4]^2 0.666666666666666 + 5.268577352792162e-16*x[4] - 2.2195756829484477e-16*x[4]^2 + x[3]*x[4] -1.536525535226813e-14 - 0.49999999999999983*x[4] + 9.01468757230224e-15*x[4]^2 + x[2]*x[4] 1.6090974574336504e-14 - 0.5000000000000008*x[4] - 1.0732798282105085e-14*x[4]^2 + x[1]*x[4] The row shifts are: MonomialBasis([]) The column basis is: MonomialBasis([])
Each coefficient is close to a rational number, after some clean up, we get the following system:
real_system = clean(M; tol = 1e-6).polynomials
6-element Vector{MultivariatePolynomials.VectorPolynomial{Float64, MultivariatePolynomials.Term{Float64, TypedPolynomials.Monomial{(x₁, x₂, x₃, x₄), 4}}}}:
x₄ + x₃
-1.0 + 2.0x₂
-1.0 + 2.0x₁
2.0 + 3.0x₃x₄
-x₄ + 2.0x₂x₄
-x₄ + 2.0x₁x₄
We have the equations 2x[1] = 1 and 2x[2] = 1 from which we get
x[1] = 1/2 and x[2] = 1/2.
We then have the equations x[3] = -x[4] and 3x[3] * x[4] = -2 from
which we get x[3] = ±√(2/3) and x[4] = ∓√(2/3).
Notice that the complex solutions are not solutions of the system anymore.
The Macaulay solver indeed find these real solutions direction at degree 2.
solve_system(real_system, column_maxdegree = 2)
[ Info: Added 18 rows to complete columns up to degree 2 [ Info: Nullspace of dimensions (15, 2) computed from Macaulay matrix of dimension (18, 15) in 9.5116e-5 seconds. [ Info: Found 2 real solutions
2-element Vector{Vector{Float64}}:
[0.4999999999999999, 0.5000000000000001, 0.8164965809277269, -0.8164965809277265]
[0.5, 0.4999999999999998, -0.8164965809277267, 0.8164965809277263]
This notebook was generated using Literate.jl.