Defective univariate example
using LinearAlgebra
using TypedPolynomials
using MacaulayMatrix
using MultivariateMoments
import ClarabelConsider the following system with a root of multiplicity 4:
@polyvar x
system = [(x - 4)^4]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:
solve_system(system, column_maxdegree = 4)1-element Vector{Vector{Float64}}:
[4.00000000000001]Let's find the max rank PSD hankel from from the degree 4 Macaulay nullspace:
ν = moment_matrix(system, Clarabel.Optimizer, 4)[ Info: Nullspace of dimensions (5, 4) computed from Macaulay matrix of dimension (1, 5) in 3.792e-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 = 1.86ms
[ Info: Terminated with ALMOST_OPTIMAL (ALMOST_SOLVED) in 0.001863842 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) : 1.86384e-03
│ Barrier iterations : 13
└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:25We find the following PSD hankel
ν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.31107602137627Looking at its singular values, the rank seems to be 1:
svd(value_matrix(ν)).S3-element Vector{Float64}:
272.2861108836006
0.001268176595828647
3.684569158988992e-9The rank of the truncation is also 1:
svd(value_matrix(truncate(ν, 1))).S2-element Vector{Float64}:
16.976239739447575
6.329568793473016e-5We conclude that the PSD hankel is flat. The root can be obtained using:
atomic_measure(ν, FixedRank(1))Atomic measure on the variables x with 1 atoms:
at [3.99742282946471] with weight 0.9998785045530533The solution is also apparent from the equation in the nullspace of the moment matrix:
nullspace(ν, ShiftNullspace(), FixedRank(1)).polynomials1-element Vector{MultivariatePolynomials.VectorPolynomial{Float64, MultivariatePolynomials.Term{Float64, TypedPolynomials.Monomial{(x,), 1}}}}:
-3.99742282946471 + xThis page was generated using Literate.jl.