ProjectedHypersurfaces.jl

This package implements numerical elimination techniques for real algebraic hypersurfaces that arise as the Zariski closure of the projection of a variety with known defining equations.

The key concept from numerical algebraic geometry that underlies the package is that of pseudo-witness sets. For an accessible introduction, as well as some new results on which this package builds, see the recent preprint Elimination Without Eliminating: Computing Complements of Real Hypersurfaces Using Pseudo-Witness Sets by Paul Breiding, John Cobb, Aviva Englander, Nayda Farnsworth, Jon Hauenstein, Oskar Henriksson, David Johnson, Jordy Lopez Garcia, and Deepak Mundayur (2026).

Note

This package is under active early development. Interfaces and behavior may change as the package evolves. In particular, no formal versioned releases are not available yet. For reproducible use, consider pinning to a specific git commit. If you encounter bugs, have questions or want to propose features, feel free to raise an issue or reach out directly to one of the authors.

Installation

You can install the package directly from the github repository as follows:

julia> using Pkg
julia> Pkg.add(url="https://github.com/oskarhenriksson/ProjectedHypersurfaces.jl")

To use the package, make sure that you have activated a Julia environment where the package is added.

You can then load the package in a Julia session by running the following command:

julia> using ProjectedHypersurfaces

Examples of usage

As a case study, suppose that we want to study the complement of the discriminant for the quadratic polynomial

$$f_{a,b}(x)=x^2+ax+b$$

with parameters $a$ and $b$.

We start by setting up the incidence variety $\{(a,b,x)\in ℂ^3\mid f_{a,b}(x)=f′_{a,b}(x)=0\}$ of the discriminant, which we use to form a ProjectedHypersurface that represents the discrimiminant via a pseudo-witness set.

julia> @var a b x;
julia> F = System([x^2 + a * x + b, 2x + a], variables = [a, b, x]);
julia> h = ProjectedHypersurface(F, [a, b])
Projected hypersurface of degree 2 in ambient dimension 2

Degree

We can extract the degree of the hypersurface as follows:

julia> degree(h)
2

Trace test

To verify the completeness of the pseudowitness set (and hence the correctness of the degree), we can run a trace test. Theoretically, this value is zero if and only if the pseudowitness set is complete. Hence, a value close to machine precision is strong evidence (albeit not a certificate) of completeness.

julia> trace_test(h)
1.4101715336057762e-18

Membership test

The pseduo-witness set constitutes a powerful implicit representation of the hypersurface. For instance, we can test membership by moving the pseudowitness line so that it passes through the candidate point, and check if the pseudowitness points converge to the candidate.

julia> contains(h, [2, 1])
true

julia> containts(h, [1, 1])
false

Sampling

By moving around the pseudowitness line, we can easilly obtain a large set of sample points from the hypersurface.

julia> sample_points(h, 10)
10-element Vector{Vector{ComplexF64}}:
 [-0.44365770524121323 + 0.2686256389642566im, 0.031168106377736024 - 0.05958891727591836im]
 [11.842559600352642 - 4.795419017670839im, 29.312543583216335 - 28.395017762715714im]
 [-0.8996583995828127 + 0.8399604833359556im, 0.025962905593483965 - 0.37783875207541584im]
 [12.298560294694239 - 5.3667538620425415im, 30.613134576620286 - 33.00167297955668im]
 [-0.4922927872348897 - 0.2919171871604691im, 0.03928413605095399 + 0.07185436285449809im]
 [11.89119468234632 - 4.234876191546116im, 30.86658365393431 - 25.1788686246541im]
 [0.7673031891362762 - 0.5085351244614289im, 0.08253655281192462 - 0.1951003113935338im]
 [10.631598705975149 - 4.018258254245153im, 24.2211229117708 - 21.36025462805337im]
 [-0.40224428636217896 - 0.29158223834385893im, 0.019195066048350917 + 0.05864364468925615im]
 [11.801146181473607 - 4.2352111403627255im, 30.332509448264137 - 24.990172888413024im]

Interpolation

Based on a large sample, we can attempt to interpolate a defining polynomial for the hypersurface. This works best if the degree of h is low.

julia> interpolate(h)
Interpolation result for projected hypersurface
===============================================
 Smallest singular value: 9.2901e-17
 Ratio of next-smallest to smallest singular value: 1.0873e16
 Residual: 1.1115e-16
-----------------------------------------------
 Variables: a, b
 Polynomial: -4*b + a^2

Irreducible decomposition

We can do a numerical irreducible decomposition of the hypersurfaces. For instance:

julia> @var a b x;
julia> f = b + a*x + 2*a*x^4 + a^2*x^3 + 2*b*x^3 + b^2*x + 2*a*b*x^2 + x^2 + x^5;
julia> F = System([f, differentiate(f, x)], variables=[a, b, x]);
julia> h = ProjectedHypersurface(F, [a, b])
Projected hypersurface of degree 6 in ambient dimension 2

julia> irreducible_components = decompose(h)
2-element Vector{ProjectedHypersurface}:
 Projected hypersurface of degree 4 in ambient dimension 2
 Projected hypersurface of degree 2 in ambient dimension 2

julia> polynomial.(interpolate.(irreducible_components))
2-element Vector{Expression}:
                            -4*b + a^2
 27 - 18*a*b - a^2*b^2 + 4*a^3 + 4*b^3

Evaluation

We can use h to evaluate (up to a constant) the logarithm of the defining polynomial of the discriminant, as well as the gradient and Hessian.

julia> p = [1, 1];

julia> h(p) # the value depends on the direction of the pseudo-witness line
1.5362619674238103

julia> gradient(h, p)
2-element Vector{ComplexF64}:
 -0.6666666666666665 + 4.440892098500626e-16im
  1.3333333333333335 - 2.220446049250313e-16im

julia> hessian(h, p) 
2×2 Matrix{ComplexF64}:
 -1.11111-9.99201e-16im  0.888889+4.44089e-16im
 0.888889+7.77156e-16im  -1.77778+9.71445e-16im

Gradient roadmaps

We use h to form a routing function as follows. (If we don't specify the center c for the denominator, it is chosen randomly.)

julia> r = RoutingFunction(h; c=[13, 2])
Routing function for projected hypersurface
===========================================
 Variables: a, b
 Numerator: Projected hypersurface of degree 2 in ambient dimension 2
 Denominator: (1 + (-13 + a)^2 + (-2 + b)^2)^2

We find the critical points via the critical_points function:

julia> routing_result = critical_points(r)
Routing points result with 4 routing point(s)

julia> pts = routing_points(routing_result)
4-element Vector{Vector{Float64}}:
 [13.040296300414134, 1.993819726256856]
 [3.2168112092392143, 8.082538361382136]
 [-3.9180890683992504, -6.635887940807433]
 [-12.339018441254092, -2.1071368134982262]

Finally, we connect the critical points that belong to the same component of the complement, to obtain a gradient roadmap of the complement of the hypersurface.

julia> roadmap = gradient_roadmap(r, routing_result)
Gradient roadmap of a hypersurface complement
=============================================
• 2 connected component(s)
• 4 routing point(s)
• return_code → :success

The gradient roadmap represents each connected component as a Region object.

julia> connected_comps = regions(roadmap)
2-element Vector{Region}:
 Region 1
========
• 3 routing point(s)
• morse_indices → [0, 1, 0]
• χ → 1
 Region 2
========
• 1 routing point(s)
• morse_indices → [0]
• χ → 1

We can obtain the routing points for each region via the routing_points command:

julia> C1 = connected_comps[1];

julia> routing_points(C1)
3-element Vector{Vector{Float64}}:
 [13.040296300414134, 1.993819726256856]
 [-3.9180890683992433, -6.635887940807435]
 [-12.339018441254096, -2.1071368134982267]

julia> C2 = connected_comps[2];

julia> routing_points(C2)
1-element Vector{Vector{Float64}}:
 [3.2168112092392134, 8.082538361382136]

We see that the first, third and fourth critical points belong to the same connected component, and that the second one belongs to its own component. We can obtain this partition via the partition command:

julia> partition(roadmap)
2-element Vector{Vector{Int64}}:
 [1, 3, 4]
 [2]

The roadmap also allows a membership test for points in the complement by tracing them via gradient flow and determining which region's routing points they converge to.

julia> membership(roadmap, [1, 1])
Region 2
========
• 1 routing point(s)
• morse_indices → [0]
• χ → 1

The resulting roadmap is illustrated by the following picture.

Further examples

The following pictures show gradient roadmaps for other examples of discriminants. Code for computing the roadmaps with ProjectedHypersurfaces.jl can be found in the examples directory.

Dependencies

The code relies on the following Julia packages:

• HomotopyContinuation.jl (for numerical algebraic geometry)
• OrdinaryDiffEq.jl (for gradient flow)
• LightGraphs.jl (for building the connectivity graph).

Statement on AI use

This repository was developed with the assistance of AI tools, including large language models such as Codex and and GitHub Copilot, which have been used to assist with tasks such as code review, memory allocation optimization, documentaiton, and minor code generation for routine tasks. The authors have reviewed, edited, and verified the outputs from these tools, and take full responsibility for the correctness.