Computing Arrangements of Hypersurfaces

Arrangements of real hyperplanes are ubiquitous in combinatorics, algebra and geometry. Their regions (i.e. connected components) are convex polyhedra, either bounded or unbounded, and their numbers are invariants of the underlying matroid [16]. The number of bounded regions agrees with the Euler characteristic of the complex arrangement complement [6, 7], and also with the maximum likelihood degree (ML degree). The latter is the number of complex critical points of the associated models in statistics and physics [17].

Arrangements of nonlinear hypersurfaces are equally important, but they are studied much less. Polynomials of higher degree create features that are not seen when dealing with hyperplanes. The number of regions is no longer a combinatorial invariant, but it depends in a subtle way on the coefficients. Moreover, the regions are generally not contractible.

We present a practical software tool, called HypersurfaceRegions.jl and implemented in the programming language Julia [1], whose input consists of $k$ polynomials in $n$ variables:

The output is a list of all regions $C$ of the $n$-dimensional manifold

The list is grouped according to sign vectors $\sigma\in\{-1,+1\}^{k}$, where $\sigma_{i}$ is the sign of $f_{i}$ on $C$. Unlike in the case of hyperplane arrangements, each sign vector $\sigma$ typically corresponds to multiple regions. For each region $C$ we find the Euler characteristic via a Morse function.

Our software HypersurfaceRegions.jl also features heuristics for deciding whether a region is bounded or unbounded, and which of the regions are fused when the hyperplane at infinity is added. In our Section 4, titled How to Use the Software, we explain how this works.

The method we present is an adaptation of the algorithm by Cummings, Hauenstein, Hong and Smyth [8], which in turn is inspired by [11]. Their setting is more general in that they allow semialgebraic sets defined by both equations and inequalities. We restrict ourselves to arrangement complements, defined by inequations $f_{1}\not=0,\ldots,f_{k}\not=0$. This enables us to offer tools in Julia [1] that are easy to use and widely applicable. Our implementation is based on the numerical algebraic geometry software HomotopyContinuation.jl [5] and on the software DifferentialEquations.jl [14] for solving differential equations.

This paper is organized as follows. In Section 2 we introduce a Morse function for our problem. This is a modified version of the log-likelihood function associated with $f_{1},\ldots,f_{k}$. Every region of $\mathcal{U}$ contains at least one critical point. We compute the critical points using HomotopyContinuation.jl . These points determine the Euler characteristic of each region.

In Section 3 we turn to the Mountain Pass Theorem [2, 12], which ensures that we obtain a connected graph in each region when tracking paths from index $1$ critical points to index $0$ critical points. This tracking procedure is realized in our software by solving an ODE using DifferentialEquations.jl. In short, our software is built on flows in Morse theory.

Section 4 explains how to run HypersurfaceRegions.jl. It is aimed at beginners with no prior experience with numerical software. We demonstrate that HypersurfaceRegions.jl can be used for a wide range of scenarios from the mathematical spectrum. In Section 5 we report on test runs of our software on generic instances, obtained by sampling random hypersurfaces and random spectrahedra. Our presentation emphasizes ease and simplicity.

The algorithm in [8] rests on Morse theory. We review some basics from the textbook [13]. Let $\mathcal{M}$ be an $n$-dimensional manifold. A smooth function $g:\mathcal{M}\rightarrow\mathbb{R}$ is a Morse function if the Hessian matrix $\bigl{(}\partial^{2}g/\partial x_{i}\partial x_{j}\bigr{)}$ is invertible at every critical point of $g$. The number of positive eigenvalues of the Hessian matrix is the index of the critical point. Let $g$ be a Morse function which is exhaustive, i.e. the superlevel set $\{g\geq c\}=\{u\in\mathcal{M}:g(u)\geq c\}$ is compact for each $c\in\mathbb{R}$. If an interval $[a,b]\subset\mathbb{R}$ contains no critical value of $g$ then the superlevel sets $\{g\geq c\}$ are diffeomorphic for $c\in[a,b]$. By contrast, suppose $c$ is a critical value of $g$, corresponding to a unique critical point of index $\ell$. Then, for sufficiently small $\epsilon$, the superlevel set $\{g\geq c-\epsilon\}$ is diffeomorphic to $\{g\geq c+\epsilon\}$ with one $\ell$-handle of dimension $n$ attached. Since $\ell$-handles can be shrunk to $\ell$-cells, one arrives at the following result.

