We define a polyhedral version of a stratification for projection maps that applies to any complex or real toric variety and show that it yields similarly desirable properties to the classical map stratification of a proper map. Our results are constructive and give rise to a method for associating the Whitney strata of the projection to the faces of the polytope of the corresponding toric variety. For all the examples we consider, our resulting algorithm outperforms known general purpose methods, e.g., Helmer and Nanda (FoCM, 2022), and Đinh and Jelonek (DCG, 2021), for computing map stratifications.

We study bihomogeneous systems defining, non-zero dimensional, biprojective varieties for which the projection onto the first group of variables results in a finite set of points. To compute (with) the 0-dimensional projection and the corresponding quotient ring, we introduce linear maps that greatly extend the classical multiplication maps for zero-dimensional systems, but are not those associated to the elimination ideal; we also call them multiplication maps. We construct them using linear algebra on the restriction of the ideal to a carefully chosen bidegree or, if available, from an arbitrary Gröbner basis. The multiplication maps allow us to compute the elimination ideal of the projection, by generalizing FGLM algorithm to bihomogenous, non-zero dimensional, varieties. We also study their properties, like their minimal polynomials and the multiplicities of their eigenvalues, and show that we can use the eigenvalues to compute numerical approximations of the zero-dimensional projection. Finally, we establish a single exponential complexity bound for computing multiplication maps and Gröbner bases, that we express in terms of the bidegrees of the generators of the corresponding bihomogeneous ideal.

We propose novel randomized geometric tools to detect low-volatility anomalies in stock markets; a principal problem in financial economics. Our modeling of the (detection) problem results in sampling and estimating the (relative) volume of geodesically non-convex and non-connected spherical patches that arise by intersecting a non-standard simplex with a sphere. To sample, we introduce two novel Markov Chain Monte Carlo (MCMC) algorithms that exploit the geometry of the problem and employ state-of-the-art continuous geometric random walks (such as Billiard walk and Hit-and-Run) adapted on spherical patches. To our knowledge, this is the first geometric formulation and MCMC-based analysis of the volatility puzzle in stock markets. We have implemented our algorithms in C++ (along with an R interface) and we illustrate the power of our approach by performing extensive experiments on real data. Our analyses provide accurate detection and new insights into the distribution of portfolios’ performance characteristics. Moreover, we use our tools to show that classical methods for low-volatility anomaly detection in finance form bad proxies that could lead to misleading or inaccurate results.

We present a method for solving two minimal problems for relative camera pose estimation from three views, which are based on three view correspondences of (i) three points and one line and the novel case of (ii) three points and two lines through two of the points. These problems are too difficult to be efficiently solved by the state of the art Gröbner basis methods. Our method is based on a new efficient homotopy continuation (HC) solver framework MINUS, which dramatically speeds up previous HC solving by specializing hc methods to generic cases of our problems. We characterize their number of solutions and show with simulated experiments that our solvers are numerically robust and stable under image noise, a key contribution given the borderline intractable degree of nonlinearity of trinocular constraints. We show in real experiments that (i) sift feature location and orientation provide good enough point-and-line correspondences for three-view reconstruction and (ii) that we can solve difficult cases with too few or too noisy tentative matches, where the state of the art structure from motion initialization fails.

We introduce Reflective Hamiltonian Monte Carlo (ReHMC), an HMC-based algorithm to sample from a log-concave distribution restricted to a convex body. The random walk is based on incorporating reflections to the Hamiltonian dynamics such that the support of the target density is the convex body. We develop an efficient open source implementation of ReHMC and perform an experimental study on various high-dimensional datasets. The experiments suggest that ReHMC outperforms Hit-and-Run and Coordinate-Hit-and-Run regarding the time it needs to produce an independent sample, introducing practical truncated sampling in thousands of dimensions.

