Jonathan D. Hauenstein

Yu Wang, Emma R. Cobian, Jubilee Lee et al.

Variational inference is an increasingly popular method in statistics and machine learning for approximating probability distributions. We developed LINFA (Library for Inference with Normalizing Flow and Annealing), a Python library for variational inference to accommodate computationally expensive models and difficult-to-sample distributions with dependent parameters. We discuss the theoretical background, capabilities, and performance of LINFA in various benchmarks. LINFA is publicly available on GitHub at https://github.com/desResLab/LINFA.

Jonathan D. Hauenstein

Given a polynomial system f, a fundamental question is to determine if f has real roots. Many algorithms involving the use of infinitesimal deformations have been proposed to answer this question. In this article, we transform an approach of Rouillier, Roy, and Safey El Din, which is based on a classical optimization approach of Seidenberg, to develop a homotopy based approach for computing at least one point on each connected component of a real algebraic set. Examples are presented demonstrating the effectiveness of this parallelizable homotopy based approach.

Dhagash Mehta, Jonathan D. Hauenstein, David J. Wales

Typically, there is no guarantee that a numerical approximation obtained using standard nonlinear equation solvers is indeed an actual solution, meaning that it lies in the quadratic convergence basin. Instead, it may lie only in the linear convergence basin, or even in a chaotic region, and hence not converge to the corresponding stationary point when further optimization is attempted. In some cases, these non-solutions could be misleading. Proving that a numerical approximation will quadratically converge to a stationary point is termed \textit{certification}. In this report, we provide details of how Smale's $α$-theory can be used to certify numerically obtained stationary points of a potential energy landscape, providing a \textit{mathematical proof} that the numerical approximation does indeed correspond to an actual stationary point, independent of the precision employed.

Heather A. Harrington, Dhagash Mehta, Helen M. Byrne et al.

Many systems in biology, physics and engineering can be described by systems of ordinary differential equation containing many parameters. When studying the dynamic behavior of these large, nonlinear systems, it is useful to identify and characterize the steady-state solutions as the model parameters vary, a technically challenging problem in a high-dimensional parameter landscape. Rather than simply determining the number and stability of steady-states at distinct points in parameter space, we decompose the parameter space into finitely many regions, the steady-state solutions being consistent within each distinct region. From a computational algebraic viewpoint, the boundary of these regions is contained in the discriminant locus. We develop global and local numerical algorithms for constructing the discriminant locus and classifying the parameter landscape. We showcase our numerical approaches by applying them to molecular and cell-network models.

Alessandra Bernardi, Noah S. Daleo, Jonathan D. Hauenstein et al.

A computationally challenging classical elimination theory problem is to compute polynomials which vanish on the set of tensors of a given rank. By moving away from computing polynomials via elimination theory to computing pseudowitness sets via numerical elimination theory, we develop computational methods for computing ranks and border ranks of tensors along with decompositions. More generally, we present our approach using joins of any collection of irreducible and nondegenerate projective varieties $X_1,\ldots,X_k\subset\mathbb{P}^N$ defined over $\mathbb{C}$. After computing ranks over $\mathbb{C}$, we also explore computing real ranks. Various examples are included to demonstrate this numerical algebraic geometric approach.

Jonathan D. Hauenstein, Jose Israel Rodriguez, Frank Sottile

The Galois/monodromy group of a family of geometric problems or equations is a subtle invariant that encodes the structure of the solutions. Computing monodromy permutations using numerical algebraic geometry gives information about the group, but can only determine it when it is the full symmetric group. We give numerical methods to compute the Galois group and study it when it is not the full symmetric group. One algorithm computes generators while the other gives information on its structure as a permutation group. We illustrate these algorithms with examples using a Macaulay2 package we are developing that relies upon Bertini to perform monodromy computations.

Emilie Dufresne, Parker B. Edwards, Heather A. Harrington et al.

