Experiment 2: Intersecting curves on a Hirzebruch surface

In this notebook we use the Cox homotopy to intersect curves on the Hirzebruch surface \(\mathscr{H}_2\).

We first load our routines.

include("CoxHomotopy.jl")

trackPathRandom (generic function with 1 method)

We generate two bivariate equations.

n = 2;
@polyvar t[1:n];

f̂ = [1 + t[1] + t[2] + t[1]*t[2] + t[1]^2*t[2] + t[1]^3*t[2]; 1 + t[1] + 2*t[2] + 3*t[1]*t[2] + 4*t[1]^2*t[2] + 1*t[1]^3*t[2]];
f̂ = [1 + t[1] + t[2]*(t[1]-1)*(1+sum(map(i->randn()*t[1]^i,1:2))); 1 + t[1] + t[2]*(t[1]-1)*(1+sum(map(i->randn()*t[1]^i,1:2)))]

2-element Array{DynamicPolynomials.Polynomial{true,Float64},1}:
 -0.36515775180019705t₁³t₂ + 0.3626534954403268t₁²t₂ + 1.0025042563598703t₁t₂ + t₁ - t₂ + 1.0
 0.4374666674815006t₁³t₂ + 0.0017765144844340552t₁²t₂ + 0.5607568180340654t₁t₂ + t₁ - t₂ + 1.0

We compute the polytope information.

P, Fᵀ, a = computePolytope(f̂);

We solve the system using the Cox homotopy.

@time homSols, f, toricSols, Fᵀ, x̂, K, B, data = coxHomotopy(f̂; fixSlice = true, details = true);

  0.053850 seconds (67.02 k allocations: 4.872 MiB)

We compute the residuals of the toric solutions and the solutions in their homogeneous coordinates.

display(get_residual(f̂,toricSols,t))
display(get_residual(f,homSols,variables(f)))

4-element Array{Float64,1}:
 2.1506011421870678e-16
 1.0816808319406293e-15
 2.5502337173693015e-16
 1.6503625390396785e-16

4-element Array{Float64,1}:
 2.2484317241831537e-16
 1.7830125319256865e-15
 9.223413031523224e-17
 2.08199467410775e-16

We plot the magnitude of the homogeneous coordinates of the solutions. Each solution corresponds to a column in the table and each row represents a Cox coordinate. Dark (blue) colors indicate small absolute values. The plot shows that only one solution is in the torus. The others are on the divisors \(D_1, D_2\) and \(D_3\).

using PyPlot
imshow(log10.(abs.(hcat(homSols...))))

PyObject <matplotlib.image.AxesImage object at 0x179c7b790>

To illustrate the orthogonal slicing method from Section 4.4 in our paper, we generate some equations of a higher degree on the same toric variety.

