Getting rid of root at infinity¶

Adapted from: Section~2.2.1 p. 33 of [D13]

[D13] Dreesen, Philippe. Back to the Roots: Polynomial System Solving Using Linear Algebra Ph.D. thesis (2013)

In [1]:

using LinearAlgebra
using TypedPolynomials
using MacaulayMatrix
using JuMP
using MultivariateMoments

Consider the system given in [D13, (2.3)] which corresponds to Systems/dreesen2 of macaulaylab:

In [2]:

@polyvar x y z
system = [
    x^2 - x*y + z,
    2y^3 - 2x*y^2 - 3x * y,
    z^3 - x*y*z - 2,
]

Out[2]:

3-element Vector{MultivariatePolynomials.VectorPolynomial{Int64, MultivariatePolynomials.Term{Int64, TypedPolynomials.Monomial{(x, y, z), 3}}}}:
 z - xy + x²
 -3xy + 2y³ - 2xy²
 -2 + z³ - xyz

We first try to solve the system:

In [3]:

sols = solve_system(system, column_maxdegree = 6)
nothing # hide

[ Info: Added 6 rows to complete columns up to degree 3
[ Info: Nullspace of dimensions (14, 8) computed from Macaulay matrix of dimension (6, 14) in 8.4166e-5 seconds.
[ Info: Added 12 rows to complete columns up to degree 4
[ Info: Nullspace of dimensions (33, 15) computed from Macaulay matrix of dimension (18, 33) in 0.000132997 seconds.
[ Info: Added 22 rows to complete columns up to degree 5
[ Info: Nullspace of dimensions (56, 18) computed from Macaulay matrix of dimension (40, 56) in 0.000318782 seconds.
[ Info: Added 35 rows to complete columns up to degree 6
[ Info: Nullspace of dimensions (84, 18) computed from Macaulay matrix of dimension (75, 84) in 0.000789476 seconds.
[ Info: Found 1 real solution

The real solutions are

In [4]:

sols

Out[4]:

1-element Vector{Vector{Float64}}:
 [3.8420560175113363, 4.99567150625108, 4.432255330406793]

Staicase analysis¶

In [5]:

solver = Iterator(system, MacaulayMatrix.Solver())
step!(solver)

[ Info: Added 6 rows to complete columns up to degree 3
[ Info: Nullspace of dimensions (14, 8) computed from Macaulay matrix of dimension (6, 14) in 4.1076e-5 seconds.

We can look at the solver as follows:

In [6]:

solver

Out[6]:

MacaulayMatrix matrix solver. Last iteration considered:
6×14 Macaulay matrix for polynomials:
  z - x*y + x^2
  -3*x*y + 2*y^3 - 2*x*y^2
  -2 + z^3 - x*y*z
The row shifts are:
MonomialBasis([1, z, y, x])
The column basis is:
MonomialBasis([1, z, z^2, y*z, x*z, x*y, x^2, z^3, y^3, x*y*z, x*y^2, x^2*z, x^2*y, x^3])
BorderBasis with independent rows and dependent columns in:
BasisDependence for bases:
 Standard:
 MonomialBasis([1, z, z^2, y*z, x*z, x*y, z^3, y^3])
 Trivial Standard:
 MonomialBasis([y, x, y^2, y*z^2, y^2*z, x*z^2, z^4, y*z^3, y^2*z^2, y^3*z, y^4, x*z^3])
 Corners:
 MonomialBasis([x^2, x*y*z, x*y^2])
 Trivial Independent Border:
 MonomialBasis([x*y*z^2, x*y^2*z, x*y^3, x^2*z^2])
 Dependent Border:
 MonomialBasis([x^2*z, x^2*y])

And entries in a 8×5 adjoint(::Matrix{Float64}) with eltype Float64:
  4.159065386202991e-16   -1.9999999999999998      …  -1.9069420129706838e-16
 -1.0                     -3.1080979595180984e-16      4.5135339788020976e-17
  1.1499638451679124e-17   2.176915759031188e-16      -1.2298933558687206e-16
 -1.041204178188655e-16   -3.470680593962184e-16      -1.0
 -7.832624135939215e-17   -2.4017185283265312e-17     -4.699574481563529e-16
  1.0                     -4.3281882871179737e-17  …  -1.5000000000000002
  4.8138576458356415e-17   1.0                        -3.677747436866566e-17
 -5.551115123125784e-17    3.8163916471489786e-16      1.0000000000000004
Current status is OPTIMIZE_NOT_CALLED
History of iterations:
1×3 DataFrame
 Row │ nullity  num_rows  num_cols
     │ Int64    Int64     Int64
─────┼─────────────────────────────
   1 │       8         6        14

We can see in the border dependence that even if x^2 is a corner, the multiplication matrix for x cannot be computed as y^3 (resp. x*z^2, y*z^2) is standard but x * y^3 (resp. x^2 * z^2, x * y * z^2) is indepedent.

In [7]:

using Plots
plot(saturated_dependence(solver))

Out[7]:

Let's do another step:

In [8]:

step!(solver)

[ Info: Added 12 rows to complete columns up to degree 4
[ Info: Nullspace of dimensions (33, 15) computed from Macaulay matrix of dimension (18, 33) in 0.000138978 seconds.

This time, the blocker for computing the multiplication matrix for x are z^4 and y * z^3 among others.

In [9]:

plot(saturated_dependence(solver))

Out[9]:

Let's do another step:

In [10]:

step!(solver)

[ Info: Added 22 rows to complete columns up to degree 5
[ Info: Nullspace of dimensions (56, 18) computed from Macaulay matrix of dimension (40, 56) in 0.000340843 seconds.

Now they are z^5 and y * z^4 among others.

In [11]:

plot(saturated_dependence(solver))

Out[11]:

Let's do a last step:

In [12]:

step!(solver)

[ Info: Added 35 rows to complete columns up to degree 6
[ Info: Nullspace of dimensions (84, 18) computed from Macaulay matrix of dimension (75, 84) in 0.000839189 seconds.
[ Info: Found 1 real solution

Now we see that the whole border is dependent so the four multiplication matrices can be computed.

In [13]:

plot(saturated_dependence(solver))

Out[13]:

In retrospect, we we probably should have expanded in priority towards larger exponents for z.

Saturation¶

Let's start again by starting with the same first stop.

In [14]:

solver = Iterator(system, MacaulayMatrix.Solver())
step!(solver)

[ Info: Added 6 rows to complete columns up to degree 3
[ Info: Nullspace of dimensions (14, 8) computed from Macaulay matrix of dimension (6, 14) in 5.8619e-5 seconds.

This time, let's focus on saturating the first non-saturated standard monomials:

In [15]:

step!(solver, FirstStandardNonSaturated(10))

[ Info: Added 12 rows to saturate columns `TypedPolynomials.Monomial{(x, y, z), 3}[z, y, x, z², yz, y², xz, xy, z³, yz²]`
[ Info: Nullspace of dimensions (38, 20) computed from Macaulay matrix of dimension (18, 38) in 0.000134219 seconds.

We can see that z is now saturated:

In [16]:

plot(saturated_dependence(solver))

Out[16]:

Let's saturate the next one:

In [17]:

step!(solver, FirstStandardNonSaturated(10))

[ Info: Added 20 rows to saturate columns `TypedPolynomials.Monomial{(x, y, z), 3}[y²z, y³, xz², z⁴, yz³, y²z², y³z, y⁴, z⁵, yz⁴]`
[ Info: Nullspace of dimensions (70, 32) computed from Macaulay matrix of dimension (38, 70) in 0.00034996 seconds.

We can see that y went from trivial to saturated:

In [18]:

plot(saturated_dependence(solver))

Out[18]:

Let's saturate the next one:

In [19]:

step!(solver, FirstStandardNonSaturated(10))

[ Info: Added 20 rows to saturate columns `TypedPolynomials.Monomial{(x, y, z), 3}[y⁵, z⁶, yz⁵, y⁶, z⁷, yz⁶, y⁷, z⁸, yz⁷, y⁸]`
[ Info: Nullspace of dimensions (103, 45) computed from Macaulay matrix of dimension (58, 103) in 0.000688619 seconds.

We can see that x went from trivial to saturated:

In [20]:

plot(saturated_dependence(solver))

Out[20]:

Let's saturate the next one:

In [21]:

step!(solver, FirstStandardNonSaturated(10))

[ Info: Added 20 rows to saturate columns `TypedPolynomials.Monomial{(x, y, z), 3}[z⁹, yz⁸, y⁹, z¹⁰, yz⁹, y¹⁰, z¹¹, yz¹⁰, y¹¹, z¹²]`
[ Info: Nullspace of dimensions (136, 58) computed from Macaulay matrix of dimension (78, 136) in 0.001223162 seconds.

We can see that z^2 is now saturated:

In [22]:

plot(saturated_dependence(solver))

Out[22]:

We can see that this does not seem to work. By saturating a monomial, we generate all the rows that have a nonzero entry for the corresponding column in the MacaulayMatrix matrix. If the monomial is standard, it means that the column is a linear combination of other columns (which are dependent in the MacaulayMatrix Nullspace) of the MacaulayMatrix matrix. If we saturate these other columns, the linear combination may not work anymore. But we are only saturating the standard monomials, we don't saturate the rows that are dependent in the MacaulayMatrix Nullspace which explains the result we have seen here.

Moment approach¶

With moment matrix of degree 3:

In [23]:

import Clarabel
solver = Clarabel.Optimizer

M = moment_matrix(system, solver, 3)
nothing # hide

atomic_measure(M, 1e-4, ShiftNullspace())

[ Info: Nullspace of dimensions (14, 8) computed from Macaulay matrix of dimension (6, 14) in 7.4087e-5 seconds.
┌ Warning: Missing monomials TypedPolynomials.Monomial{(x, y, z), 3}[y, x, y²] in Macaulay nullspace, leaving them free to take any value in the moment matrix then.
└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:65
-------------------------------------------------------------
           Clarabel.jl v0.9.0  -  Clever Acronym              
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 11
  constraints   = 11
  nnz(P)        = 0
  nnz(A)        = 67
  cones (total) = 2
    : Zero        = 1,  numel = 1
    : PSDTriangle = 1,  numel = 10

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  3.88e-01  6.59e-01  1.00e+00  1.24e+00   ------   
  1   0.0000e+00  -6.5201e-02  6.52e-02  3.19e-02  7.78e-02  7.16e-03  1.25e-01  9.13e-01  
  2   0.0000e+00  -5.8123e-04  5.81e-04  3.43e-04  9.79e-04  1.89e-04  1.46e-03  9.89e-01  
  3   0.0000e+00  -5.8069e-06  5.81e-06  3.43e-06  9.79e-06  1.89e-06  1.46e-05  9.90e-01  
  4   0.0000e+00  -5.8068e-08  5.81e-08  3.43e-08  9.79e-08  1.89e-08  1.46e-07  9.90e-01  
  5   0.0000e+00  -5.8068e-10  5.81e-10  3.43e-10  9.79e-10  1.89e-10  1.46e-09  9.90e-01  
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time =  933μs
[ Info: Terminated with OPTIMAL (SOLVED) in 0.000933484 seconds.

With moment matrix of degree 4:

In [24]:

M = moment_matrix(system, solver, 4)
nothing # hide

atomic_measure(M, 1e-4, ShiftNullspace())

[ Info: Nullspace of dimensions (33, 15) computed from Macaulay matrix of dimension (18, 33) in 0.000128759 seconds.
┌ Warning: Missing monomials TypedPolynomials.Monomial{(x, y, z), 3}[y², y²z²] in Macaulay nullspace, leaving them free to take any value in the moment matrix then.
└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:65
-------------------------------------------------------------
           Clarabel.jl v0.9.0  -  Clever Acronym              
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 17
  constraints   = 56
  nnz(P)        = 0
  nnz(A)        = 784
  cones (total) = 2
    : Zero        = 1,  numel = 1
    : PSDTriangle = 1,  numel = 55

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  7.02e-01  1.07e+00  1.00e+00  4.05e+00   ------   
  1   0.0000e+00   2.1244e+00  2.12e+00  2.49e-01  2.23e-01  2.83e+00  8.73e-01  8.09e-01  
  2   0.0000e+00   1.0127e+01  1.01e+01  4.85e-02  4.27e-02  1.08e+01  2.18e-01  8.63e-01  
  3   0.0000e+00   2.9803e+01  2.98e+01  7.30e-03  7.27e-03  3.02e+01  4.34e-02  8.59e-01  
  4   0.0000e+00   1.4744e+01  1.47e+01  1.30e-03  1.40e-03  1.49e+01  9.54e-03  8.22e-01  
  5   0.0000e+00   3.0716e+00  3.07e+00  2.61e-04  3.30e-04  3.11e+00  3.23e-03  8.96e-01  
  6   0.0000e+00   8.9652e-02  8.97e-02  1.35e-05  1.74e-05  9.14e-02  1.89e-04  9.44e-01  
  7   0.0000e+00   9.2794e-04  9.28e-04  1.36e-07  1.76e-07  9.46e-04  1.95e-06  9.90e-01  
  8   0.0000e+00   9.2781e-06  9.28e-06  1.36e-09  1.76e-09  9.46e-06  1.95e-08  9.90e-01  
  9   0.0000e+00   9.2783e-08  9.28e-08  1.35e-11  1.76e-11  9.46e-08  1.95e-10  9.90e-01  
 10   0.0000e+00   1.0031e-09  1.00e-09  2.78e-09  1.90e-13  1.02e-09  2.09e-12  9.89e-01  
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 3.27ms
[ Info: Terminated with OPTIMAL (SOLVED) in 0.0032739120000000003 seconds.

With moment matrix of degree 5:

In [25]:

M = moment_matrix(system, solver, 5)
nothing # hide

atomic_measure(M, 1e-4, ShiftNullspace())

[ Info: Nullspace of dimensions (56, 18) computed from Macaulay matrix of dimension (40, 56) in 0.000325575 seconds.
-------------------------------------------------------------
           Clarabel.jl v0.9.0  -  Clever Acronym              
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 18
  constraints   = 56
  nnz(P)        = 0
  nnz(A)        = 1008
  cones (total) = 2
    : Zero        = 1,  numel = 1
    : PSDTriangle = 1,  numel = 55

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  7.25e-01  1.01e+00  1.00e+00  4.89e+00   ------   
  1   0.0000e+00   3.7044e+00  3.70e+00  2.92e-01  1.85e-01  4.54e+00  8.83e-01  8.53e-01  
  2   0.0000e+00   2.4003e+02  2.40e+02  2.57e-02  1.15e-02  2.43e+02  5.85e-02  9.57e-01  
  3   0.0000e+00   1.4580e+03  1.46e+03  3.15e-03  1.16e-03  1.46e+03  5.54e-03  9.48e-01  
  4   0.0000e+00   7.8960e+00  7.90e+00  4.79e-04  2.89e-04  8.10e+00  1.88e-03  7.49e-01  
  5   0.0000e+00   2.2070e-01  2.21e-01  1.50e-05  9.30e-06  2.28e-01  6.24e-05  9.67e-01  
  6   0.0000e+00   1.5989e-02  1.60e-02  4.61e-07  2.90e-07  1.62e-02  1.98e-06  9.69e-01  
  7   0.0000e+00   3.2459e-04  3.25e-04  1.08e-08  6.90e-09  3.30e-04  4.76e-08  9.76e-01  
  8   0.0000e+00   3.2459e-04  3.25e-04  1.08e-08  6.90e-09  3.30e-04  4.76e-08  0.00e+00  
---------------------------------------------------------------------------------------------
Terminated with status = numerical error
solve time = 2.62ms
[ Info: Terminated with NUMERICAL_ERROR (NUMERICAL_ERROR) in 0.002618824 seconds.
┌ Warning: * Solver : Clarabel
│ 
│ * Status
│   Result count       : 1
│   Termination status : NUMERICAL_ERROR
│   Message from the solver:
│   "NUMERICAL_ERROR"
│ 
│ * Candidate solution (result #1)
│   Primal status      : OTHER_RESULT_STATUS
│   Dual status        : OTHER_RESULT_STATUS
│   Objective value    : -0.00000e+00
│   Dual objective value : -3.24587e-04
│ 
│ * Work counters
│   Solve time (sec)   : 2.61882e-03
│   Barrier iterations : 8
└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:25

Out[25]:

Atomic measure on the variables x, y, z with 1 atoms:
 at [3.84205516042132, 4.99567052245059, 4.432254383584636] with weight 0.9999998706390161

This notebook was generated using Literate.jl.