Topological data analysis (TDA) provides a growing body of tools for computing geometric and topological information about spaces from a finite sample of points. We present a new adaptive algorithm for finding provably dense samples of points on real algebraic varieties given a set of defining polynomials. The algorithm utilizes methods from numerical algebraic geometry to give formal guarantees about the density of the sampling and it also employs geometric heuristics to reduce the size of the sample. As TDA methods consume significant computational resources that scale poorly in the number of sample points, our sampling minimization makes applying TDA methods more feasible. We provide a software package that implements the algorithm and also demonstrate the implementation with several examples.

Tulay Ayyildiz Akoglu, Jonathan D. Hauenstein, Agnes Szanto

This paper is concerned with certifying that a given point is near an exact root of an overdetermined or singular polynomial system with rational coefficients. The difficulty lies in the fact that consistency of overdetermined systems is not a continuous property. Our certification is based on hybrid symbolic-numeric methods to compute the exact "rational univariate representation" (RUR) of a component of the input system from approximate roots. For overdetermined polynomial systems with simple roots, we compute an initial RUR from approximate roots. The accuracy of the RUR is increased via Newton iterations until the exact RUR is found, which we certify using exact arithmetic. Since the RUR is well-constrained, we can use it to certify the given approximate roots using alpha-theory. To certify isolated singular roots, we use a determinantal form of the "isosingular deflation", which adds new polynomials to the original system without introducing new variables. The resulting polynomial system is overdetermined, but the roots are now simple, thereby reducing the problem to the overdetermined case. We prove that our algorithms have complexity that are polynomial in the input plus the output size upon successful convergence, and we use worst case upper bounds for termination when our iteration does not converge to an exact RUR. Examples are included to demonstrate the approach.

Jonathan D. Hauenstein, Viktor Levandovskyy

Smale's alpha-theory certifies that Newton iterations will converge quadratically to a solution of a square system of analytic functions based on the Newton residual and all higher order derivatives at the given point. Shub and Smale presented a bound for the higher order derivatives of a system of polynomial equations based in part on the degrees of the equations. For a given system of polynomial-exponential equations, we consider a related system of polynomial-exponential equations and provide a bound on the higher order derivatives of this related system. This bound yields a complete algorithm for certifying solutions to polynomial-exponential systems, which is implemented in alphaCertified. Examples are presented to demonstrate this certification algorithm.

Danielle A. Brake, Jonathan D. Hauenstein, Frank-Olaf Schreyer et al.

Singular value decompositions of matrices are widely used in numerical linear algebra with many applications. In this paper, we extend the notion of singular value decompositions to finite complexes of real vector spaces. We provide two methods to compute them and present several applications.

Jonathan D. Hauenstein, Frank Sottile

Smale's alpha-theory uses estimates related to the convergence of Newton's method to give criteria implying that Newton iterations will converge quadratically to solutions to a square polynomial system. The program alphaCertified implements algorithms based on alpha-theory to certify solutions to polynomial systems using both exact rational arithmetic and arbitrary precision floating point arithmetic. It also implements an algorithm to certify whether a given point corresponds to a real solution to a real polynomial system, as well as algorithms to heuristically validate solutions to overdetermined systems. Examples are presented to demonstrate the algorithms.

Jonathan D. Hauenstein, Margaret H. Regan

Three aspects of applying homotopy continuation, which is commonly used to solve parameterized systems of polynomial equations, are investigated. First, for parameterized systems which are homogeneous, we investigate options for performing computations on an adaptively chosen affine coordinate patch. Second, for parameterized systems which are overdetermined, we investigate options for adaptively selecting a well-constrained subsystem to restore numerical stability. Finally, since one is typically interested in only computing real solutions for parameterized problems which arise from applications, we investigate a scheme for heuristically identifying solution paths which appear to be ending at nonreal solutions and truncating them. We demonstrate these three aspects on two problems arising in computer vision.

Jonathan D. Hauenstein, Frank Sottile

We present two algorithms that compute the Newton polytope of a polynomial defining a hypersurface H in C^n using numerical computation. The first algorithm assumes that we may only compute values of f - this may occur if f is given as a straight-line program, as a determinant, or as an oracle. The second algorithm assumes that H is represented numerically via a witness set. That is, it computes the Newton polytope of H using only the ability to compute numerical representatives of its intersections with lines. Such witness set representations are readily obtained when H is the image of a map or is a discriminant. We use the second algorithm to compute a face of the Newton polytope of the Lüroth invariant, as well as its restriction to that face.

