Getting rid of root at infinity
Adapted from: Section~2.2.1 p. 33 of [D13]
[D13] Dreesen, Philippe. Back to the Roots: Polynomial System Solving Using Linear Algebra Ph.D. thesis (2013)
using LinearAlgebra
using TypedPolynomials
using MacaulayMatrix
using JuMP
using MultivariateMomentsConsider the system given in [D13, (2.3)] which corresponds to Systems/dreesen2 of macaulaylab:
@polyvar x y z
system = [
x^2 - x*y + z,
2y^3 - 2x*y^2 - 3x * y,
z^3 - x*y*z - 2,
]3-element Vector{MultivariatePolynomials.VectorPolynomial{Int64, MultivariatePolynomials.Term{Int64, TypedPolynomials.Monomial{(x, y, z), 3}}}}:
z - xy + x²
-3xy + 2y³ - 2xy²
-2 + z³ - xyzWe first try to solve the system:
sols = solve_system(system, column_maxdegree = 6)[ Info: Added 6 rows to complete columns up to degree 3
[ Info: Nullspace of dimensions (14, 8) computed from Macaulay matrix of dimension (6, 14) in 6.1113e-5 seconds.
[ Info: Added 12 rows to complete columns up to degree 4
[ Info: Nullspace of dimensions (33, 15) computed from Macaulay matrix of dimension (18, 33) in 0.000101689 seconds.
[ Info: Added 22 rows to complete columns up to degree 5
[ Info: Nullspace of dimensions (56, 18) computed from Macaulay matrix of dimension (40, 56) in 0.000320595 seconds.
[ Info: Added 35 rows to complete columns up to degree 6
[ Info: Nullspace of dimensions (84, 18) computed from Macaulay matrix of dimension (75, 84) in 0.000838557 seconds.
[ Info: Found 1 real solutionThe real solutions are
sols1-element Vector{Vector{Float64}}:
[3.8420560175113305, 4.995671506251074, 4.432255330406786]Staicase analysis
solver = Iterator(system, MacaulayMatrix.Solver())
step!(solver)[ Info: Added 6 rows to complete columns up to degree 3
[ Info: Nullspace of dimensions (14, 8) computed from Macaulay matrix of dimension (6, 14) in 6.1424e-5 seconds.We can look at the solver as follows:
solverMacaulayMatrix matrix solver. Last iteration considered:
6×14 Macaulay matrix for polynomials:
z - x*y + x^2
-3*x*y + 2*y^3 - 2*x*y^2
-2 + z^3 - x*y*z
The row shifts are:
MonomialBasis([1, z, y, x])
The column basis is:
MonomialBasis([1, z, z^2, y*z, x*z, x*y, x^2, z^3, y^3, x*y*z, x*y^2, x^2*z, x^2*y, x^3])
BorderBasis with independent rows and dependent columns in:
BasisDependence for bases:
Standard:
MonomialBasis([1, z, z^2, y*z, x*z, x*y, z^3, y^3])
Trivial Standard:
MonomialBasis([y, x, y^2, y*z^2, y^2*z, x*z^2, z^4, y*z^3, y^2*z^2, y^3*z, y^4, x*z^3])
Corners:
MonomialBasis([x^2, x*y*z, x*y^2])
Trivial Independent Border:
MonomialBasis([x*y*z^2, x*y^2*z, x*y^3, x^2*z^2])
Dependent Border:
MonomialBasis([x^2*z, x^2*y])
And entries in a 8×5 adjoint(::Matrix{Float64}) with eltype Float64:
4.159065386202991e-16 -1.9999999999999998 … -1.9069420129706838e-16
-1.0 -3.1080979595180984e-16 4.5135339788020976e-17
1.1499638451679124e-17 2.176915759031188e-16 -1.2298933558687206e-16
-1.041204178188655e-16 -3.470680593962184e-16 -1.0
-7.832624135939215e-17 -2.4017185283265312e-17 -4.699574481563529e-16
1.0 -4.3281882871179737e-17 … -1.5000000000000002
4.8138576458356415e-17 1.0 -3.677747436866566e-17
-5.551115123125784e-17 3.8163916471489786e-16 1.0000000000000004
Current status is OPTIMIZE_NOT_CALLED
History of iterations:
1×3 DataFrame
Row │ nullity num_rows num_cols
│ Int64 Int64 Int64
─────┼─────────────────────────────
1 │ 8 6 14We can see in the border dependence that even if x^2 is a corner, the multiplication matrix for x cannot be computed as y^3 (resp. x*z^2, y*z^2) is standard but x * y^3 (resp. x^2 * z^2, x * y * z^2) is indepedent.
using Plots
plot(saturated_dependence(solver))
Let's do another step:
step!(solver)[ Info: Added 12 rows to complete columns up to degree 4
[ Info: Nullspace of dimensions (33, 15) computed from Macaulay matrix of dimension (18, 33) in 0.00013418 seconds.This time, the blocker for computing the multiplication matrix for x are z^4 and y * z^3 among others.
plot(saturated_dependence(solver))
Let's do another step:
step!(solver)[ Info: Added 22 rows to complete columns up to degree 5
[ Info: Nullspace of dimensions (56, 18) computed from Macaulay matrix of dimension (40, 56) in 0.00037123 seconds.Now they are z^5 and y * z^4 among others.
plot(saturated_dependence(solver))
Let's do a last step:
step!(solver)[ Info: Added 35 rows to complete columns up to degree 6
[ Info: Nullspace of dimensions (84, 18) computed from Macaulay matrix of dimension (75, 84) in 0.000836954 seconds.
[ Info: Found 1 real solutionNow we see that the whole border is dependent so the four multiplication matrices can be computed.
plot(saturated_dependence(solver))
In retrospect, we we probably should have expanded in priority towards larger exponents for z.
Saturation
Let's start again by starting with the same first stop.
solver = Iterator(system, MacaulayMatrix.Solver())
step!(solver)[ Info: Added 6 rows to complete columns up to degree 3
[ Info: Nullspace of dimensions (14, 8) computed from Macaulay matrix of dimension (6, 14) in 6.0632e-5 seconds.This time, let's focus on saturating the first non-saturated standard monomials:
step!(solver, FirstStandardNonSaturated(10))[ Info: Added 12 rows to saturate columns `TypedPolynomials.Monomial{(x, y, z), 3}[z, y, x, z², yz, y², xz, xy, z³, yz²]`
[ Info: Nullspace of dimensions (38, 20) computed from Macaulay matrix of dimension (18, 38) in 0.00010251 seconds.We can see that z is now saturated:
plot(saturated_dependence(solver))
Let's saturate the next one:
step!(solver, FirstStandardNonSaturated(10))[ Info: Added 20 rows to saturate columns `TypedPolynomials.Monomial{(x, y, z), 3}[y²z, y³, xz², z⁴, yz³, y²z², y³z, y⁴, z⁵, yz⁴]`
[ Info: Nullspace of dimensions (70, 32) computed from Macaulay matrix of dimension (38, 70) in 0.000368775 seconds.We can see that y went from trivial to saturated:
plot(saturated_dependence(solver))
Let's saturate the next one:
step!(solver, FirstStandardNonSaturated(10))[ Info: Added 20 rows to saturate columns `TypedPolynomials.Monomial{(x, y, z), 3}[y⁵, z⁶, yz⁵, y⁶, z⁷, yz⁶, y⁷, z⁸, yz⁷, y⁸]`
[ Info: Nullspace of dimensions (103, 45) computed from Macaulay matrix of dimension (58, 103) in 0.000711572 seconds.We can see that x went from trivial to saturated:
plot(saturated_dependence(solver))
Let's saturate the next one:
step!(solver, FirstStandardNonSaturated(10))[ Info: Added 20 rows to saturate columns `TypedPolynomials.Monomial{(x, y, z), 3}[z⁹, yz⁸, y⁹, z¹⁰, yz⁹, y¹⁰, z¹¹, yz¹⁰, y¹¹, z¹²]`
[ Info: Nullspace of dimensions (136, 58) computed from Macaulay matrix of dimension (78, 136) in 0.001374643 seconds.We can see that z^2 is now saturated:
plot(saturated_dependence(solver))
We can see that this does not seem to work. By saturating a monomial, we generate all the rows that have a nonzero entry for the corresponding column in the MacaulayMatrix matrix. If the monomial is standard, it means that the column is a linear combination of other columns (which are dependent in the MacaulayMatrix Nullspace) of the MacaulayMatrix matrix. If we saturate these other columns, the linear combination may not work anymore. But we are only saturating the standard monomials, we don't saturate the rows that are dependent in the MacaulayMatrix Nullspace which explains the result we have seen here.
Moment approach
With moment matrix of degree 3:
import Clarabel
solver = Clarabel.Optimizer
M = moment_matrix(system, solver, 3)
atomic_measure(M, 1e-4, ShiftNullspace())[ Info: Nullspace of dimensions (14, 8) computed from Macaulay matrix of dimension (6, 14) in 5.952e-5 seconds.
┌ Warning: Missing monomials TypedPolynomials.Monomial{(x, y, z), 3}[y, x, y²] in Macaulay nullspace, leaving them free to take any value in the moment matrix then.
└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:65
-------------------------------------------------------------
Clarabel.jl v0.9.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 11
constraints = 11
nnz(P) = 0
nnz(A) = 67
cones (total) = 2
: Zero = 1, numel = 1
: PSDTriangle = 1, numel = 10
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 3.88e-01 6.59e-01 1.00e+00 1.24e+00 ------
1 0.0000e+00 -6.5201e-02 6.52e-02 3.19e-02 7.78e-02 7.16e-03 1.25e-01 9.13e-01
2 0.0000e+00 -5.8123e-04 5.81e-04 3.43e-04 9.79e-04 1.89e-04 1.46e-03 9.89e-01
3 0.0000e+00 -5.8069e-06 5.81e-06 3.43e-06 9.79e-06 1.89e-06 1.46e-05 9.90e-01
4 0.0000e+00 -5.8068e-08 5.81e-08 3.43e-08 9.79e-08 1.89e-08 1.46e-07 9.90e-01
5 0.0000e+00 -5.8068e-10 5.81e-10 3.43e-10 9.79e-10 1.89e-10 1.46e-09 9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 976μs
[ Info: Terminated with OPTIMAL (SOLVED) in 0.000976112 seconds.With moment matrix of degree 4:
M = moment_matrix(system, solver, 4)
atomic_measure(M, 1e-4, ShiftNullspace())[ Info: Nullspace of dimensions (33, 15) computed from Macaulay matrix of dimension (18, 33) in 0.000106117 seconds.
┌ Warning: Missing monomials TypedPolynomials.Monomial{(x, y, z), 3}[y², y²z²] in Macaulay nullspace, leaving them free to take any value in the moment matrix then.
└ @ MacaulayMatrix ~/.julia/packages/MacaulayMatrix/2RDTD/src/hankel.jl:65
-------------------------------------------------------------
Clarabel.jl v0.9.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 17
constraints = 56
nnz(P) = 0
nnz(A) = 784
cones (total) = 2
: Zero = 1, numel = 1
: PSDTriangle = 1, numel = 55
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 7.02e-01 1.07e+00 1.00e+00 4.05e+00 ------
1 0.0000e+00 2.1244e+00 2.12e+00 2.49e-01 2.23e-01 2.83e+00 8.73e-01 8.09e-01
2 0.0000e+00 1.0127e+01 1.01e+01 4.85e-02 4.27e-02 1.08e+01 2.18e-01 8.63e-01
3 0.0000e+00 2.9803e+01 2.98e+01 7.30e-03 7.27e-03 3.02e+01 4.34e-02 8.59e-01
4 0.0000e+00 1.4744e+01 1.47e+01 1.30e-03 1.40e-03 1.49e+01 9.54e-03 8.22e-01
5 0.0000e+00 3.0716e+00 3.07e+00 2.61e-04 3.30e-04 3.11e+00 3.23e-03 8.96e-01
6 0.0000e+00 8.9652e-02 8.97e-02 1.35e-05 1.74e-05 9.14e-02 1.89e-04 9.44e-01
7 0.0000e+00 9.2794e-04 9.28e-04 1.36e-07 1.76e-07 9.46e-04 1.95e-06 9.90e-01
8 0.0000e+00 9.2781e-06 9.28e-06 1.36e-09 1.76e-09 9.46e-06 1.95e-08 9.90e-01
9 0.0000e+00 9.2783e-08 9.28e-08 1.35e-11 1.76e-11 9.46e-08 1.95e-10 9.90e-01
10 0.0000e+00 1.0031e-09 1.00e-09 2.78e-09 1.90e-13 1.02e-09 2.09e-12 9.89e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 3.22ms
[ Info: Terminated with OPTIMAL (SOLVED) in 0.003217203 seconds.With moment matrix of degree 5:
M = moment_matrix(system, solver, 5)
atomic_measure(M, 1e-4, ShiftNullspace())Atomic measure on the variables x, y, z with 1 atoms:
at [3.84205516042132, 4.99567052245059, 4.432254383584636] with weight 0.9999998706390161This page was generated using Literate.jl.