{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Positive dimensional roots at infinity\n",
"**Adapted from**: Example 6.14 p. 102 of [D13]\n",
"\n",
"[D13] Dreesen, Philippe.\n",
"*Back to the Roots: Polynomial System Solving Using Linear Algebra*\n",
"Ph.D. thesis (2013)"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using LinearAlgebra\n",
"using TypedPolynomials\n",
"using MacaulayMatrix\n",
"using JuMP\n",
"using MultivariateMoments"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"Consider the system given in [D13, Example 6.14] which corresponds to `Systems/dreesen10` of [macaulaylab](http://www.macaulaylab.net/):"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "4-element Vector{MultivariatePolynomials.VectorPolynomial{Float64, MultivariatePolynomials.Term{Float64, TypedPolynomials.Monomial{(x₁, x₂, x₃, x₄), 4}}}}:\n -1.0 + x₂ + x₁\n x₂x₄ + x₁x₃\n -0.6666666666666666 + x₂x₄² + x₁x₃²\n x₂x₄³ + x₁x₃³"
},
"metadata": {},
"execution_count": 2
}
],
"cell_type": "code",
"source": [
"@polyvar x[1:4]\n",
"system = [\n",
" x[1] + x[2] - 1,\n",
" x[1] * x[3] + x[2] * x[4],\n",
" x[1] * x[3]^2 + x[2] * x[4]^2 - 2/3,\n",
" x[1] * x[3]^3 + x[2] * x[4]^3,\n",
"]"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"With the classical MacaulayMatrix approach, the nullity increases by 2 at every degree\n",
"because of the positive dimensional solution set at infinity."
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ Info: Added 56 rows to complete columns up to degree 4\n",
"[ Info: Nullspace of dimensions (65, 15) computed from Macaulay matrix of dimension (56, 65) in 0.000884172 seconds.\n",
"[ Info: Added 69 rows to complete columns up to degree 5\n",
"[ Info: Nullspace of dimensions (120, 17) computed from Macaulay matrix of dimension (125, 120) in 0.001944609 seconds.\n",
"[ Info: Added 121 rows to complete columns up to degree 6\n",
"[ Info: Nullspace of dimensions (203, 18) computed from Macaulay matrix of dimension (246, 203) in 0.007104496 seconds.\n",
"[ Info: Found 2 real solutions\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "2-element Vector{Vector{Float64}}:\n [0.4999999999999999, 0.5000000000000001, 0.8164965809277263, -0.8164965809277261]\n [0.5000000000000001, 0.5, -0.816496580927726, 0.8164965809277256]"
},
"metadata": {},
"execution_count": 3
}
],
"cell_type": "code",
"source": [
"solve_system(system, column_maxdegree = 6)"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"Let's compute the moment matrix of degree 4 from the nullspace of\n",
"the Macaulay matrix of degree 5 using the Clarabel SDP solver:"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ Info: Nullspace of dimensions (120, 17) computed from Macaulay matrix of dimension (125, 120) in 0.001968263 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 = 17\n",
" constraints = 121\n",
" nnz(P) = 0\n",
" nnz(A) = 2025\n",
" cones (total) = 2\n",
" : Zero = 1, numel = 1\n",
" : PSDTriangle = 1, numel = 120\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.01e-01 3.42e+00 1.00e+00 3.22e+00 ------ \n",
" 1 0.0000e+00 -1.1323e-03 1.13e-03 5.81e-02 4.66e-01 1.19e-01 4.64e-01 8.61e-01 \n",
" 2 0.0000e+00 2.8687e-04 2.87e-04 1.56e-03 1.23e-02 3.48e-03 1.27e-02 9.74e-01 \n",
" 3 0.0000e+00 2.8697e-06 2.87e-06 1.56e-05 1.23e-04 3.48e-05 1.28e-04 9.90e-01 \n",
" 4 0.0000e+00 2.8697e-08 2.87e-08 1.56e-07 1.23e-06 3.48e-07 1.28e-06 9.90e-01 \n",
" 5 0.0000e+00 2.8697e-10 2.87e-10 1.56e-09 1.23e-08 3.48e-09 1.28e-08 9.90e-01 \n",
" 6 0.0000e+00 2.8697e-10 2.87e-10 1.56e-09 1.23e-08 3.48e-09 1.28e-08 0.00e+00 \n",
"---------------------------------------------------------------------------------------------\n",
"Terminated with status = solved (reduced accuracy)\n",
"solve time = 12.2ms\n",
"[ Info: Terminated with ALMOST_OPTIMAL (ALMOST_SOLVED) in 0.012238786000000001 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 : -2.86973e-10\n",
"│ \n",
"│ * Work counters\n",
"│ Solve time (sec) : 1.22388e-02\n",
"│ Barrier iterations : 6\n",
"└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:25\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "MomentMatrix with row/column basis:\n MonomialBasis([1, x[4], x[3], x[2], x[1], x[4]^2, x[3]*x[4], x[3]^2, x[2]*x[4], x[2]*x[3], x[2]^2, x[1]*x[4], x[1]*x[3], x[1]*x[2], x[1]^2])\nAnd entries in a 15×15 MultivariateMoments.SymMatrix{Float64}:\n 1.0 3.725744452537372e-10 … 0.2500000002434462\n 3.725744452537372e-10 0.6666666666666679 9.315270776966145e-11\n -3.7257428739390086e-10 -0.6666666666666662 -9.313425725077096e-11\n 0.5000000000000245 1.8628763028688544e-10 0.12499999987823852\n 0.499999999999976 1.8628704828715925e-10 0.1250000003652074\n 0.6666666666666679 2.4838786583103456e-10 … 0.16666666650438633\n -0.6666666666666662 -2.4838329223259015e-10 -0.16666666682896436\n 0.6666666666666669 2.4838344055144734e-10 0.1666666665043526\n 1.8628763028688544e-10 0.333333333333317 4.6547795999607455e-11\n -1.8628716624835562e-10 -0.3333333333333493 -4.6587067537018356e-11\n 0.25000000024349434 9.315274246413097e-11 … 0.2555649744273872\n 1.8628704828715925e-10 0.33333333333335063 4.660461599970134e-11\n -1.8628748023330477e-10 -0.333333333333317 -4.654744515178444e-11\n 0.2499999997565302 9.313424337298315e-11 -0.13056497454914884\n 0.2500000002434462 9.315270776966145e-11 0.2555649749143563"
},
"metadata": {},
"execution_count": 4
}
],
"cell_type": "code",
"source": [
"import Clarabel\n",
"ν = moment_matrix(system, Clarabel.Optimizer, 5)"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"The nullspace of this moment matrix is a Macaulay matrix"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "0×0 Macaulay matrix for polynomials:\n 5.897711921664038e-14 + 0.9999999999999994*x[4] + x[3] - 1.4682525325425826e-14*x[4]^2\n -0.5000000000001729 + 3.2596772937398524e-16*x[4] + x[2] + 1.0484472896961149e-13*x[4]^2\n -0.49999999999983163 + 6.156031452731937e-17*x[4] + x[1] - 1.0394986062844766e-13*x[4]^2\n 0.666666666666666 + 5.268577352792162e-16*x[4] - 2.2195756829484477e-16*x[4]^2 + x[3]*x[4]\n -1.536525535226813e-14 - 0.49999999999999983*x[4] + 9.01468757230224e-15*x[4]^2 + x[2]*x[4]\n 1.6090974574336504e-14 - 0.5000000000000008*x[4] - 1.0732798282105085e-14*x[4]^2 + x[1]*x[4]\nThe row shifts are:\nMonomialBasis([])\nThe column basis is:\nMonomialBasis([])\n"
},
"metadata": {},
"execution_count": 5
}
],
"cell_type": "code",
"source": [
"M = nullspace(ν, ShiftNullspace())"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"Each coefficient is close to a rational number, after some clean up,\n",
"we get the following system:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "6-element Vector{MultivariatePolynomials.VectorPolynomial{Float64, MultivariatePolynomials.Term{Float64, TypedPolynomials.Monomial{(x₁, x₂, x₃, x₄), 4}}}}:\n x₄ + x₃\n -1.0 + 2.0x₂\n -1.0 + 2.0x₁\n 2.0 + 3.0x₃x₄\n -x₄ + 2.0x₂x₄\n -x₄ + 2.0x₁x₄"
},
"metadata": {},
"execution_count": 6
}
],
"cell_type": "code",
"source": [
"real_system = clean(M; tol = 1e-6).polynomials"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"We have the equations `2x[1] = 1` and `2x[2] = 1` from which we get\n",
"`x[1] = 1/2` and `x[2] = 1/2`.\n",
"We then have the equations `x[3] = -x[4]` and `3x[3] * x[4] = -2` from\n",
"which we get `x[3] = ±√(2/3)` and `x[4] = ∓√(2/3)`.\n",
"Notice that the complex solutions are not solutions of the system anymore.\n",
"The Macaulay solver indeed find these real solutions direction at degree 2."
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ Info: Added 18 rows to complete columns up to degree 2\n",
"[ Info: Nullspace of dimensions (15, 2) computed from Macaulay matrix of dimension (18, 15) in 9.5116e-5 seconds.\n",
"[ Info: Found 2 real solutions\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": "2-element Vector{Vector{Float64}}:\n [0.4999999999999999, 0.5000000000000001, 0.8164965809277269, -0.8164965809277265]\n [0.5, 0.4999999999999998, -0.8164965809277267, 0.8164965809277263]"
},
"metadata": {},
"execution_count": 7
}
],
"cell_type": "code",
"source": [
"solve_system(real_system, column_maxdegree = 2)"
],
"metadata": {},
"execution_count": 7
},
{
"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
}