Defective univariate example¶

In [1]:

using LinearAlgebra
using TypedPolynomials
using MacaulayMatrix
using MultivariateMoments
import Clarabel

Consider the following system with a root of multiplicity 4:

In [2]:

@polyvar x
system = [(x - 4)^4]

Out[2]:

1-element Vector{MultivariatePolynomials.VectorPolynomial{Int64, MultivariatePolynomials.Term{Int64, TypedPolynomials.Monomial{(x,), 1}}}}:
 256 - 256x + 96x² - 16x³ + x⁴

The Macaulay framework finds the root with degree 4:

In [3]:

solve_system(system, column_maxdegree = 4)

[ Info: Added 1 rows to complete columns up to degree 4
[ Info: Nullspace of dimensions (5, 4) computed from Macaulay matrix of dimension (1, 5) in 0.044666183 seconds.
[ Info: Found 1 real solution

Out[3]:

1-element Vector{Vector{Float64}}:
 [4.00000000000001]

Let's find the max rank PSD hankel from from the degree 4 Macaulay nullspace:

In [4]:

ν = moment_matrix(system, Clarabel.Optimizer, 4)
nothing #hide

[ Info: Nullspace of dimensions (5, 4) computed from Macaulay matrix of dimension (1, 5) in 4.1157e-5 seconds.
-------------------------------------------------------------
           Clarabel.jl v0.9.0  -  Clever Acronym              
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 4
  constraints   = 7
  nnz(P)        = 0
  nnz(A)        = 28
  cones (total) = 2
    : Zero        = 1,  numel = 1
    : PSDTriangle = 1,  numel = 6

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.73e-01  5.03e-01  1.00e+00  1.80e+00   ------   
  1   0.0000e+00   3.7385e-01  3.74e-01  1.03e-01  8.16e-02  6.35e-01  3.24e-01  8.37e-01  
  2   0.0000e+00   1.7155e+01  1.72e+01  2.30e-02  1.02e-02  1.83e+01  5.02e-02  9.71e-01  
  3   0.0000e+00  -2.6489e-01  2.65e-01  6.79e-03  5.26e-03  8.93e-02  3.12e-02  5.18e-01  
  4   0.0000e+00   6.0742e-03  6.07e-03  5.04e-04  3.93e-04  3.57e-02  2.45e-03  9.23e-01  
  5   0.0000e+00   2.2712e-01  2.27e-01  1.53e-04  8.38e-05  2.36e-01  4.68e-04  9.25e-01  
  6   0.0000e+00   2.5411e-01  2.54e-01  2.71e-05  2.49e-05  2.59e-01  1.76e-04  8.26e-01  
  7   0.0000e+00   6.3020e-02  6.30e-02  2.77e-06  2.37e-06  6.36e-02  1.61e-05  9.22e-01  
  8   0.0000e+00   1.4433e-02  1.44e-02  5.99e-07  3.95e-07  1.45e-02  2.38e-06  9.04e-01  
  9   0.0000e+00   1.4781e-03  1.48e-03  3.45e-08  2.27e-08  1.48e-03  1.37e-07  9.47e-01  
 10   0.0000e+00   1.3987e-05  1.40e-05  1.60e-09  1.11e-09  1.43e-05  6.89e-09  9.65e-01  
 11   0.0000e+00   7.8573e-07  7.86e-07  1.60e-10  1.21e-10  8.14e-07  6.34e-10  8.83e-01  
 12   0.0000e+00   1.7100e-07  1.71e-07  3.04e-11  9.66e-12  1.74e-07  6.93e-11  9.48e-01  
 13   0.0000e+00   1.5429e-07  1.54e-07  9.27e-07  1.18e-11  1.57e-07  2.92e-11  9.87e-02  
---------------------------------------------------------------------------------------------
Terminated with status = solved (reduced accuracy)
solve time =  397ms
[ Info: Terminated with ALMOST_OPTIMAL (ALMOST_SOLVED) in 0.397181296 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 : -1.71002e-07
│ 
│ * Work counters
│   Solve time (sec)   : 3.97181e-01
│   Barrier iterations : 13
└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:25

We find the following PSD hankel

In [5]:

ν

Out[5]:

MomentMatrix with row/column basis:
 MonomialBasis([1, x, x^2])
And entries in a 3×3 MultivariateMoments.SymMatrix{Float64}:
  1.0000000000000075   3.996902364627232   15.9763030351355
  3.996902364627232   15.9763030351355     63.86432262811344
 15.9763030351355     63.86432262811344   255.31107602137627

Looking at its singular values, the rank seems to be 1:

In [6]:

svd(value_matrix(ν)).S

Out[6]:

3-element Vector{Float64}:
 272.2861108836006
   0.001268176595828647
   3.684569158988992e-9

The rank of the truncation is also 1:

In [7]:

svd(value_matrix(truncate(ν, 1))).S

Out[7]:

2-element Vector{Float64}:
 16.976239739447575
  6.329568793473016e-5

We conclude that the PSD hankel is flat. The root can be obtained using:

In [8]:

atomic_measure(ν, FixedRank(1))

Out[8]:

Atomic measure on the variables x with 1 atoms:
 at [3.99742282946471] with weight 0.9998785045530533

The solution is also apparent from the equation in the nullspace of the moment matrix:

In [9]:

nullspace(ν, ShiftNullspace(), FixedRank(1)).polynomials

Out[9]:

1-element Vector{MultivariatePolynomials.VectorPolynomial{Float64, MultivariatePolynomials.Term{Float64, TypedPolynomials.Monomial{(x,), 1}}}}:
 -3.99742282946471 + x

This notebook was generated using Literate.jl.