Daniel K. Molzahn, Matthew Niemerg, Dhagash Mehta et al.

The power flow equations model the steady-state relationship between the power injections and voltage phasors in an electric power system. By separating the real and imaginary components of the voltage phasors, the power flow equations can be formulated as a system of quadratic polynomials. Only the real solutions to these polynomial equations are physically meaningful. This paper focuses on the maximum number of real solutions to the power flow equations. An upper bound on the number of real power flow solutions commonly used in the literature is the maximum number of complex solutions. There exist two- and three-bus systems for which all complex solutions are real. It is an open question whether this is also the case for larger systems. This paper investigates four-bus systems using techniques from numerical algebraic geometry and conjectures a negative answer to this question. In particular, this paper studies lossless, four-bus systems composed of PV buses connected by lines with arbitrary susceptances. Computing the Galois group, which is degenerate, enables conversion of the problem of counting the number of real solutions to the power flow equations into counting the number of positive roots of a univariate sextic polynomial. From this analysis, it is conjectured that the system has at most 16 real solutions, which is strictly less than the maximum number of complex solutions, namely 20. We also provide explicit parameter values where this system has 16 real solutions so that the conjectured upper bound is achievable.

Jonathan D. Hauenstein, Margaret H. Regan

The monodromy group is an invariant for parameterized systems of polynomial equations that encodes structure of the solutions over the parameter space. Since the structure of real solutions over real parameter spaces are of interest in many applications, real monodromy action is investigated here. A naive extension of monodromy action from the complex numbers to the real numbers is shown to be very restrictive. Therefore, we define a real monodromy structure which need not be a group but contains tiered characteristics about the real solutions. This real monodromy structure is applied to an example in kinematics which summarizes all the ways performing loops parameterized by leg lengths can cause a mechanism to change poses.

Shaunak D. Bopardikar, Dhagash Mehta, Jonathan D. Hauenstein

We revisit the static coverage control problem for placement of vehicles with simple motion on the real line, under the assumption that the cost is a polynomial function of the locations of the vehicles. The main contribution of this paper is to demonstrate the use of tools from numerical algebraic geometry, in particular, a numerical polynomial homotopy continuation method that guarantees to find all solutions of polynomial equations, in order to characterize the \emph{global minima} for the coverage control problem. The results are then compared against a classic distributed approach involving the use of Lloyd descent, which is known to converge only to a local minimum under certain technical conditions.

Daniel A. Brake, Jonathan D. Hauenstein, Cynthia Vinzant

Tropical varieties capture combinatorial information about how coordinates of points in a classical variety approach zero or infinity. We present algorithms for computing the rays of a complex and real tropical curve defined by polynomials with constant coefficients. These algorithms rely on homotopy continuation, monodromy loops, and Cauchy integrals. Several examples are presented which are computed using an implementation that builds on the numerical algebraic geometry software Bertini.

Sophia J. Abraham, Kehelwala D. G. Maduranga, Jeffery Kinnison et al.

Machine learning has achieved remarkable success over the past couple of decades, often attributed to a combination of algorithmic innovations and the availability of high-quality data available at scale. However, a third critical component is the fine-tuning of hyperparameters, which plays a pivotal role in achieving optimal model performance. Despite its significance, hyperparameter optimization (HPO) remains a challenging task for several reasons. Many HPO techniques rely on naive search methods or assume that the loss function is smooth and continuous, which may not always be the case. Traditional methods, like grid search and Bayesian optimization, often struggle to quickly adapt and efficiently search the loss landscape. Grid search is computationally expensive, while Bayesian optimization can be slow to prime. Since the search space for HPO is frequently high-dimensional and non-convex, it is often challenging to efficiently find a global minimum. Moreover, optimal hyperparameters can be sensitive to the specific dataset or task, further complicating the search process. To address these issues, we propose a new hyperparameter optimization method, HomOpt, using a data-driven approach based on a generalized additive model (GAM) surrogate combined with homotopy optimization. This strategy augments established optimization methodologies to boost the performance and effectiveness of any given method with faster convergence to the optimum on continuous, discrete, and categorical domain spaces. We compare the effectiveness of HomOpt applied to multiple optimization techniques (e.g., Random Search, TPE, Bayes, and SMAC) showing improved objective performance on many standardized machine learning benchmarks and challenging open-set recognition tasks.

Paul Breiding, John Cobb, Aviva K. Englander et al.

Many hypersurfaces in algebraic geometry, such as discriminants, arise as the projection of another variety. The real complement of such a hypersurface partitions its ambient space into open regions. In this paper, we propose a new method for computing these regions. Existing methods for computing regions require the explicit equation of the hypersurface as input. However, computing this equation by elimination can be computationally demanding or even infeasible. Our approach instead derives from univariate interpolation by computing the intersection of the hypersurface with a line. Such an intersection can be done using so-called pseudo-witness sets without computing a defining equation for the hypersurface - we perform elimination without actually eliminating. We implement our approach in a forthcoming Julia package and demonstrate, on several examples, that the resulting algorithm accurately recovers all regions of the real complement of a hypersurface.

Sophia J. Abraham, Jin Huang, Brandon RichardWebster et al.

We introduce a unique semantic segmentation dataset of 6,096 high-resolution aerial images capturing indigenous and invasive grass species in Bega Valley, New South Wales, Australia, designed to address the underrepresented domain of ecological data in the computer vision community. This dataset presents a challenging task due to the overlap and distribution of grass species, which is critical for advancing models in ecological and agronomical applications. Our study features a homotopy-based multi-objective fine-tuning approach that balances segmentation accuracy and contextual consistency, applicable to various models. By integrating DiceCELoss for pixel-wise classification and a smoothness loss for spatial coherence, this method evolves during training to enhance robustness against noisy data. Performance baselines are established through a case study on the Segment Anything Model (SAM), demonstrating its effectiveness. Our annotation methodology, emphasizing pen size, zoom control, and memory management, ensures high-quality dataset creation. The dataset and code will be made publicly available, aiming to drive research in computer vision, machine learning, and ecological studies, advancing environmental monitoring and sustainable development.

Sophia J. Abraham, Jonathan D. Hauenstein, Walter J. Scheirer

Face parsing is a fundamental task in computer vision, enabling applications such as identity verification, facial editing, and controllable image synthesis. However, existing face parsing models often lack fairness and robustness, leading to biased segmentation across demographic groups and errors under occlusions, noise, and domain shifts. These limitations affect downstream face synthesis, where segmentation biases can degrade generative model outputs. We propose a multi-objective learning framework that optimizes accuracy, fairness, and robustness in face parsing. Our approach introduces a homotopy-based loss function that dynamically adjusts the importance of these objectives during training. To evaluate its impact, we compare multi-objective and single-objective U-Net models in a GAN-based face synthesis pipeline (Pix2PixHD). Our results show that fairness-aware and robust segmentation improves photorealism and consistency in face generation. Additionally, we conduct preliminary experiments using ControlNet, a structured conditioning model for diffusion-based synthesis, to explore how segmentation quality influences guided image generation. Our findings demonstrate that multi-objective face parsing improves demographic consistency and robustness, leading to higher-quality GAN-based synthesis.

Emma R. Cobian, Jonathan D. Hauenstein, Fang Liu et al.

Approximating probability distributions can be a challenging task, particularly when they are supported over regions of high geometrical complexity or exhibit multiple modes. Annealing can be used to facilitate this task which is often combined with constant a priori selected increments in inverse temperature. However, using constant increments limit the computational efficiency due to the inability to adapt to situations where smooth changes in the annealed density could be handled equally well with larger increments. We introduce AdaAnn, an adaptive annealing scheduler that automatically adjusts the temperature increments based on the expected change in the Kullback-Leibler divergence between two distributions with a sufficiently close annealing temperature. AdaAnn is easy to implement and can be integrated into existing sampling approaches such as normalizing flows for variational inference and Markov chain Monte Carlo. We demonstrate the computational efficiency of the AdaAnn scheduler for variational inference with normalizing flows on a number of examples, including density approximation and parameter estimation for dynamical systems.

