{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Defective univariate example"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using LinearAlgebra\n",
"using TypedPolynomials\n",
"using MacaulayMatrix\n",
"using MultivariateMoments\n",
"import Clarabel"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"Consider the following system with a root of multiplicity 4:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "1-element Vector{MultivariatePolynomials.VectorPolynomial{Int64, MultivariatePolynomials.Term{Int64, TypedPolynomials.Monomial{(x,), 1}}}}:\n 256 - 256x + 96x² - 16x³ + x⁴"
},
"metadata": {},
"execution_count": 2
}
],
"cell_type": "code",
"source": [
"@polyvar x\n",
"system = [(x - 4)^4]"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"The Macaulay framework finds the root with degree 4:"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ Info: Added 1 rows to complete columns up to degree 4\n",
"[ Info: Nullspace of dimensions (5, 4) computed from Macaulay matrix of dimension (1, 5) in 0.044666183 seconds.\n",
"[ Info: Found 1 real solution\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "1-element Vector{Vector{Float64}}:\n [4.00000000000001]"
},
"metadata": {},
"execution_count": 3
}
],
"cell_type": "code",
"source": [
"solve_system(system, column_maxdegree = 4)"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"Let's find the max rank PSD hankel from from the degree 4\n",
"Macaulay nullspace:"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ Info: Nullspace of dimensions (5, 4) computed from Macaulay matrix of dimension (1, 5) in 4.1157e-5 seconds.\n",
"-------------------------------------------------------------\n",
" Clarabel.jl v0.9.0 - Clever Acronym \n",
" (c) Paul Goulart \n",
" University of Oxford, 2022 \n",
"-------------------------------------------------------------\n",
"\n",
"problem:\n",
" variables = 4\n",
" constraints = 7\n",
" nnz(P) = 0\n",
" nnz(A) = 28\n",
" cones (total) = 2\n",
" : Zero = 1, numel = 1\n",
" : PSDTriangle = 1, numel = 6\n",
"\n",
"settings:\n",
" linear algebra: direct / qdldl, precision: Float64\n",
" max iter = 200, time limit = Inf, max step = 0.990\n",
" tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,\n",
" static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32\n",
" dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07\n",
" iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12, \n",
" max iter = 10, stop ratio = 5.0\n",
" equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04\n",
" max iter = 10\n",
"\n",
"iter pcost dcost gap pres dres k/t μ step \n",
"---------------------------------------------------------------------------------------------\n",
" 0 0.0000e+00 -0.0000e+00 0.00e+00 4.73e-01 5.03e-01 1.00e+00 1.80e+00 ------ \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
" 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 \n",
"---------------------------------------------------------------------------------------------\n",
"Terminated with status = solved (reduced accuracy)\n",
"solve time = 397ms\n",
"[ Info: Terminated with ALMOST_OPTIMAL (ALMOST_SOLVED) in 0.397181296 seconds.\n",
"┌ Warning: * Solver : Clarabel\n",
"│ \n",
"│ * Status\n",
"│ Result count : 1\n",
"│ Termination status : ALMOST_OPTIMAL\n",
"│ Message from the solver:\n",
"│ \"ALMOST_SOLVED\"\n",
"│ \n",
"│ * Candidate solution (result #1)\n",
"│ Primal status : NEARLY_FEASIBLE_POINT\n",
"│ Dual status : NEARLY_FEASIBLE_POINT\n",
"│ Objective value : -0.00000e+00\n",
"│ Dual objective value : -1.71002e-07\n",
"│ \n",
"│ * Work counters\n",
"│ Solve time (sec) : 3.97181e-01\n",
"│ Barrier iterations : 13\n",
"└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:25\n"
]
}
],
"cell_type": "code",
"source": [
"ν = moment_matrix(system, Clarabel.Optimizer, 4)\n",
"nothing #hide"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"We find the following PSD hankel"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "MomentMatrix with row/column basis:\n MonomialBasis([1, x, x^2])\nAnd entries in a 3×3 MultivariateMoments.SymMatrix{Float64}:\n 1.0000000000000075 3.996902364627232 15.9763030351355\n 3.996902364627232 15.9763030351355 63.86432262811344\n 15.9763030351355 63.86432262811344 255.31107602137627"
},
"metadata": {},
"execution_count": 5
}
],
"cell_type": "code",
"source": [
"ν"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"Looking at its singular values, the rank seems to be 1:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "3-element Vector{Float64}:\n 272.2861108836006\n 0.001268176595828647\n 3.684569158988992e-9"
},
"metadata": {},
"execution_count": 6
}
],
"cell_type": "code",
"source": [
"svd(value_matrix(ν)).S"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"The rank of the truncation is also 1:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "2-element Vector{Float64}:\n 16.976239739447575\n 6.329568793473016e-5"
},
"metadata": {},
"execution_count": 7
}
],
"cell_type": "code",
"source": [
"svd(value_matrix(truncate(ν, 1))).S"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"We conclude that the PSD hankel is flat.\n",
"The root can be obtained using:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "Atomic measure on the variables x with 1 atoms:\n at [3.99742282946471] with weight 0.9998785045530533"
},
"metadata": {},
"execution_count": 8
}
],
"cell_type": "code",
"source": [
"atomic_measure(ν, FixedRank(1))"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"The solution is also apparent from the equation in the nullspace of the moment matrix:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "1-element Vector{MultivariatePolynomials.VectorPolynomial{Float64, MultivariatePolynomials.Term{Float64, TypedPolynomials.Monomial{(x,), 1}}}}:\n -3.99742282946471 + x"
},
"metadata": {},
"execution_count": 9
}
],
"cell_type": "code",
"source": [
"nullspace(ν, ShiftNullspace(), FixedRank(1)).polynomials"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"---\n",
"\n",
"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.0"
},
"kernelspec": {
"name": "julia-1.11",
"display_name": "Julia 1.11.0",
"language": "julia"
}
},
"nbformat": 4
}