ACSVHomotopy¶

Setup¶

Import and activate the locally stored ACSVHomotopy package. The package reexports @polyvar along with functionality from HomotopyContinuation for working with polynomials.

In [1]:
# If you don't already have HomotopyContinuation.jl installed this might take a while
using Pkg
Pkg.activate(".")
using ACSVHomotopy

  Activating project at `~/Desktop/Projects/ACSVHomotopy`

Example Usage¶

Consider first the rational function $$ F(x,y,z) = \frac{1}{1-(1+z)(x+y-xy)},$$ whose main diagonal is related to one of Apéry's irrationality proofs.

In [2]:
# Define the polynomial
@polyvar x y z
H = 1-(1+z)*(x+y-x*y)
Out[2]:
$$ xyz + xy - xz - yz - x - y + 1 $$
In [3]:
# critical points obtained by solving the smooth critical point system.
# Note: This step is not necessary to use our methods.
function get_crits(h)
    vars = variables(h)
    n = length(vars)
    @polyvar λ
    sys = System(
        [h; vars .* differentiate(h, vars) - λ], 
        variables=[vars; λ])
    return [sol[1:end-1] for sol in solutions(solve(sys))]
end

# In this case, we get three critical points. But which (if any) is minimal?
crits = get_crits(H)
Out[3]:
3-element Vector{Vector{ComplexF64}}:
 [0.9999999999999999 + 0.0im, 0.9999999999999999 + 4.591774807899561e-41im, 6.293253624448641e-17 - 2.2958874039497803e-41im]
 [2.618033988749895 + 1.401298464324817e-45im, 2.618033988749895 + 1.7516230804060213e-45im, -1.618033988749895 + 8.758115402030107e-47im]
 [0.38196601125010504 - 1.7516230804060213e-46im, 0.38196601125010504 - 9.85287982728387e-47im, 0.6180339887498951 + 3.503246160812043e-46im]
In [4]:
# Knowing that F is combinatorial,
# we can get the minimal critical point as follows
min_cp = find_min_crits_comb(H)
Out[4]:
1-element Vector{Vector{ComplexF64}}:
 [0.3819660112501051 + 5.877471754111438e-39im, 0.3819660112501051 + 5.877471754111438e-39im, 0.6180339887498949 + 2.350988701644575e-38im]
In [5]:
# Knowing the minimal critical point, we can get the leading asymptotic term
# The exponential growth implies irrationality of zeta(2)
leading_asymptotics(1,H,min_cp)
Out[5]:
"(0.09016994374947422 + 6.20501202900196e-39im)^(-n) n^(-1.0) ((0.4767251395617113 - 9.688970898543992e-38im))"
In [6]:
# It's not obvious from the definition that F is combinatorial.
# If we can't assume this, then we can use the find_min_crits command.
# WARNING: This cell takes 10 - 20 MINUTES to run
# Set show_progress=true to see the progess of the homotopy solver
find_min_crits(H; show_progress=true)

Tracking 36720 paths... 100%|███████████████████████████| Time: 0:07:31
  # paths tracked:                  36720
  # non-singular solutions (real):  16 (0)
  # singular endpoints (real):      60 (0)
  # total solutions (real):         76 (0)
Out[6]:
1-element Vector{Vector{ComplexF64}}:
 [0.38196601125010515 + 0.0im, 0.38196601125010515 + 1.8367099231598242e-40im, 0.6180339887498947 - 7.346839692639297e-40im]
In [7]:
# The general code is very expensive.
# If we want to use the numeric approximation approach of the paper,
# we can obtain the same point in a less rigorous manner much faster.
find_min_crits(H, approx_crit=true)
Out[7]:
1-element Vector{Vector{ComplexF64}}:
 [0.3819660112501052 - 2.311115933264683e-33im, 0.3819660112501052 + 0.0im, 0.6180339887498947 + 9.14816723583937e-34im]
In [8]:
# We can also use the monodromy approach from the paper
find_min_crits(H, monodromy=true)
Out[8]:
1-element Vector{Vector{ComplexF64}}:
 [0.38196601125010515 - 6.244813738743402e-39im, 0.38196601125010515 - 4.775445800215543e-39im, 0.6180339887498947 + 1.4693679385278594e-38im]

Benchmarking¶

In order to properly bench mark the functions, the each cell should be run twice. The first run of the cell with compile and run the function, the @time macro will measure the compilation time plus the running time. On the second run of the cell, the function is compiled and therefore only the running time is measured.