We shall apply Theorem 3 to the manifold $\mathcal{U}$ defined in (2). This requires an exhaustive Morse function. On each region of $\,\mathcal{U}$, this role will be played by the rational function

For the denominator $q(x)$ we take a generic polynomial of degree $2$ that is strictly positive on $\mathbb{R}^{n}$. The exponents $s_{1},\ldots,s_{k},t$ are arbitrary positive integers which satisfy the inequality

The Morse function $g(x)$ in (4) is called master function in the theory of arrangements [7] and likelihood function in algebraic statistics [17]. Its logarithm is the log-likelihood function

The critical points of $g(x)$ are the solutions in $\mathbb{C}^{n}$ of the equation $\nabla\,{\rm log}(g(x))=0$. The generic choice of the quadric $q(x)$ implies that all critical points have distinct critical values. The number of critical values is the Euler characteristic of the very affine complex variety

The construction of $g(x)$ ensures that each region of $\mathcal{U}$ contains at least one critical point. Our algorithm computes all real critical points, and it connects them via the Mountain Pass Theorem. This will be explained in Section 3. We now first return to our three ellipsoids.

The signed Euler characteristic of $\mathcal{U}$ is known as the maximum likelihood degree (ML degree) in algebraic statistics [6, 17]. This is the number of critical points we must compute. Explicitly, the logarithmic derivative of $g(x)$ translates into the rational function equations

The following a priori bound on the number of solutions was established in [6, Theorem 1].

Since every region of our manifold $\mathcal{U}$ contains at least one critical point, we conclude:

The bound in Corollary 8 is tight for linear polynomials. This fact, and many more results on ML degrees, will be proved in a forthcoming article by Berhard Reinke and Kexin Wang. We close this section by sketching a proof for hyperplanes in general position.

Suppose $k\geq n$ and let $d_{1}=\cdots=d_{k}=1$. Then the generating function (8) becomes

We find that the coefficient of $z^{n}$ in the Taylor expansion of this rational function equals

This agrees with the number of all regions in a general arrangement of $k$ hyperplanes in $\mathbb{R}^{n}$. Thus, all complex critical points are real, and each critical point lies in its own region.

The Mountain Pass Theorem guarantees that all critical points in a region will be connected. This theorem originates in the Calculus of Variations. We learned about its importance for numerical computations in real algebraic geometry from the work of Cummings et al. [8].

We fix the Morse function $g$ on $\mathcal{U}$ that is given in (4). Since the quadric $q$ is generic, the function $g$ has only finitely many critical points, and $g$ takes distinct values at these critical points. Given any starting point $x_{0}\in\mathcal{U}$, we will reach one of these critical points by numerical path tracking. This is done by solving the ODE that describes the gradient flow

If $x_{0}$ is chosen at random then the gradient flow will reach a critical point of index $0$. Loosely speaking, hill climbing in the steepest direction will probably lead us to a mountain peak.

Let $p_{1},p_{2},\ldots,p_{d}\in\mathcal{U}$ denote the real critical points of $g$. Our aim is to build a graph with vertex set $\{p_{1},p_{2},\ldots,p_{d}\}$ whose connected components correspond to the regions of $\mathcal{U}$. For each critical point $p_{i}$, we compute the eigenvalues and eigenvectors of the Hessian $H_{p_{i}}$. This is the symmetric $n\times n$ matrix of second derivatives of $g(x)$ evaluated at $x=p_{i}$.

Note that ${\rm det}(H_{p_{i}})\not=0$ because $g$ was constructed to be a Morse function. An eigenvector $v$ of $H_{p_{i}}$ is called unstable if the corresponding eigenvalue is positive. If this holds then

If $x_{0}=p_{i}+\epsilon v$ is the starting point in (9) then the ODE will lead us to a critical point $p_{j}$ of $g$ that is different from $p_{i}$. Whenever this happens, our graph acquires the edge $p_{i}\rightarrow p_{j}$.

We now explain why this method connects all critical points in a fixed region. First of all, each $p_{i}$ is connected to some critical point of index $0$. Suppose we start at a critical point $p_{i_{1}}$ with positive index. There is a path from $p_{i_{1}}$ that limits to some critical point $p_{i_{2}}$ with $g(p_{i_{2}})>g(p_{i_{1}})$. If $p_{i_{2}}$ has positive index, then the solution path from $p_{i_{2}}$ limits to some critical point $p_{i_{3}}$ with $g(p_{i_{3}})>g(p_{i_{2}})$. Since there are only finitely many critical points, the process must terminate and the last critical point in the sequence $p_{i_{1}},p_{i_{2}},\ldots$ has index $0$.

We are left to show that all index $0$ critical points in the same region will be connected. To this end, we introduce one more tool from Morse theory. The stable manifold of $p_{i}$ is

The dimension of $M(p_{i})$ is the number of stable eigenvectors of $H_{p_{i}}$. The stable manifolds for critical points of index $0$ have full dimension $n$ in $\mathcal{U}$. So, for each region $C$ of $\mathcal{U}$, we have

where the closure is taken in $\mathcal{U}$. Consider any two critical points $p_{i_{\alpha}}$ and $p_{i_{\beta}}$ of index $0$ in $C$. Since $C$ is connected, (10) implies that there is a sequence of index $0$ critical points $p_{i_{1}},\ldots,p_{i_{s}}$ such that the intersections $\overline{M(p_{i_{\alpha}})}\cap\overline{M(p_{i_{1}})}$ and $\overline{M(p_{i_{s}})}\cap\overline{M(p_{i_{\beta}})}$ are non-empty, and also $\,\overline{M(p_{i_{j}})}\cap\overline{M(p_{i_{j+1}})}\,$ is non-empty for $j=1,\ldots,s-1$. Corollary 10 below tells us that each of these intersections contains at least one critical point of index $1$. Thus any pair of index $0$ critical points in $C$ is connected through index $1$ critical points. Therefore, the connectivity of the finite graph we are building for $C$ will be ensured by Corollary 10.

The Morse function $g:C\rightarrow\mathbb{R}$ satisfies the Palais-Smale condition. This means that every sequence $\{x_{j}\}$ in $C$ which satisfies $\,\lim_{j\rightarrow\infty}g(x_{j})=\alpha\in\mathbb{R}\,$ and $\,\lim_{j\rightarrow\infty}\nabla\,g(x_{j})=0\,$ has a convergent subsequence. The limit of such a convergent subsequence is a critical point in $C$ with critical value $\alpha$. The Palais-Smale condition holds in our situation because $g$ is positive on $C$, it tends to zero on the boundary, and $g^{-1}([\alpha_{1},\alpha_{2}])$ is compact for $0<\alpha_{1}<\alpha_{2}$.

A path between two points $a$ and $b$ in $C$ is a continuous function $\gamma:[0,1]\mapsto C$ with $\gamma(0)=a$ and $\gamma(1)=b$. We write $\Gamma_{a,b}$ for the set of all such paths. A set $S\subset C$ separates $a$ and $b$ if every path $\gamma\in\Gamma_{a,b}$ contains some point in $S$. We write $g(\gamma)=\{g(\gamma(t)):t\in[0,1]\}$. The Palais-Smale condition for $g$ is a hypothesis for the next result, which is [12, Theorem 3].

We have shown that the regions $C$ of $\mathcal{U}$ can be identified by gradient ascent from index $1$ critical points. Here we start in the two directions determined by the unique unstable eigenvector. The resulting graph on the index $0$ critical points is guaranteed to be connected. We now summarize our approach. Algorithm 1 forms the basis of HypersurfaceRegions.jl.

Frequently one is interested in the regions of an arrangement in projective space $\mathbb{R}\mathbb{P}^{n}$. Algorithm 2 addresses this point. To begin with, we take the same input as before, namely $f_{1},f_{2},\ldots,f_{k}\in\mathbb{R}[x_{1},\ldots,x_{n}]$. We further assume that Algorithm 1 has already been run.

In our implementation of Algorithm 2 we discard critical points that are close to the hypersurface $\mathbb{R}^{n-1}\backslash\mathcal{U}_{\infty}$. These points cause numerical problems when solving the corresponding ODE. Hence, that part of our software is based on a heuristic.

We now explain the distinction between bounded and undecided, alluded to in Step 7. It is a feature of numerical computations that tangencies and singularities are hard to detect.

We identify bounded regions as follows. We fix a small parameter $\delta>0$. A region in Step 7 is declared bounded if it is contained in $-1/\delta<x_{1}<1/\delta$. In practice, we compute critical points of our Morse function on the hyperplanes $x_{1}=1/\delta$ and $x_{1}=-1/\delta$. We process these critical points similar to what is done in Algorithm 2. Namely, for each critical point $q$ on the hyperplane $x_{1}=1/\delta$, we find the largest value among all zeros of the polynomials $f_{j}(t,tq)$ for $j=1,\ldots,k$ that are smaller than $1/\delta$. Fix a real number $\lambda$ larger than that value, but smaller than $1/\delta$. We then solve the ODE (9) with starting points $\lambda\cdot(1,q)$ and record the limiting critical points. We proceed similarly for the other hyperplane $x_{1}=-1/\delta$.

If a region in $\mathbb{R}^{n}$ is not visited during this process, then it does not intersect either hyperplane. If its critical points satisfy $-1/\delta<x_{1}<1/\delta$, then it is bounded. If a region is not bounded or unbounded, then it is declared undecided. Thus, undecided regions are close to being tangent to the hyperplane at infinity, where “close” refers to the tolerance $\delta$.

Our software is easy to use, even for beginners. We here offer a step-by-step introduction. For a complete overview over our implementation we refer to the online documentation at

The first step is to download the programming language Julia. One navigates to the website https://julialang.org/downloads/ for the latest version. After it is downloaded and installed, we start Julia and enter the package manager by pressing the ] button. Once in the package manager, type add HypersurfaceRegions and hit enter. This installs HypersurfaceRegions.jl. One leaves the package manager by pressing the back space button. To load our software into the current Julia session, type the following command:

Now, we are ready to use the implementation. Let’s compute our first arrangement!!

Being ambitious, we start with an example in $\mathbb{R}^{3}$. Consider the four surfaces in Figure 2. The yellow plane is defined by $z=0$. The three quadrics are invariant under rotation around the $z$-axis. The two paraboloids are defined by $z=x^{2}+y^{2}-1$ and $z=-x^{2}-y^{2}+1$. The hyperboloid is defined by $x^{2}+y^{2}-\frac{3}{2}z^{2}=\frac{1}{4}$. We input this arrangement into our software:

The output for this input is shown on the left in Figure 3. We see that $\mathcal{U}$ has $12$ regions, and each of them is uniquely identified by its sign pattern $\sigma\in\{-,+\}^{4}$. Six regions lie outside the hyperboloid, which means that $\sigma_{4}=+$. Each of them is homotopy equivalent to a circle, so $\chi=0$. The other six regions lie inside the hyperboloid, which means $\sigma_{4}=-$.

Two of the regions inside the hyperboloid are unbounded and also circular ($\chi=0$). Finally, there are two bounded regions and two regions tangent to the hyperplane at infinity. They are all contractible ($\chi=1$). We compute this information by setting a flag as follows:

The output of HypersurfaceRegions.jl for this input is shown on the right in Figure 3. This also reminds us that the number of real critical points depends on $q(x)$, which is random.

One application of our method is the study of discriminantal arrangments. Such arrangements arise whenever one is interested in the topology of real varieties that depend on parameters. The one-dimensional case of this is known as Real Root Classification, where one expresses the number of real roots of a zero-dimensional system as a function of parameters.

To use the current version of HypersurfaceRegions.jl in this context, it is assumed that the discriminant is known and that it breaks up into smaller factors. Here is an example.

You are now invited to run the following code and to match its output with Figure 1:

A range of additional features is available in HypersurfaceRegions.jl. For instance, one finds the critical points in the $i$-th region of the output by running the following command:

One can also locate the region to which a given point (e.g. the origin) belongs, as follows:

We next discuss a few implementation details. We use monodromy [9], implemented in HomotopyContinuation.jl [5], for computing the critical points of $g(x)$. The ODE solver from DifferentialEquations.jl [14] is used for the Morse flows (9) between critical points.

To find critical points, we solve the rational function equations $\nabla\,{\rm log}(g(x))=0$ in (4). These are $n$ equations in $n$ variables $x_{1},\ldots,x_{n}$ and $k+1$ parameters $u=(s_{1},\ldots,s_{k},t)$. For the monodromy, we need a start pair $(x_{0},u_{0})$ such that $x_{0}$ is a solution for $\nabla\,{\rm log}(g(x))=0$ with parameters $u_{0}$. Since $\nabla\,{\rm log}(g(x))$ is linear in $u$, we can use linear algebra to get such a start pair. Indeed, we can sample a random point $x_{0}$ and compute $u_{0}$ by solving linear equations. However, this works only if $k+1>n$, and if $\nabla\,{\rm log}(g(x_{0}))$ has full rank. To avoid these issues, we use the following trick. If $k+1\leq n$, we define new auxiliary parameters $v=(v_{1},\ldots,v_{n-k})$. Using monodromy, we solve the system $\nabla\,{\rm log}(g(x))-(Au+Bv)=0$, where $A\in\mathbb{C}^{n\times(k+1)}$ and $B\in\mathbb{C}^{n\times(n-k)}$ are random matrices. Afterwards, we track the homotopy $\nabla\,{\rm log}(g(x))-\lambda\cdot(Au+Bv)=0$ from $\lambda=1$ to $\lambda=0$ using the solutions we have just computed. The parameter continuation theorem [15, 4] asserts that this approach works.

