Computing the volume of compact semi-algebraic sets
Let be a compact basic semi-algebraic set defined as the real solution set of multivariate polynomial inequalities with rational coefficients. We design an algorithm which takes as input a polynomial system defining and an integer and returns the -dimensional volume of at absolute precision .
Our algorithm relies on the relationship between volumes of semi-algebraic sets and periods of rational integrals. It makes use of algorithms computing the Picard-Fuchs differential equation of appropriate periods, properties of critical points, and high-precision numerical integration of differential equations.
The algorithm runs in essentially linear time with respect to . This improves upon the previous exponential bounds obtained by Monte-Carlo or moment-based methods. Assuming a conjecture of Dimca, the arithmetic cost of the algebraic subroutines for computing Picard-Fuchs equations and critical points is singly exponential in and polynomial in the maximum degree of the input.
Semi-algebraic sets are the subsets of which are finite unions of real solution sets to polynomial systems of equations and inequalities with coefficients in . Starting from Tarski’s algorithm for quantifier elimination (Tarski, 1998) improved by Collins through the Cylindrical Algebraic Decomposition algorithm (Collins, 1975), effective real algebraic geometry yields numerous algorithmic innovations and asymptotically faster routines for problems like deciding the emptiness of semi-algebraic sets, answering connectivity queries or computing Betti numbers (e.g., Basu et al., 2006; Safey El Din and Schost, 2003; Canny, 1988; Safey El Din and Schost, 2017; Bürgisser et al., 2019). The output of all these algorithms is algebraic in nature. In this paper, we study the problem of computing the volume of a (basic) compact semi-algebraic set defined over . The output may be transcendental: for instance, the area of the unit circle in is .
Volumes of semi-algebraic sets actually lie in a special class of real numbers, for they are closely related to Kontsevich-Zagier periods introduced in (Kontsevich and Zagier, 2001). A (real) period is the value of an absolutely convergent integral of a rational function with rational coefficients over a semi-algebraic set defined by polynomials with rational coefficients. For example, algebraic numbers are periods, as are , , . Since , volumes of semi-algebraic sets defined over are periods. Conversely, interpreting an integral as a “volume under a graph” shows that periods are differences of volumes of semi-algebraic sets defined over . In (Viu-Sos, 9 03), it is further shown that periods are differences of volumes of compact semi-algebraic sets defined over .
The problem we consider in this paper is thus a basic instance of the more general problem of integrating an algebraic function over a semi-algebraic set; it finds applications in numerous areas of engineering sciences. Performing these computations at high precision (hundreds to thousands of digits) is also relevant in experimental mathematics, especially for discovering formulas, as explained, for example, in (Bailey and Borwein, 2011). Most of the examples featured in this reference are periods, sometimes in disguise.
The simplest semi-algebraic sets one can consider are polytopes. The computation of their volume has been extensively studied, with a focus on the complexity with respect to the dimension. It is known that even approximating the volume of a polytope deterministically is #P-hard (Dyer and Frieze, 1988; Khachiyan, 1989). The celebrated probabilistic approximation algorithm in (Dyer et al., 1991), which applies to more general convex sets, computes an -approximation in time polynomial in the dimension of the set and . A key ingredient for this algorithm is a Monte Carlo method for efficiently sampling points from a convex set. Since then, Monte Carlo schemes have been adopted as the framework of several volume estimation algorithms.
In contrast, we deal here with compact semi-algebraic sets which can be non-convex and even non-connected. Additionally, while volumes of polytopes are rational, the arithmetic nature of volumes of semi-algebraic sets is much less clear, as unclear as the nature of periods. This raises the question of the computational complexity of a volume, even taken as a single real number.
A simple Monte Carlo technique applies in our setting as well: one samples points uniformly in a box containing and estimates the probability that they lie in . This method is of practical interest at low precision but requires samples to achieve an error bounded by with high probability. We refer to (Koiran, 1995) which deals with definable sets, a class which encompasses semi-algebraic sets.
In a different direction, numerical approximation schemes based on the moment problem and semi-definite programming have been designed in (Henrion et al., 2009). They are also of practical interest at low precision, and can provide rigorous error bounds, but the convergence is worse than exponential with respect to (Korda and Henrion, 2018).
Another line of research, going back to the nineteenth century, is concerned with the computation of periods of algebraic varieties. In particular, we build on work (Chudnovsky and Chudnovsky, 1990) on the high-precision numerical solution of ODEs with polynomial coefficients which was motivated, among other things, by applications to periods of Abelian integrals (see Chudnovsky and Chudnovsky, 1990, p. 133).
We describe a new strategy for computing volumes of semi-algebraic sets, at the crossroads of effective algebraic and real algebraic geometry, symbolic integration, and rigorous numerical computing. Our approach effectively reduces the volume computation to the setting of (Chudnovsky and Chudnovsky, 1990). It yields an algorithm that approximates the volume of a fixed, bounded basic semi-algebraic set in almost linear time with respect to the precision. More precisely, we prove the following bit complexity estimate.
Theorem 1 ().
Let be polynomials in , and let be the semi-algebraic set defined by . Assume that is compact. There exists an algorithm which computes, on input and , an approximation of the volume of with . When are fixed, the algorithm runs in time (for any ) as .
The algorithm recursively computes integrals of volumes of sections of . Let denote the -dimensional volume of , for some nonzero linear projection . In our setting, is a piecewise analytic function and, except at finitely many , is solution of a linear differential equation with polynomial coefficients known as a Picard-Fuchs equation.
The problem points belong to the critical locus of the restriction of the projection to a certain hypersurface containing the boundary of and are found by solving appropriate polynomial systems. (Compare (Khachiyan, 1993) in the case of polytopes.) The Picard-Fuchs equation for is produced by algorithms from symbolic integration, in particular (Bostan et al., 2013; Lairez, 2016). To obtain the volume of , it then suffices to compute with a rigorous numerical ODE solving algorithm, starting from values at suitable points obtained through recursive calls.
The complexity with respect to the dimension of the ambient space and the number , maximum degree , and coefficient size of the polynomials is harder to analyze. We will see, though, that under reasonable assumptions, the “algebraic” steps (computing the critical loci and of the Picard-Fuchs equations) take at most arithmetic operations in .
The idea of the method is well illustrated by the example of a torus , here of major radius and minor radius . Let
The area (2-dimensional volume) of a section defines a function (see Figure 1).
It is analytic, except maybe at the critical values and where the real locus of the curve is singular. On each interval on which is analytic, it satisfies the Picard-Fuchs equation
We know some special values of , namely , and . Additionally, we have as . These properties characterize the analytic function in the -dimensional space of analytic solutions of the differential equation (1) on , and similarly for . (Our algorithm actually uses recursive calls at generic points instead of these ad hoc conditions.) The rigorous ODE solver part of the Sage package ore_algebra (Mezzarobba, 2016) determines in less than a second that
And indeed, it is not hard to see in this case that . We can obtain 1000 digits in less than a minute.
The remainder of this article is organized as follows. In Section 2, we give a high-level description of the main algorithm. As sketched above, the algorithm relies on the computation of critical points, Picard-Fuchs equations, and numerical solutions of these equations. In Section 3, we discuss the computation of Picard-Fuchs equations and critical points, relating these objects with analyticity properties of the “section volume” function. Then, in Section 4, we describe the numerical solution process and study its complexity with respect to the precision. Finally, in Section 5, we conclude the proof of Theorem 1 and state partial results on the complexity of the algorithm with respect to , , and .
We would like to thank the anonymous reviewers for their careful reading and valuable comments.
2. Volumes of semi-algebraic sets
We start by designing an algorithm which deals with the case of a union of connected components of a semi-algebraic set defined by a single inequality. Next, we will use a deformation technique to handle semi-algebraic sets defined by several inequalities.
2.1. Sets defined by a single inequality
Let and be the semi-algebraic set
Let be the projection on the -coordinate. We want to compute the volume of a union of connected components of starting from the volumes of suitable fibers . For technical reasons, we first consider the slightly more general situation where is a union of connected components of for some open interval . From a computational point of view, we assume that is described by a semi-algebraic formula , that is,
where is a finite disjunction of conjunctions of polynomial inequalities with (in our setting) rational coefficients.
For , let and . Let (we will often omit the subscript ) be the set of exceptional values
Thus, when is square-free, exceptional values are either critical values of the restriction to the hypersurface of the projection , or images of singular points of . By definition of , for any , the zero set of is a smooth submanifold of .
Further, we say that assumption (R) holds for if
The mainstay of the method is the next result, to be proved in §3. Let denote the set of Fuchsian linear differential operators with coefficients in whose local exponents at singular points are rational (see §4 for reminders on Fuchsian operators and their exponents).
Theorem 2 ().
If is bounded and , then the function is solution of a computable differential equation of the form , where depends only on .
We will also use the following proposition, which summarizes the results of Proposition 12 and Lemma 13 in §4. The complete definition of “good initial conditions” is given there as well. Up to technical details, this simply means a system of linear equations of the form that suffices to characterize a particular solution among the solutions of . An -approximation of is made of the same equations with each right-hand side replaced by an enclosure of diameter .
Proposition 2 ().
Let have order , and let be a real interval with algebraic endpoints. Let be a solution of with a finite limit at and be a system of good initial conditions for on defining .
Given , a precision and a -approximation of , one can compute an interval of width (as for fixed , , and ) containing .
Given , one can compute such that the form a system of good initial conditions for on .
Assume now that is a bounded union of connected components of (i.e., that we can take above), and that (R) holds for . The algorithm is recursive. Starting with input , and , it first computes the set of exceptional values so as to decompose into intervals over which the function satisfies the differential equation given by Theorem 2. Since is bounded, one has
Fix and consider the interval . Since is annihilated by , its anti-derivative vanishing at is annihilated by the operator , which belongs to since does. Additionally, if is a system of good initial conditions for that defines , then is a system of good initial conditions for defining (see Lemma 12 in §4). Thus, by Proposition 2, to compute to absolute precision , it suffices to compute , , to precision .
By definition of , since , there is no solution to the system
which means that (R) holds for . Additionally, is a bounded union of connected components of . Hence, the values can be obtained by recursive calls to the algorithm with instantiated to .
The process terminates since each recursive call handles one less variable. In the base case, we are left with the problem of computing the length of a union of real intervals encoded by a semi-algebraic formula. This is classically done using basic univariate polynomial arithmetic and real root isolation (Basu et al., 2006, Chap. 10).
The complete procedure is formalized in Algorithm 1. The quantities denoted with a tilde in the pseudo-code are understood to be represented by intervals, and the operations involving them follow the semantics of interval arithmetic. Additionally, we assume that we have at our disposal the following subroutines:
, which returns an encoding for a finite set of real algebraic numbers containing the exceptional values associated to , sorted in increasing order;
where and is a semi-algebraic formula describing a union of connected components of , which returns an interval of width containing .
The following result summarizes the above discussion.
2.2. Sets defined by several inequalities
Now, we show how to compute the volume of a basic semi-algebraic set defined by
assuming that is compact.
We set , and consider the semi-algebraic set defined by . Observe that the polynomial satisfies (R) because . We can hence choose an interval with that contains no element of . Let and be the projection on the -coordinate. For fixed , the set can be viewed as a bounded subset of , whose volume tends to as .
The set itself is bounded and the formula
defines in . In addition, is a union of connected components of . Indeed, for any point with , it holds that . This implies that where is the interior of . Therefore, is both relatively closed (as the trace of ) and open (as that of ) in .
We are hence in the setting of the previous subsection. Since by definition of , Theorem 2 applies, and the function is annihilated by an operator which is computed using the routine PicardFuchs introduced earlier. By Proposition 2, one can choose rational points such that the values of at these points characterize it among the solutions of , and, given sufficiently precise approximations of , one can compute to any desired accuracy.
The “initial conditions” are computed by calls to Algorithm 1 with and specialized to . In the notation of §2.1, this corresponds to taking and . Thus, is compact, and, since no can change sign on a connected component of for , it is the union of those connected components of where . Additionally, (R) holds for since . Therefore, the assumptions of Theorem 3 are satisfied.
Theorem 4 ().
Let . Let be the semi-algebraic set defined by . Assume that is bounded. Then, given and a working precision , Algorithm 2 (Volume) computes an interval containing of width as for fixed
Remark 5 ().
In case has empty interior, Algorithm 2 returns zero. When is contained in a linear subspace of dimension , one could in principle obtain the -volume of by computing linear equations defining the subspace (using quantifier elimination as in (Khachiyan and Porkolab, 2000; Safey El Din and Zhi, 2010)) and eliminating variables. The new system would in general have algebraic instead of rational coefficients, though.
Lastly, we note that a more direct symbolic computation of integrals on general semi-algebraic sets depending on a parameter is possible with Oaku’s algorithm (Oaku, 2013), based on the effective theory of -modules.
3. Periods depending on a parameter
Let us now discuss in more detail the main black boxes used by the volume computation algorithm. In this section, we study how the volume of a section varies with the parameter .
3.1. Picard-Fuchs equations
Let be a rational function. A period of the parameter-dependent rational integral is an analytic function , for some open subset of or such that for any there is an -cycle and a neighborhood of such that for any , is disjoint from the poles of and
Recall that an -cycle is a compact -dimensional real submanifold of and that such an integral is invariant under a continuous deformation of the integration domain as long as it stays away from the poles of , as a consequence of Stokes’ theorem. It is also well known that such a function depends analytically on , by Morera’s theorem for example.
For instance, algebraic functions are periods: if satisfies a nontrivial relation , with square-free , then is a period by the residue theorem applied to
where encloses and no other root of . Indeed, the integrand decomposes as , where the functions parametrize the roots of , and, w.l.o.g., .
Periods of rational functions are solutions of Fuchsian linear differential equations with polynomial coefficients known as Picard-Fuchs equations. This was proved in (Picard, 1902) in the case of three variables at most and a parameter and generalized later, using either the finiteness of the algebraic De Rham cohomology (e.g. Grothendieck, 1966; Monsky, 1972; Christol, 1985) or the theory of D-finite functions (Lipshitz, 1988). The regularity of Picard-Fuchs equations is due to Griffiths (see Katz, 1971).
Theorem 6 ().
If is the period of a rational integral then is solution of a nontrivial linear differential equation with polynomial coefficients , where the operator belongs to the class introduced in §2.1.
Theorem 7 ((Bostan et al., 2013)).
A period of the form (3) is solution of a differential equation of order at most where is the degree of ; and one can compute such an equation in operations in .
Note however that the algorithm underlying this result might not return the equation of minimal order, but rather a left multiple of the Picard-Fuchs equation. So there is no guarantee that the computed operator belongs to . On the other hand, Lairez’s algorithm (Lairez, 2016) can compute a sequence of operators with non-increasing order which eventually stabilizes to the minimal order operator. In particular, as long as the computed operator is not in , we can compute the next one, with the guarantee that this procedure terminates. A conjecture of Dimca (Dimca, 1991) ensures that it terminates after at most steps, leading to a complexity bound as in Theorem 7.
3.2. Volume of a section and proof of Theorem 2
Theorem 8 ().
If and if is bounded then the function is a period of the rational integral
Let . By Stokes’ formula,
where is the boundary of . Due to the regularity assumption , the gradient does not vanish on the real zero locus of , denoted . Because is a union of connected components of , it follows that is a compact -dimensional submanifold of contained in .
For , let be the Leray tube defined by
This is an -dimensional submanifold of . We choose small enough that : this is possible because does not vanish on which is compact.
Let ; observe that does not cancel the denominator of . Leray’s residue theorem (Leray, 1959) shows that
(In Pham’s (Pham, 2011, Thm. III.2.4) notation, we have , , , and .)
To match the definition of a period and conclude the proof, it is enough to prove that, locally, the integration domain can be made independent of . And indeed, since is a union of connected components of , we have . Therefore, since is connected and , the restriction of the projection defines a submersive map from onto . Additionally, is compact, hence this map is proper. Ehresmann’s theorem then implies that there exists a continuous map such that induces a homeomorphism for any . In particular, we have
This formulation makes it clear that deforms continuously into as varies. Since does not intersect the polar locus of , neither does when and are close enough, by compactness of and continuity of the deformation. Therefore, given any , we have for close enough to . ∎
The choice of as a primitive of in Theorem 8 is arbitrary, but of little consequence, since the Picard-Fuchs equation only depends on the cohomology class of the integrand.
3.3. Critical values
Theorem 2 does not guarantee that satisfies the Picard-Fuchs equation on the whole domain where the equation is nonsingular. It could happen that the solutions extend analytically across an exceptional point, or that some of them have singularities between two consecutive exceptional points. As a consequence, we need to explicitly compute .
Lemma 9 ().
There exists an algorithm which, given on input a polynomial of degree satisfying (R), computes a polynomial of degree whose set of real roots contains , using operations in .
Recall that, when (R) holds, the set is finite. Our goal is to write as the root set of a univariate polynomial . Consider the polynomial We start by computing at least one point in each connected component of the real algebraic set defined by using (Basu et al., 2006, Algorithm 13.3). By (Basu et al., 2006, Theorem 13.22), this algorithm uses operations. It returns a rational parametrization: polynomials , , in of degree such that is square-free and the set of points
meets every connected component of the zero set of . In particular, . As a polynomial , we take the resultant with respect to of and : its set of roots contains . Since and have degree , this last step also uses operations in (von zur Gathen and Gerhard, 1999). ∎
Let us turn to the numerical part of the main algorithm. It is known (Chudnovsky and Chudnovsky, 1990; van der Hoeven, 2001) that Fuchsian differential equations with coefficients in can be solved numerically in quasi-linear time w.r.t. the precision. Yet, some minor technical points must be addressed to apply the results of the literature to our setting. We start with reminders on the theory of linear ODEs in the complex domain (e.g. Poole, 1936; Hille, 1976). Consider a linear differential operator
of order with coefficients in .
Recall that is a singular point of when the leading coefficient of vanishes at . A point that is not a singular point is called ordinary. Singular points are traditionally classified in two categories: a singular point is a regular singular point of if, for , its multiplicity as a pole of is at most , and an irregular singular point otherwise. The point at infinity in is said to be ordinary, singular, etc., depending on the nature of after the change of variable . An operator with no irregular singular point in is called Fuchsian.
Fix a simply connected domain containing only ordinary points of , and let be the space of analytic solutions of the differential equation . According to the Cauchy existence theorem for linear analytic ODEs, is a complex vector space of dimension . A particular solution is determined by the initial values at any point .
At a singular point, there may not be any nonzero analytic solution. Yet, if is a regular singular point, the differential equation still admits linearly independent solutions defined in the slit disk for small enough and each of the form
where , , and for exactly one (Poole, 1936, §16). The functions are analytic for (including at ). The algebraic numbers are called the exponents of at .
Suppose now that is either an ordinary point of lying in the topological closure of , or a regular singular point of situated on the boundary of . As a result of the previous discussion, we can choose a distinguished basis of in which each is characterized by the leading monomial111 More precisely, denoting in (5), there are computable pairs such that, for all , we have , for , and whenever . of its local expansion (5) at . At an ordinary point for instance, the coefficients of the decomposition of a solution on are , that is, essentially the classical initial values. Observe that when no two exponents have the same imaginary part, the elements of all have distinct asymptotic behaviours as . In particular, at most one of them tends to a nonzero finite limit. As Picard-Fuchs operators have real exponents according to Theorem 6, this observation applies to them.
Let be a second point subject to the same restrictions as . Let be the transformation matrix from to . The key to the quasi-linear complexity of our algorithm is that the entries of this matrix can be computed efficiently, by solving the ODE with a Taylor method in which sums of Taylor series are computed by binary splitting (Beeler et al., 1972, item 178), (Chudnovsky and Chudnovsky, 1990). The exact result we require is due to van der Hoeven (van der Hoeven, 2001, Theorems 2.4 and 4.1); see also (Mezzarobba, 2011) for a detailed algorithm and some further refinements. Denote by the complexity of -bit integer multiplication.
Theorem 10 ((van der Hoeven, 2001)).
For a fixed operator and fixed algebraic numbers as above, one can compute the matrix with an entry-wise error bounded by in operations.
Since is linear, this result suffices to implement the procedure DSolve required by the main algorithm. More precisely, suppose that returns a matrix of complex intervals of width that encloses entry-wise.
Definition 11 ().
A system of good initial conditions for on , denoted , is a finite family of pairs where and is a linear form that belongs to the dual basis of for some algebraic point (which may depend on ), with the property that span the dual space of .
A system of good initial conditions on is a system of good initial conditions on for some .
In other words, a system of good initial conditions is a choice of coefficients of local decompositions of a solution of whose values determine at most one solution, and of prescribed values for these coefficients. When the system is compatible, we say that it defines the unique solution of that satisfies all the constraints. Let us note in passing the following fact, which was used in §2.2.
Lemma 12 ().
Let be ordinary points of such that is a system of good initial conditions for on , and let . Then is a system of good initial conditions for on .
The derivative maps the solution space of to that of , and its kernel consists exactly of the constant functions. By assumption, a solution of is completely defined by its values at , hence a solution of is characterized by the values of its derivative at the same points, along with its limit at . Because has order at least (otherwise, would not be a system of good initial conditions), the conditions are of the form with belonging to the dual basis of some , as required. So is the condition since has solutions with a nonzero finite limit at . ∎
Algorithm 3 evaluates the solution of an operator given by a system of good initial conditions. Note that the algorithm is allowed to fail. It fails if the intervals are not accurate enough for the linear algebra step on line 10 to succeed, or if the linear system, which is in general over-determined, has no solution. The following proposition assumes a large enough working precision to ensure that this does not happen. Additionally, we only require that the output be accurate to within , so as to absorb any loss of precision resulting from numerical stability issues or from the use of interval arithmetic.
Proposition 12 ().
Suppose that the operator is Fuchsian with real exponents. Let be real algebraic numbers, and let be a real analytic solution of on the interval such that tends to a finite limit as . Let be a system of good initial conditions for on that defines .
Given the operator , the point , a large enough working precision , and an approximation of where is an interval of width at most containing , (Algorithm 3) computes a real interval of width containing in time .
At the end of the loop, we have , , and the entries of are intervals of width . The coefficients of the decomposition of in the basis satisfy for all , where is the th row of . As is a system of good initial conditions, the linear system has no other solution. Step 10 hence succeeds in solving the interval version as soon as the and the entries of the are thin enough intervals. It then returns intervals of width .
We assumed that tends to a finite limit at . It follows that the decomposition of on only involves the basis elements with a finite limit at . Either contains an element that tends to , in which case , or every solution that converges tends to zero, and then the limit is zero. Since, by assumption, is real, we can ignore the imaginary part of the computed value. In both cases, the algorithm, when it succeeds, returns a real interval of width containing .
As for the complexity analysis, all including are algebraic, hence Theorem 10 applies and shows that each call to TransitionMatrix runs in time . The matrix multiplications at step 8 take operations. The cost of solving the linear system (which is of bounded size) is as well. The cost of the remaining steps is independent of