In [43]:
# Define polynomial variables for examples
@polyvar w x y z;
In [19]:
# Binomial coefficients
H = 1 - x - y
@time find_min_crits_comb(H);
@time find_min_crits(H);
@time find_min_crits(H; approx_crit=true);
@time find_min_crits(H; monodromy=true);

  0.005278 seconds (16.86 k allocations: 893.383 KiB)
  0.040406 seconds (113.55 k allocations: 5.896 MiB)
  0.019672 seconds (74.59 k allocations: 3.543 MiB)
  2.294933 seconds (12.08 M allocations: 523.298 MiB, 6.37% gc time)
In [20]:
# Melczer Salvy example of two positive real critical points
H = (1-x-y)*(20-x-40*y)-1 
@time find_min_crits_comb(H);
@time find_min_crits(H);
@time find_min_crits(H; approx_crit=true);
@time find_min_crits(H; monodromy=true);

  0.028840 seconds (35.90 k allocations: 1.691 MiB)
  4.137446 seconds (503.76 k allocations: 22.608 MiB)
  0.336834 seconds (615.59 k allocations: 27.538 MiB)
  2.837442 seconds (12.14 M allocations: 496.231 MiB, 4.03% gc time)
In [25]:
# Related to asymptotics of sqrt(1-z)
H = 1 - x*y - x*y^2 - 2*x^2*y aa
@time find_min_crits_comb(H);
@time find_min_crits(H);
@time find_min_crits(H; approx_crit=true);
@time find_min_crits(H; monodromy=true);

  0.010525 seconds (26.59 k allocations: 1.265 MiB)
 29.567463 seconds (1.75 M allocations: 62.426 MiB, 0.10% gc time)
  0.717018 seconds (665.70 k allocations: 29.158 MiB)
 14.880803 seconds (31.64 M allocations: 1.217 GiB, 1.89% gc time)
In [22]:
# Gillis-Reznick-Zeilberger function
H = 1 - (x+y+z) + 5*x*y*z
@time find_min_crits(H);
@time find_min_crits(H; approx_crit=true);
@time find_min_crits(H; monodromy=true);

Computing mixed cells... 1499 	 Time: 0:00:00
  mixed_volume:  13068

236.717911 seconds (9.71 M allocations: 248.729 MiB, 0.03% gc time)
  3.582242 seconds (953.89 k allocations: 39.776 MiB)
  3.846443 seconds (4.12 M allocations: 168.022 MiB, 1.85% gc time)
In [23]:
# Apéry zeta(2)
H = 1-(1+z)*(x+y-x*y)
@time find_min_crits(H);
@time find_min_crits(H; approx_crit=true);
@time find_min_crits(H; monodromy=true);

Computing mixed cells... 2333 	 Time: 0:00:02
  mixed_volume:  36720

669.286671 seconds (21.70 M allocations: 588.746 MiB, 0.02% gc time)
  3.812803 seconds (685.75 k allocations: 29.921 MiB)
  8.568219 seconds (45.00 M allocations: 1.871 GiB, 5.10% gc time)
In [24]:
# Random poly of degree 4
H = 1 - (72*x^3*z + 97*y*z^3 + 53*x*z^2 + 47*x*y + 39*z^2 + 71*x)
@time find_min_crits_comb(H);
@time find_min_crits(H; approx_crit=true);
@time find_min_crits(H; monodromy=true);

  0.086118 seconds (70.67 k allocations: 3.302 MiB)
189.421242 seconds (13.67 M allocations: 522.053 MiB, 0.04% gc time)
583.118787 seconds (682.06 M allocations: 24.238 GiB, 0.90% gc time)
In [52]:
# 2D lattice path example
H = 1 - z*x^2*y - z*x*y^2 - z*x - z*y
@time find_min_crits_comb(H);
@time find_min_crits(H; approx_crit=true);
@time find_min_crits(H; monodromy=true);

  0.032641 seconds (51.61 k allocations: 2.303 MiB)
 15.298483 seconds (1.30 M allocations: 52.189 MiB, 0.18% gc time)
 31.895528 seconds (174.84 M allocations: 6.555 GiB, 4.73% gc time)
In [27]:
# 3D lattice path example
H = 1 - w*x^2*y*z - w*x*y^2*z - w*x*y*z^2 - w*x*y - w*x*z - w*y*z
@time find_min_crits_comb(H);

  0.078442 seconds (69.97 k allocations: 3.193 MiB)