Systems of polynomial equations arise frequently in computer vision, especially in multiview geometry problems. Traditional methods for solving these systems typically aim to eliminate variables to reach a univariate polynomial, e.g., a tenth-order polynomial for 5-point pose estimation, using clever manipulations, or more generally using Grobner basis, resultants, and elimination templates, leading to successful algorithms for multiview geometry and other problems. However, these methods do not work when the problem is complex and when they do, they face efficiency and stability issues. Homotopy Continuation (HC) can solve more complex problems without the stability issues, and with guarantees of a global solution, but they are known to be slow. In this paper we show that HC can be parallelized on a GPU, showing significant speedups up to 56 times on polynomial benchmarks. We also show that GPU-HC can be generically applied to a range of computer vision problems, including 4-view triangulation and trifocal pose estimation with unknown focal length, which cannot be solved with elimination template but they can be efficiently solved with HC. GPU-HC opens the door to easy formulation and solution of a range of computer vision problems.

Isolating the real roots of univariate polynomials is a fundamental problem in symbolic computation and it is arguably one of the most important problems in computational mathematics. The problem has a long history decorated with numerous ingenious algorithms and furnishes an active area of research. However, the worst-case analysis of root-finding algorithms does not correlate with their practical performance. We develop a smoothed analysis framework for polynomials with integer coefficients to bridge the gap between the complexity estimates and the practical performance. In this setting, we derive that the expected bit complexity of Descartes solver to isolate the real roots of a polynomial, with coefficients uniformly distributed, is ÕB(d2 + dτ), where d is the degree of the polynomial and τ the bitsize of the coefficients.

Motivated by applications in combinatorial geometry, we consider the following question: Let $\lambda=(\lambda_1,\lambda_2,\ldots,\lambda_m)$ be an $m$-partition of a positive integer $n$, $S_i \subseteq \mathbb{C}^{\lambda_i}$ be finite sets, and let $S:=S_1 \times S_2 \times \cdots \times S_m \subset \mathbb{C}^n$ be the multigrid defined by $S_i$. Suppose $p$ is an $n$-variate degree $d$ polynomial. How many zeros does $p$ have on $S$? We first develop a multivariate generalization of the combinatorial nullstellensatz that certifies existence of a point $t \in S$ so that $p(t) \neq 0$. Then we show that a natural multivariate generalization of the DeMillo--Lipton--Schwartz--Zippel lemma holds, except for a special family of polynomials that we call $\lambda$-reducible. This yields a simultaneous generalization of the Szemerédi--Trotter theorem and the Schwartz--Zippel lemma into higher dimensions, and has applications in incidence geometry. Finally, we develop a symbolic algorithm that identifies certain $\lambda$-reducible polynomials. More precisely, our symbolic algorithm detects polynomials that include a Cartesian product of hypersurfaces in their zero set. It is likely that using Chow forms the algorithm can be generalized to handle arbitrary $\lambda$-reducible polynomials, which we leave as an open problem.

status: Accepted

We present a method for solving two minimal problems for relative camera pose estimation from three views, which are based on three view correspondences of (i) three points and one line and (ii) three points and two lines through two of the points. These problems are too difficult to be efficiently solved by the state of the art Grobner basis methods. Our method is based on a new efficient homotopy continuation (HC) solver, which dramatically speeds up previous HC solving by specializing HC methods to generic cases of our problems. We show in simulated experiments that our solvers are numerically robust and stable under image noise. We show in real experiment that (i) SIFT features provide good enough point-and-line correspondences for three-view reconstruction and (ii) that we can solve difficult cases with too few or too noisy tentative matches where the state of the art structure from motion initialization fails.

We consider the problem of computing the topology and describing the geometry of a parametric curve in Rn. We present an algorithm, PTOPO, that constructs an abstract graph that is isotopic to the curve in the embedding space. Our method exploits the benefits of the parametric representation and does not resort to implicitization.

The condition-based complexity analysis framework is one of the gems of modern numerical algebraic geometry and theoretical computer science. Among the challenges that it poses is to expand the currently limited range of random polynomials that we can handle. Despite important recent progress, the available tools cannot handle random sparse polynomials and Gaussian polynomials, that is polynomials whose coefficients are i.i.d. Gaussian random variables. We initiate a condition-based complexity framework based on the norm of the cube that is a step in this direction. We present this framework for real hypersurfaces and univariate polynomials. We demonstrate its capabilities in two problems, under very mild probabilistic assumptions. On the one hand, we show that the average run-time of the Plantinga-Vegter algorithm is polynomial in the degree for random sparse (alas a restricted sparseness structure) polynomials and random Gaussian polynomials. On the other hand, we study the size of the subdivision tree for Descartes' solver and run-time of the solver by Jindal and Sagraloff (2017). In both cases, we provide a bound that is polynomial in the size of the input (size of the support plus the logarithm of the degree) not only for the average but also for all higher moments.

