# Experiment 3: Solving equations on a weighted projective space

We generate a system of polynomial equations that have solutions of large modulus. Many solvers prematurely truncate these paths and do not compute these solutions. These equations homogenize to equations on a weighted projective space, where we compute all solutions.

We activate the project’s environment.

:

using Pkg
Pkg.activate(".")

 Activating environment at ~/Documents/Projects/cox-homotopies/code_cox_homotopies_fix/Project.toml


We include our methods below.

:

include("CoxHomotopy.jl")

polymake version 4.2
Ewgenij Gawrilow, Michael Joswig, and the polymake team
Technische Universität Berlin, Germany
https://polymake.org

This is free software licensed under GPL; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


:

trackPathRandom (generic function with 1 method)


The following code generates a system that homogenizes to a weighted projective space of the specified weights.

This first block generates the polytope corresponding to our weighted projective space.

:

weights = [1 2 2 2 4];
n = length(weights);
verts = zeros(Int64,n,n+1);
for i = 1:n
verts[i,i] = lcm(weights)/weights[i]
end
P = convexHull(3*verts);
Fᵀ,a = facetRepresentation(P);

polymake: used package ppl
The Parma Polyhedra Library ([[wiki:external_software#PPL]]): A C++ library for convex polyhedra
and other numerical abstractions.
http://www.cs.unipr.it/ppl/



We create a degenerate system $$\hat{f}$$ below. The degeneracy comes from solutions lying on the intersection of divisors $$D_5\cap D_6$$.

:

k = n+1;
@polyvar x[1:k] t[1:n];
latpts = getLatticePoints(P,x);
mons = [prod(t.^latpts[ℓ,:]) for ℓ=1:size(latpts,1)];

M5 = (Fᵀ[5,:]'*latpts')[:];
M6 = (Fᵀ[6,:]'*latpts')[:];
indsmonsonF5 = findall(ℓ -> ℓ == minimum(M5), M5);
indsmonsonF6 = findall(ℓ -> ℓ == minimum(M6), M6);
indsmonsonF56 = intersect(indsmonsonF5,indsmonsonF6);
indsmonsnotonF56 = setdiff(1:length(mons),indsmonsonF56);
g = convert(Array{ComplexF64},rand(1:10,length(indsmonsonF56)))'*mons[indsmonsonF56];
f̂ = [convert(Array{ComplexF64},rand(1:10,length(mons)))'*mons for ℓ=1:n-3];
f̂ = vcat([g + convert(Array{ComplexF64},rand(1:10,length(indsmonsnotonF56)))'*mons[indsmonsnotonF56] for i = 1:3],f̂);


We perturb the degenerate system $$\hat{f}$$ to move all solutions to the torus.

:

ε = 1e-7
for i = 1:3
f̂[i] = f̂[i] + ε*exp.(randn(length(mons))*im)'*mons
end


We solve the system using the Cox homotopy.

:

@time homSols, f, toricSols, homTargetResult = coxHomotopy(f̂, P, Fᵀ, a);

Tracking 7776 paths... 100%|████████████████████████████| Time: 0:05:19
# paths tracked:                  7776
# non-singular solutions (real):  7776 (0)
# singular endpoints (real):      0 (0)
# total solutions (real):         7776 (0)
1332.423907 seconds (263.71 M allocations: 22.131 GiB, 0.75% gc time)

polymake: used package cdd
cddlib
Implementation of the double description method of Motzkin et al.
http://www-oldurls.inf.ethz.ch/personal/fukudak/cdd_home/



Here’s a plot of the magnitude of the computed solutions.

:

using PyPlot
maxcoords = [maximum(abs.(sol)) for sol ∈ toricSols]
semilogy(sort(maxcoords),"bo")
xlabel("solution index");
ylabel("solution norm"); All computed solutions (including the large ones) have a small backward error.

:

residuals = get_residual(f̂,toricSols,t)
semilogy(residuals,"ro") :

1-element Array{PyCall.PyObject,1}:
PyObject <matplotlib.lines.Line2D object at 0x1777fe610>


Due to premature truncation of seemingly diverging paths, the standard polyhedral homotopy misses about 200 solutions.

:

solve(f̂)

Tracking 7776 paths... 100%|████████████████████████████| Time: 0:04:05
# paths tracked:                  7776
# non-singular solutions (real):  7561 (24)
# singular endpoints (real):      0 (0)
# total solutions (real):         7561 (24)

:

Result with 7561 solutions
==========================
• 7776 paths tracked
• 7561 non-singular solutions (24 real)
• random_seed: 0xf28bb5a0
• start_system: :polyhedral


:

length(findall(r->r<1e-10,residuals))

:

7776

[ ]:

nruns = 10
nsols_HC_list = []
for ℓ = 4:16
resℓ = []
for s = 1:nruns
k = n+1;
@polyvar x[1:k] t[1:n];
latpts = getLatticePoints(P,x);
mons = [prod(t.^latpts[ℓ,:]) for ℓ=1:size(latpts,1)];

M5 = (Fᵀ[5,:]'*latpts')[:];
M6 = (Fᵀ[6,:]'*latpts')[:];
indsmonsonF5 = findall(ℓ -> ℓ == minimum(M5), M5);
indsmonsonF6 = findall(ℓ -> ℓ == minimum(M6), M6);
indsmonsonF56 = intersect(indsmonsonF5,indsmonsonF6);
indsmonsnotonF56 = setdiff(1:length(mons),indsmonsonF56);
g = convert(Array{ComplexF64},rand(1:10,length(indsmonsonF56)))'*mons[indsmonsonF56];
f̂ = [convert(Array{ComplexF64},rand(1:10,length(mons)))'*mons for ℓ=1:n-3];
f̂ = vcat([g + convert(Array{ComplexF64},rand(1:10,length(indsmonsnotonF56)))'*mons[indsmonsnotonF56] for i = 1:3],f̂);
ε = 10.0^(-ℓ)
for i = 1:3
f̂[i] = f̂[i] + ε*exp.(randn(length(mons))*im)'*mons
end
#@time homSols, f, toricSols, homTargetResult = coxHomotopy(f̂, P, Fᵀ, a);
R = solve(f̂)
residuals_HC = get_residual(f̂,solutions(R),t)
nsols_HC = length(findall(r->r<1e-10,residuals_HC))
push!(resℓ,nsols_HC)
end
push!(nsols_HC_list,sum(nsols_HC)/nruns)
end
nsols_HC_list

Tracking 7776 paths...  10%|██▉                         |  ETA: 0:03:38
# paths tracked:                  813
# non-singular solutions (real):  812 (0)
# singular endpoints (real):      0 (0)
# total solutions (real):         812 (0)

[ ]: