Faster evaluation

[PBrdng]

Summary

This PR was written by Codex.

This PR substantially accelerates repeated gradient and Hessian evaluations of routing functions associated with projected algebraic hypersurfaces.

The main improvements are:

• caching and moving pseudo-witness fibers between nearby parameter values;
• replacing forward linear solves with an adjoint solve for gradient evaluation;
• evaluating Hessian contractions directly instead of constructing full derivative tensors;
• using the available Julia threads for path tracking;
• implementing a first-order Taylor update for routing gradients;
• adding cache controls and diagnostics.

Mathematical formulation

Let the square system defining a fiber point be

$$G(z;p,b)=F\left(p+\frac{b}{s},u\right)=0,$$

where $z=(s,u)$, $p$ is the evaluation point, and $b$ is the slicing direction. Write

$$J=G_z.$$

For each regular solution $z=z(p,b)$, implicit differentiation gives

$$Jz_{p_a}+G_{p_a}=0,\qquad Jz_{b_b}+G_{b_b}=0.$$

Thus, the forward sensitivities satisfy

$$z_{p_a}=-J^{-1}G_{p_a},\qquad z_{b_b}=-J^{-1}G_{b_b}.$$

The routing function depends on the first coordinate $s$ of every fiber solution.

Adjoint gradient evaluation

The previous implementation recovered derivatives of $s$ by solving the forward systems for every parameter direction. For the gradient, however, only the first component of each sensitivity vector is needed.

Let $e_1$ denote the first coordinate vector and introduce the adjoint variable $\lambda$, defined by

$$J^T\lambda=e_1.$$

For any parameter direction $q$, we then have

$$\frac{\partial s}{\partial q}=e_1^Tz_q=-e_1^TJ^{-1}G_q=-\lambda^TG_q.$$

In particular,

$$-\nabla_b s=G_b^T\lambda.$$

Therefore, the gradient contribution of one fiber point can be computed from the single transposed linear system

$$J^T\lambda=e_1,$$

instead of solving one forward system for every component of $b$. This reduces the number of linear right-hand sides required for a gradient evaluation from $k$ to $1$.

Contracted Hessian evaluation

Differentiating the implicit system a second time gives

$$Jz_{p_ab_b}+R_{ab}=0,$$

where

$$R_{ab}=G_{p_ab_b}+G_{zp_a}z_{b_b}+G_{zb_b}z_{p_a}+G_{zz}[z_{p_a},z_{b_b}].$$

Here, $G_{zz}[z_{p_a},z_{b_b}]$ denotes the contraction of the second derivative of $G$ with the two first-order sensitivity vectors.

Consequently,

$$z_{p_ab_b}=-J^{-1}R_{ab}.$$

Only the first component of this vector enters the Hessian of the routing function. Using the same adjoint variable, we obtain

$$\frac{\partial^2s}{\partial p_a,\partial b_b}=e_1^Tz_{p_ab_b}=-e_1^TJ^{-1}R_{ab}=-\lambda^TR_{ab}.$$

Hence, the required contribution is obtained directly as

$$-\frac{\partial^2s}{\partial p_a,\partial b_b}=\lambda^TR_{ab}.$$

The new implementation evaluates these scalar contractions directly. It no longer constructs the complete second-order sensitivity vectors $z_{p_ab_b}$, nor does it solve $k^2$ additional linear systems for them.

The first-order sensitivities appearing in $R_{ab}$ are still computed together using the block linear system

$$J\begin{bmatrix}z_p&z_b\end{bmatrix}=-\begin{bmatrix}G_p&G_b\end{bmatrix}.$$

After this solve, the Hessian entries are evaluated through the compiled contractions

$$\lambda^TR_{ab}.$$

Compiled derivative systems

The following symbolic expressions are now compiled once when the gradient cache is constructed:

• $G_b^T\lambda$, used for adjoint gradient evaluation;
• $\lambda^TR_{ab}$, used for contracted Hessian evaluation.

This avoids repeatedly assembling large derivative tensors and contracting them in Julia during every routing-function evaluation.

The previous runtime representations of the full tensors $G_{zz}$, $G_{zp}$, $G_{zb}$, and $G_{pb}$, together with their temporary storage and tensor-unpacking logic, are no longer required.

Moving-fiber cache

Repeated evaluations during monodromy usually occur at nearby parameter values. The implementation now caches the most recently computed pseudo-witness fiber and transports it to the next parameter value instead of reconstructing every fiber from the original witness set.

Cache updates are transactional: a transported fiber replaces the cached fiber only after the tracking operation succeeds. If transport from the cached fiber fails, evaluation falls back to the original pseudo-witness construction.

The cache supports:

• exact reuse when the same parameter is evaluated again;
• transport from the most recently accepted fiber;
• fallback to the original pseudo-witness construction;
• cache-hit, cache-miss, and fallback statistics;
• an option to disable moving-fiber reuse.

Taylor prediction

The first-order Taylor update for RoutingGradient now uses

$$\nabla r(p+\Delta p)\approx\nabla r(p)+\mathrm{Hess}r(p),\Delta p.$$

This provides monodromy tracking with a meaningful local prediction for the routing gradient instead of reusing a zeroth-order approximation.

Parallel path tracking

Pseudo-witness transport now uses the available Julia threads. Mutable tracking state is copied for each tracker so that concurrently tracked paths do not share unsafe state.

Correctness

The optimized evaluator was checked against the known symbolic discriminant of the cubic test problem,

$$\Delta(a,b)=4a^3-a^2b^2-18ab+4b^3+27.$$

The pseudo-witness-set gradient and Hessian agree with the symbolic gradient and Hessian at multiple real and complex evaluation points.

The observed errors were below:

• $10^{-11}$ for the gradient;
• $10^{-10}$ for the Hessian.

Performance

The benchmarks exclude object construction and report the median of five runs.

Avoidance example

This example has $k=4$, degree $4$, and uses 300 evaluations.

Evaluation	main	This PR	Speedup
Gradient	46.72 ms	6.53 ms	$7.16\times$
Gradient and Hessian	51.37 ms	9.93 ms	$5.17\times$

The corresponding per-evaluation times are:

Evaluation	main	This PR
Gradient	155.7 μs	21.8 μs
Gradient and Hessian	171.2 μs	33.1 μs

3-RPR example

This example has $k=3$, degree $12$, and uses 100 evaluations.

Evaluation	main	This PR	Speedup
Gradient	45.52 ms	7.26 ms	$6.27\times$
Gradient and Hessian	49.13 ms	9.23 ms	$5.33\times$

The corresponding per-evaluation times are:

Evaluation	main	This PR
Gradient	455.2 μs	72.6 μs
Gradient and Hessian	491.3 μs	92.3 μs

These measurements cover the repeated evaluator workload used by monodromy, rather than a complete end-to-end monodromy run.

Setup-time tradeoff

Compiling the additional adjoint and contracted derivative systems slightly increases GradientCache construction time:

Example	main	This PR	Increase
Avoidance	0.88 s	1.02 s	0.14 s
3-RPR	0.81 s	1.09 s	0.29 s

The complete pseudo-witness-set setup for these examples takes approximately 30 seconds. The additional compilation cost is therefore below $1%$ of the total setup time.

The extra setup cost is recovered after approximately:

• 1,030 gradient-and-Hessian evaluations for the avoidance example;
• 720 gradient-and-Hessian evaluations for the 3-RPR example.

This tradeoff is negligible for monodromy computations, where the evaluator is called many times.

[PBrdng]

Note: I neither checked the new code nor the math behind it. This should be done before merging.

[PBrdng]

@oskarhenriksson @dave-k-johnson @johndcobb

[johndcobb]

The math looks good to me! As before, we should expect the code to act incorrectly at singular psuedowitness set intersections, but this will not happen generically (i.e. see the #13). My own codex only found a few things, mostly related to parallelization:

• ntrackers did not actually control how many workers were present because it was overwritten by HomotopyContinuation.
• fiber_tracking_stats does not aggregate total tracking activity from each worker
  I'll push the recommended changes here.