Grö bner bases is one the most powerful tools in algorithmic nonlinear algebra. Their computation is an intrinsically hard problem with a complexity at least single exponential in the number of variables. However, in most of the cases, the polynomial systems coming from applications have some kind of structure. We consider sparse systems where the input polynomials have a few non-zero terms. Our approach to exploit sparsity is to embed the systems in a semigroup algebra and to compute Grö bner bases over this algebra. Up to now, the algorithms that follow this approach benefit from the sparsity only in the case where all the polynomials have the same sparsity structure, that is the same Newton polytope. We introduce the first algorithm that overcomes this restriction. Under regularity assumptions, it performs no redundant computations. Further, we extend this algorithm to compute Grö bner basis in the standard algebra and solve sparse polynomials systems over the torus (\mathbbC ^*)^n. The complexity of the algorithm depends on the Newton polytopes.

The number of embeddings of minimally rigid graphs in RD is (by definition) finite, modulo rigid transformations, for every generic choice of edge lengths. Even though various approaches have been proposed to compute it, the gap between upper and lower bounds is still enormous. Specific values and its asymptotic behavior are major and fascinating open problems in rigidity theory. Our work considers the maximal number of real embeddings of minimally rigid graphs in R3. We modify a commonly used parametric semi-algebraic formulation that exploits the Cayley-Menger determinant to minimize the a priori number of complex embeddings, where the parameters correspond to edge lengths. To cope with the huge dimension of the parameter space and find specializations of the parameters that maximize the number of real embeddings, we introduce a method based on coupler curves that makes the sampling feasible for spatial minimally rigid graphs. Our methodology results in the first full classification of the number of real embeddings of graphs with 7 vertices in R3, which was the smallest open case. Building on this and certain 8-vertex graphs, we improve the previously known general lower bound on the maximum number of real embeddings in R3.

A fundamental problem in computational algebraic geometry is the computation of the resultant. A central question is when and how to compute it as the determinant of a matrix whose elements are the coefficients of the input polynomials up-to sign. This problem is well understood for unmixed multihomogeneous systems, that is for systems consisting of multihomogeneous polynomials with the same support. However, little is known for mixed systems, that is for systems consisting of polynomials with different supports. We consider the computation of the multihomogeneous resultant of bilinear systems involving two different supports. We present a constructive approach that expresses the resultant as the exact determinant of a Koszul resultant matrix, that is a matrix constructed from maps in the Koszul complex. % We exploit the resultant matrix to propose an algorithm to solve such systems. In the process we extend the classical eigenvalues and eigenvectors criterion to a more general setting. Our extension of the eigenvalues criterion applies to a general class of matrices, including the Sylvester-type and the Koszul-type ones.