Pbig = @pm polytope.Polytope(FACETS = hcat(5*a, Fᵀ))
exps = getLatticePoints(Pbig,t)
mons = [prod(t.^exps[ℓ,:]) for ℓ = 1:size(exps,1)]
cmatrix = randn(2,length(mons))
f̂ = [cmatrix[ℓ,:]'*mons for ℓ = 1:2];

We compute a generic start system and track the solution with index “solidx” (chosen randomly) using 10 different random, fixed linear spaces.

ĝ,pt = computeStartSystem(f̂)
homStartsols, p, fullHomStartsols = solveGeneric(ĝ, pt, Fᵀ)
γ = randn(ComplexF64)
(k,n) = size(Fᵀ)
@polyvar x[1:4]
f, αf = homogenize(f̂, Fᵀ,x)
g, αg = homogenize(ĝ, Fᵀ,x)
Δt = 0.001 #discretization step for naive homotopy algorithm

numexps = 10
solidx = 37
κᵣ = fill([],numexps)
torsols = fill([],numexps)
linslices = fill([],numexps)
for i = 1:numexps
    println("random slice number " *string(i))
    zᵣ, κᵣ[i], linslices[i] = trackPathRandom(f, g, x, fullHomStartsols[solidx], n, k, Δt,γ)
    if i == 1
        semilogy(1:-Δt:Δt,κᵣ[i],color = "#14edc5", label="random slice")
    else
        semilogy(1:-Δt:Δt,κᵣ[i],color = "#14edc5")
    end
    println("residual: " * string(get_residual(f,[zᵣ],x)))
    torsols[i] = dehomogenizeSolutions([zᵣ], Fᵀ)
    println("computed solution: " * string(torsols[i]))
end
legend()
xlabel("τ")
ylabel("κ")

meanκ = 10.0.^(sum([log10.(ℓ) for ℓ ∈ κᵣ])/numexps);

Tracking 100 paths... 100%|█████████████████████████████| Time: 0:00:03
  # paths tracked:                  100
  # non-singular solutions (real):  100 (0)
  # singular endpoints (real):      0 (0)
  # total solutions (real):         100 (0)
random slice number 1
residual: [6.998545843998793e-16]
computed solution: Any[Complex{Float64}[-0.6117540699006732 + 1.1173490185265766im, 0.10427879561967845 - 0.48453334544885374im]]
random slice number 2
residual: [7.913460749832348e-16]
computed solution: Any[Complex{Float64}[-0.6117540699006729 + 1.1173490185265764im, 0.10427879561967861 - 0.48453334544885385im]]
random slice number 3
residual: [1.272576505061576e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006734 + 1.1173490185265766im, 0.1042787956196785 - 0.4845333454488539im]]
random slice number 4
residual: [6.932547917400776e-16]
computed solution: Any[Complex{Float64}[-0.6117540699006732 + 1.1173490185265766im, 0.10427879561967854 - 0.48453334544885385im]]
random slice number 5
residual: [1.0122187521958066e-17]
computed solution: Any[Complex{Float64}[-0.6117540699006729 + 1.1173490185265766im, 0.10427879561967854 - 0.48453334544885374im]]
random slice number 6
residual: [1.338705067677462e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006732 + 1.1173490185265766im, 0.10427879561967861 - 0.4845333454488538im]]
random slice number 7
residual: [2.518462780271932e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006734 + 1.1173490185265766im, 0.10427879561967832 - 0.4845333454488536im]]
random slice number 8
residual: [2.7865124190615667e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006732 + 1.1173490185265764im, 0.10427879561967845 - 0.4845333454488537im]]
random slice number 9
residual: [2.333360640446741e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006736 + 1.1173490185265766im, 0.1042787956196782 - 0.4845333454488536im]]
random slice number 10

residual: [1.4710097636599483e-16]
computed solution: Any[Complex{Float64}[-0.6117540699006733 + 1.1173490185265764im, 0.10427879561967829 - 0.48453334544885374im]]

We track the same path using orthogonal slicing.

NF = smith(Fᵀ);
PP = round.(Int,inv(NF.S));
@time zₒ, κₒ = trackPathOrthogonal(f, g, x, fullHomStartsols[solidx], PP, n, k, Δt, γ);

torsolorth = dehomogenizeSolutions([zₒ], Fᵀ);
println("computed solution: " * string(torsolorth))

  5.568242 seconds (54.64 M allocations: 4.059 GiB, 18.42% gc time)
computed solution: Array{Complex{Float64},1}[[-0.6117519112530363 + 1.1173471479487864im, 0.10427739646186379 - 0.4845344267016235im]]

We plot the results.

semilogy(1:-Δt:Δt,κₒ,color = "#0f87bf", label = "orthogonal slice")
semilogy(1:-Δt:Δt,meanκ, color = "#ff0000", label = "mean random slice")
xlabel("τ")
ylabel("κ")
#legend(["orthogonal slice", "random slice"])
legend()

PyObject <matplotlib.legend.Legend object at 0x179343c40>