Edgar A. Bernal, Jonathan D. Hauenstein, Dhagash Mehta et al.

Parameterized systems of polynomial equations arise in many applications in science and engineering with the real solutions describing, for example, equilibria of a dynamical system, linkages satisfying design constraints, and scene reconstruction in computer vision. Since different parameter values can have a different number of real solutions, the parameter space is decomposed into regions whose boundary forms the real discriminant locus. This article views locating the real discriminant locus as a supervised classification problem in machine learning where the goal is to determine classification boundaries over the parameter space, with the classes being the number of real solutions. For multidimensional parameter spaces, this article presents a novel sampling method which carefully samples the parameter space. At each sample point, homotopy continuation is used to obtain the number of real solutions to the corresponding polynomial system. Machine learning techniques including nearest neighbor and deep learning are used to efficiently approximate the real discriminant locus. One application of having learned the real discriminant locus is to develop a real homotopy method that only tracks the real solution paths unlike traditional methods which track all~complex~solution~paths. Examples show that the proposed approach can efficiently approximate complicated solution boundaries such as those arising from the equilibria of the Kuramoto model.

Jonathan D. Hauenstein, Samantha N. Sherman

Synthesis problems for linkages in kinematics often yield large structured parameterized polynomial systems which generically have far fewer solutions than traditional upper bounds would suggest. This paper describes statistical models for estimating the generic number of solutions of such parameterized polynomial systems. The new approach extends previous work on success ratios of parameter homotopies to using monodromy loops as well as the addition of a trace test that provides a stopping criterion for validating that all solutions have been found. Several examples are presented demonstrating the method including Watt I six-bar motion generation problems.

Dhagash Mehta, Tianran Chen, Tingting Tang et al.

By using the viewpoint of modern computational algebraic geometry, we explore properties of the optimization landscapes of the deep linear neural network models. After clarifying on the various definitions of "flat" minima, we show that the geometrically flat minima, which are merely artifacts of residual continuous symmetries of the deep linear networks, can be straightforwardly removed by a generalized $L_2$ regularization. Then, we establish upper bounds on the number of isolated stationary points of these networks with the help of algebraic geometry. Using these upper bounds and utilizing a numerical algebraic geometry method, we find all stationary points of modest depth and matrix size. We show that in the presence of the non-zero regularization, deep linear networks indeed possess local minima which are not the global minima. Our computational results clarify certain aspects of the loss surfaces of deep linear networks and provide novel insights.

Daniel A. Brake, Jonathan D. Hauenstein, Alan C. Liddell

Computing the real solutions to a system of polynomial equations is a challenging problem, particularly verifying that all solutions have been computed. We describe an approach that combines numerical algebraic geometry and sums of squares programming to test whether a given set is "complete" with respect to the real solution set. Specifically, we test whether the Zariski closure of that given set is indeed equal to the solution set of the real radical of the ideal generated by the given polynomials. Examples with finitely and infinitely many real solutions are provided, along with an example having polynomial inequalities.

Daniel J. Bates, Jonathan D. Hauenstein, Matthew E. Niemerg et al.

We give a Descartes'-like bound on the number of positive solutions to a system of fewnomials that holds when its exponent vectors are not in convex position and a sign condition is satisfied. This was discovered while developing algorithms and software for computing the Gale transform of a fewnomial system, which is our main goal. This software is a component of a package we are developing for Khovanskii-Rolle continuation, which is a numerical algorithm to compute the real solutions to a system of fewnomials.

Jonathan D. Hauenstein, Nickolas Hein, Frank Sottile

Formulating a Schubert problem as the solutions to a system of equations in either Plücker space or in the local coordinates of a Schubert cell typically involves more equations than variables. We present a novel primal-dual formulation of any Schubert problem on a Grassmannian or flag manifold as a system of bilinear equations with the same number of equations as variables. This formulation enables numerical computations in the Schubert calculus to be certified using algorithms based on Smale's α-theory.