We modified the default options in HomotopyContinuation.jl for the monodromy loops. In HypersurfaceRegions.jl, a computation finishes if the last $10$ new monodromy loops have not provided any new solutions. The number $10$ is a conservative parameter here. Its purpose is to increase the probability of finding all critical points. The tolerance parameter $\delta$ for bounded regions, as described in the end of Section 3, is set to the default value $\delta=10^{-5}$.

This final section is dedicated to systematic experiments with our software. We experimented by running HypersurfaceRegions.jl on random instances, drawn from two classes:

1. 1.

   arrangements defined by random polynomials;

2. 2.

   arrangements defined by random spectrahedra.

The first class is self-explanatory. To create (1), we choose $f_{i}$ at random from some probability distribution on the space of inhomogeneous polynomials of degree $d$ in $n$ variables.

The second class requires an explanation. We fix positive integers $n$ and $m$, and we draw $n+1$ samples $A_{0},A_{1},\ldots,A_{n}$ from the $\binom{m+1}{2}$-dimensional space of symmetric $m\times m$ matrices. We consider the arrangement defined by the $k=2^{m}-1$ principal minors of the matrix

The distinguished region where all $k$ principal minors are positive is the spectrahedron.

The stratification of symmetric matrices by signs of principal minors was studied by Boege, Selover and Zubkov in [3]. We explore random low-dimensional slices of the Boege-Selover-Zubkov stratification. How many regions get created, and what is their topology?

We now describe our experiments with arrangements of generic hypersurfaces. We run HypersurfaceRegions.jl on $k$ polynomials in $n$ variables of degree $d=2$ and $d=3$, where the coefficients are drawn independently from the standard Gaussian. Each experiment is carried out $N$ times, where $N$ is inverse proportional to the running time for each instance.

Our findings are presented in Table 1 for $d=2$ and in Table 2 for $d=3$. For each row we fix the parameters $n$ and $k$. The time is the average running time per instance. This is measured in seconds. The fourth column concerns the total number of regions. We report the minimum number and the maximum number across all instances. The fifth column similarly reports the range for the number of realizable sign vectors $\sigma$. The sixth column gives the maximal number of regions per fixed sign pattern. And, finally, in the last column we report the minimum and maximum observed for the Euler characteristics $\chi$ of the regions.

Let us discuss the first row in Table 1. The code runs $15.5$ seconds, and it identifies $\rho$ distinct regions in $\mathbb{R}^{3}$, where $46\leq\rho\leq 325$. The number of realizable sign patterns is at most $196$. This is less than the number $2^{k}=2^{8}=256$ of all sign patterns. There were up to $11$ regions per sign pattern. The Euler characteristic of any region was between $-5$ and $2$.

Table 2 presents analogous results for cubics in $\mathbb{R}^{n}$ where $n=3,4,5$. The second-last row concerns quintuples of cubics in $\mathbb{R}^{4}$. It takes less than six minutes to identify all regions, of which there are up to $119$. We observed a narrow range $\{-4,\ldots,3\}$ of Euler characteristics.

Generic instances have no tangencies to the hyperplane at infinity. To experiment with that issue, we can try $k$ paraboloids in $\mathbb{R}^{n}$. Each paraboloid contributes one undecided region touching the hyperplane at infinity. Here is an instance of $k\!=\!4$ paraboloids for $n\!=\!3$:

The output in Figure 6 consists of six regions, each with a unique sign pattern. Notice that the label “undecided” does not mean we claim this region touches the hyperplane at infinity. It means that our algorithm could not decide whether this region is bounded or unbounded.

We now turn to random spectrahedra. We sampled symmetric $m\times m$ matrices $A_{0},\ldots,A_{n}$ where the entries are independent standard Gaussians. The principal minors of the matrix $A(x)$ define an arrangement of $k=2^{m}-1$ hypersurfaces in the affine space $\mathbb{R}^{n}$. We ran HypersurfaceRegions.jl on this input. Our results are summarized in Table 3. For each row we fix the parameters $n$ and $m$. The columns have the same meaning as before. For instance, the fifth column reports the range for the number of realizable sign vectors $\sigma$.