One of the biggest open problems in computational algebra is the design of efficient algorithms for Gröbner basis computations that take into account the sparsity of the input polynomials. We can perform such computations in the case of unmixed polynomial systems, that is systems with polynomials having the same support, using the approach of Faugère, Spaenlehauer, and Svartz [ISSAC'14]. We present two algorithms for sparse Gröbner bases computations for mixed systems. The first one computes with mixed sparse systems and exploits the supports of the polynomials. Under regularity assumptions, it performs no reductions to zero. For mixed, square, and 0-dimensional multihomogeneous polynomial systems, we present a dedicated, and potentially more efficient, algorithm that exploits different algebraic properties that performs no reduction to zero. We give an explicit bound for the maximal degree appearing in the computations.

We present the first complexity analysis of the algorithm by Plantinga and Vegter for approximating real implicit curves and surfaces. This approximation algorithm certifies the topological correctness of the output using both subdivision and interval arithmetic. In practice, it has been seen to be quite efficient; our goal is to quantify this efficiency. We focus on the subdivision step (and not the approximation step) of the Plantinga and Vegter algorithm. We begin by extending the subdivision step to arbitrary dimensions. We provide a priori worst-case bounds on the complexity of this algorithm both in terms of the number of subregions constructed and the bit complexity for the construction. Then, we use continuous amortization to derive adaptive bounds on the complexity of the subdivided region. We also provide examples showing our bounds are tight.

We present explicit worst case degree and height bounds for the rational univariate representation of the isolated roots of polynomial systems based on mixed volume. We base our estimations on height bounds of resultants and we consider the case of 0-dimensional, positive dimensional, and parametric polynomial systems.

Optimal resultant formulas have been systematically constructed mostly for unmixed polynomial systems, that is, systems of polynomials which all have the same support. However, such a condition is restrictive, since mixed systems of equations arise frequently in practical problems.

Symmetric Tensor Decomposition is a major problem that arises in areas such as signal processing, statistics, data analysis and computational neuroscience. It is equivalent to write a homogeneous polynomial in $n$ variables of degree $D$ as a sum of $D$-th powers of linear forms, using the minimal number of summands. This minimal number is called the rank of the polynomial/tensor. We consider the decomposition of binary forms, that corresponds to the decomposition of symmetric tensors of dimension $2$ and order $D$. This problem has its roots in Invariant Theory, where the decompositions are known as canonical forms. As part of that theory, different algorithms were proposed for the binary forms. In recent years, those algorithms were extended for the general symmetric tensor decomposition problem. We present a new randomized algorithm that enhances the previous approaches with results from structured linear algebra and techniques from linear recurrent sequences. It achieves a softly linear arithmetic complexity bound. To the best of our knowledge, the previously known algorithms have quadratic complexity bounds. We compute a symbolic minimal decomposition in O(M(D) log(D)) arithmetic operations, where M(D) is the complexity of multiplying two polynomials of degree D. We approximate the terms of the decomposition with an error of 2-ε, in O(D log2(D) (log2(D) + log(ε))) arithmetic operations. To bound the size of the representation of the coefficients involved in the decomposition, we bound the algebraic degree of the problem by min(rank, D-rank+1). When the input polynomial has integer coefficients, our algorithm performs, up to poly-logarithmic factors, OB(D l + D4 + D3 τ) bit operations, where τ is the maximum bitsize of the coefficients and 2-l is the relative error of the terms in the decomposition.

We bound the Boolean complexity of computing isolating hyperboxes for all complex roots of systems of bilinear polynomials. The resultant of such systems admits a family of determinantal Sylvester-type formulas, which we make explicit by means of homological complexes. The computation of the determinant of the resultant matrix is a bottleneck for the overall complexity. We exploit the quasi-Toeplitz structure to reduce the problem to efficient matrix-vector products, corresponding to multivariate polynomial multiplication. For zero-dimensional systems, we arrive at a primitive element and a rational univariate representation of the roots. The overall bit complexity of our probabilistic algorithm is OB(n4 D4 + n2D4 τ), where n is the number of variables, D equals the bilinear Bezout bound, and τ is the maximum coefficient bitsize. Finally, a careful infinitesimal symbolic perturbation of the system allows us to treat degenerate and positive dimensional systems, thus making our algorithms and complexity analysis applicable to the general case.

We assume that a real square-free polynomial A has a degree d, a maximum coefficient bitsize τ and a real root lying in an isolating interval and having no nonreal roots nearby (we quantify this assumption). Then, we combine the Double Exponential Sieve algorithm (also called the Bisection of the Exponents), the bisection, and Newton iteration to decrease the width of this inclusion interval by a factor of t=2-L. The algorithm has Boolean complexity ÕB(d2 τ + d L ). Our algorithms support the same complexity bound for the refinement of r roots, for any r ≤ d.

We present algorithmic, complexity and implementation results for the problem of isolating the real roots of a univariate polynomial in Bα ∈ L[y], where L = Q(α1,..., αℓ) is an algebraic extension of the rational numbers. Our bounds are single exponential in ℓ and match the ones presented in [34] for the case ℓ = 1. We consider two approaches. The first, indirect approach, using multivariate resultants, computes a univariate polynomial with integer coefficients, among the real roots of which are the real roots of Bα. The Boolean complexity of this approach is OB(N4ℓ+4), where N is the maximum of the degrees and the coefficient bitsize of the involved polynomials. The second, direct approach, tries to solve the polynomial directly, without reducing the problem to a univariate one. We present an algorithm that generalizes Sturm algorithm from the univariate case, and modified versions of well known solvers that are either numerical or based on Descartes' rule of sign. We achieve a Boolean complexity of OB [equation], respectively. We implemented the algorithms in C as part of the core library of Mathematica and we illustrate their efficiency over various data sets.

We present algorithmic, complexity and implementation results for the problem of isolating the real roots of a univariate polynomial in Bα ∈ L[y], where L=Qα is a simple algebraic extension of the rational numbers. We revisit two approaches for the problem. In the first approach, using resultant computations, we perform a reduction to a polynomial with integer coefficients and we deduce a bound of OB(N10) for isolating the real roots of Bα, where N is an upper bound on all the quantities (degree and bitsize) of the input polynomials. In the second approach we isolate the real roots working directly on the polynomial of the input. We compute improved separation bounds for the roots and we prove that they are optimal, under mild assumptions. For isolating the real roots we consider a modified Sturm algorithm, and a modified version of Descartes' algorithm introduced by Sagraloff. For the former we prove a complexity bound of OB(N8) and for the latter a bound of OB(N7). We implemented the algorithms in C as part of the core library of Mathematica and we illustrate their efficiency over various data sets. Finally, we present complexity results for the general case of the first approach, where the coefficients belong to multiple extensions.

Shapley's discounted stochastic games, Everett's recursive games and Gillette's undiscounted stochastic games are classical models of game theory describing two-player zero-sum games of potentially infinite duration. We describe algorithms for exactly solving these games. When the number of positions of the game isbconstant, our algorithms run in polynomial time.

Tensor decompositions permit to estimate in a deterministic way the parameters in a multi-linear model. Applications have been already pointed out in antenna array processing and digital communications, among others, and are extremely attractive provided some diversity at the receiver is available. As opposed to the widely used ALS algorithm, non-iterative algorithms are proposed in this paper to compute the required tensor decomposition into a sum of rank-1 terms, when some factor matrices enjoy some structure, such as block-Hankel, triangular, band, etc.

Our probabilistic analysis sheds light to the following questions: Why do random polynomials seem to have few, and well separated real roots, on the average? Why do exact algorithms for real root isolation may perform comparatively well or even better than numerical ones?

In this paper we derive aggregate separation bounds, named after Davenport-Mahler-Mignotte (DMM), on the isolated roots of polynomial systems, specifically on the minimum distance between any two such roots. The bounds exploit the structure of the system and the height of the sparse (or toric) resultant by means of mixed volume, as well as recent advances on aggregate root bounds for univariate polynomials, and are applicable to arbitrary positive dimensional systems. We improve upon Canny's gap theorem [7] by a factor of O(dn-1), where d bounds the degree of the polynomials, and n is the number of variables. One application is to the bitsize of the eigenvalues and eigenvectors of an integer matrix, which also yields a new proof that the problem is polynomial. We also compare against recent lower bounds on the absolute value of the root coordinates by Brownawell and Yap [5], obtained under the hypothesis there is a 0-dimensional projection. Our bounds are in general comparable, but exploit sparseness; they are also tighter when bounding the value of a positive polynomial over the simplex. For this problem, we also improve upon the bounds in [2, 16]. Our analysis provides a precise asymptotic upper bound on the number of steps that subdivision-based algorithms perform in order to isolate all real roots of a polynomial system. This leads to the first complexity bound of Milne's algorithm [22] in 2D.

We revisit the problem of computing the topology and geometry of a real algebraic plane curve. The topology is of prime interest but geometric information, such as the position of singular and critical points, is also relevant. A challenge is to compute efficiently this information for the given coordinate system even if the curve is not in generic position.

We examine the problem of computing exactly the Delaunay graph (and the dual Voronoi diagram) of a set of, possibly intersecting, smooth convex pseudo-circles in the Euclidean plane, given in parametric form. Pseudo-circles are (convex) sites, every pair of which has at most two intersecting points. The Delaunay graph is constructed incrementally. Our first contribution is to propose robust end efficient algorithms for all required predicates, thus generalizing our earlier algorithms for ellipses, and we analyze their algebraic complexity, under the exact computation paradigm. Second, we focus on InCircle, which is the hardest predicate, and express it by a simple sparse 5 X 5 polynomial system, which allows for an efficient implementation by means of successive Sylvester resultants and a new factorization lemma. The third contribution is our cgal-based c++ software for the case of ellipses, which is the first exact implementation for the problem. Our code spends about 98 sec to construct the Delaunay graph of 128 non-intersecting ellipses, when few degeneracies occur. It is faster than the cgal segment Delaunay graph, when ellipses are approximated by k-gons for k > 15.

Based on precomputed Sturm–Habicht sequences, discriminants and invariants, we classify, isolate with rational points, and compare the real roots of polynomials of degree up to 4. In particular, we express all isolating points as rational functions of the input polynomial coefficients. Although the roots are algebraic numbers and can be expressed by radicals, such representation involves some roots of complex numbers. This is inefficient, and hard to handle in applications in geometric computing and quantifier elimination. We also define rational isolating points between the roots of the quintic. We combine these results with a simple version of Rational Univariate Representation to isolate all common real roots of a bivariate system of rational polynomials of total degree ≤2 and to compute the multiplicity of these roots. We present our software within library synaps and perform experiments and comparisons with several public-domain implementations. Our package is 2–10 times faster than numerical methods and exact subdivision-based methods, including software with intrinsic filtering.

We consider exact real solving of well-constrained, bivariate systems of relatively prime polynomials. The main problem is to compute all common real roots in isolating interval representation, and to determine their intersection multiplicities. We present three algorithms and analyze their asymptotic bit complexity, obtaining a bound of ÕB(N14) for the purely projection-based method, and ÕB(N12) for two subresultants-based methods: these ignore polylogarithmic factors, and N bounds the degree and the bitsize of the polynomials. The previous record bound was ÕB(N14).

This paper examines the computation of the Voronoi diagram of a set of ellipses in the Euclidean plane. We propose the first complete algorithms, under the exact computation paradigm, for the predicates of an incremental algorithm: κ1 decides which one of 2 given ellipses is closest to a given exterior point; κ2 decides the position of a query ellipse relative to an external bitangent line of 2 given ellipses; κ3 decides the position of a query ellipse relative to a Voronoi circle of 3 given ellipses; κ4 determines the type of conflict between a Voronoi edge, defined by 4 given ellipses, and a query ellipse. The paper is restricted to non-intersecting ellipses, but the extension to arbitrary ones is straightforward. The ellipses are input in parametric representation or constructively. For κ1 and κ2 we derive optimal algebraic conditions, solve them exactly and provide efficient implementations in C++. For κ3 we compute a tight bound on the number of complex tritangent circles and use the parametric form of the ellipses in order to design an exact subdivision-based algorithm, which is implemented on Maple. This approach essentially answers κ4 as well. We conclude with current work on optimizing κ3 and implementing it in C++.

Our work goes towards answering the growing need for the robust and efficient manipulation of curved objects in numerous applications. The kernel of the CGAL library provides several functionalities which are, however, mostly restricted to linear objects. We focus here on the arrangement of conic arcs in the plane. Our first contribution is the design, implementation and testing of a kernel for computing arrangements of circular arcs.A preliminary C++ implementation exists also for arbitrary conic curves. We discuss the representation and predicates of the geometric objects. Our implementation is targeted for inclusion in the CGAL library. Our second contribution concerns exact and efficient algebraic algorithms for the case of conics. They treat all inputs, including degeneracies, and they are implemented as part of the library SYNAPS 2.1.Our tools include Sturm sequences, resultants, Descartes' rule, andisolating points. Thirdly, our experiments on circular arcs show that our methods compare favorably to existing alternatives using CORE 1.6x and LEDA 4.5.