Positive dimensional roots at infinity

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)

2-element Vector{Vector{Float64}}:
 [0.4999999999999998, 0.49999999999999994, 0.8164965809277263, -0.816496580927726]
 [0.5, 0.4999999999999999, -0.816496580927726, 0.8164965809277255]

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)

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)

2-element Vector{Vector{Float64}}:
 [0.5000000000000001, 0.5000000000000003, 0.8164965809277283, -0.8164965809277278]
 [0.5000000000000002, 0.5000000000000001, -0.8164965809277281, 0.8164965809277276]

This page was generated using Literate.jl.