Univariate

using LinearAlgebra
using DynamicPolynomials
using MultivariateMoments
using JuMP
using MacaulayMatrix

Consider the following example:

@polyvar x
p = 3x^4 + 8x^3 - 6x^2 + 24x + 1
q = differentiate(p, x)

\[ 24 - 12x + 24x^{2} + 12x^{3} \]

Let's solve that system:

sols = solve_system([q], column_maxdegree = 4)

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

The real roots are:

sols

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

With PSD Hankel matrix trick no complex root!

import SCS
solver = SCS.Optimizer
psd_hankel([q], solver, 4)

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

What happened there ? First, we computed the MacaulayMatrix matrix nullspace

Z = nullspace(macaulay([q], 4))

MultivariateMoments.MacaulayNullspace{Float64, LinearAlgebra.Adjoint{Float64, Matrix{Float64}}, MultivariateBases.MonomialBasis{DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}([-0.5648877062932673 -0.5061871608167463 -0.08949732306804505; 0.35309169447067507 -0.5869837085383489 -0.3054118223094746; … ; -0.005487708235674751 0.6097339283112615 -0.15500888702836046; 0.048969435176441645 -0.13767209715388368 0.9351372741031582], MonomialBasis([1, x, x², x³, x⁴]), 5.438959822042073e-16)

The PSD hankel matrix we find is:

M = moment_matrix(Z, solver, 2)

MomentMatrix with row/column basis:
 MonomialBasis([1, x, x^2])
And entries in a 3×3 MultivariateMoments.SymMatrix{Float64}:
  1.0000000181608577   -2.6527137858823586    7.054717498400912
 -2.6527137858823586    7.054717498400912   -18.762148819005894
  7.054717498400912   -18.762148819005894    49.88444270817741

From which we get:

atomic_measure(M, 1e-6)

This page was generated using Literate.jl.