# 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 include our methods below.

:

include("CoxHomotopy.jl")

:

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);

Solving for 2 parameters... 100%|███████████████████████| Time: 0:09:31
# parameters solved:  2
# paths tracked:      15552
1202.856065 seconds (293.22 M allocations: 22.613 GiB, 0.78% gc time)

polymake: used package cdd
cddlib
Implementation of the double description method of Motzkin et al.
Copyright by Komei Fukuda.
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 0x16e8dbb80>


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

:

solve(f̂)

Tracking 7776 paths... 100%|████████████████████████████| Time: 0:03:54
# paths tracked:                  7776
# non-singular solutions (real):  7560 (7)
# singular endpoints (real):      0 (0)
# total solutions (real):         7560 (7)

:

Result with 7560 solutions
==========================
• 7776 paths tracked
• 7560 non-singular solutions (7 real)
• random_seed: 0x9c53bcc1
• start_system: :polyhedral