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.

[2]:
include("CoxHomotopy.jl")
[2]:
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.

[3]:
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\).

[4]:
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];
 = [convert(Array{ComplexF64},rand(1:10,length(mons)))'*mons for =1:n-3];
 = vcat([g + convert(Array{ComplexF64},rand(1:10,length(indsmonsnotonF56)))'*mons[indsmonsnotonF56] for i = 1:3],);

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

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

We solve the system using the Cox homotopy.

[6]:
@time homSols, f, toricSols, homTargetResult = coxHomotopy(, 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.

[15]:
using PyPlot
maxcoords = [maximum(abs.(sol)) for sol  toricSols]
semilogy(sort(maxcoords),"bo")
xlabel("solution index");
ylabel("solution norm");
../_images/CoxHomotopies_Experiment3_14_0.png

All computed solutions (including the large ones) have a small backward error.

[16]:
residuals = get_residual(,toricSols,t)
semilogy(residuals,"ro")
../_images/CoxHomotopies_Experiment3_16_0.png
[16]:
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.

[8]:
solve()
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)
[8]:
Result with 7560 solutions
==========================
• 7776 paths tracked
• 7560 non-singular solutions (7 real)
• random_seed: 0x9c53bcc1
• start_system: :polyhedral