12 Quadrature
12.1 An Overview of Quadrature
Numerical integration is a basic task which occurs ubiquitously through Flex. There are many functions to integrate, with some being more structured or well-behaved than others. To this end, we need a bevy of general purpose algorithms: there is no one-size fits all solution.
An alternative approach to numerical integration is symbolic integration, along the lines of what might be found in sympy or Mathematica. Sometimes symbolic integration is an appropriate tool, but it is all too easy to find examples of integrands which cause even the most sophisticated computer algebra system to give up.
Instead, numerical integration is nearly always accomplished by means of an (order ) quadrature rule, consisting of pairs of points (sometimes called nodes) and weights used to reduce an integral to a sum: Here, is the quadrature error, depending on both the function being integrated and the quadrature rule. There are a wide variety of classical quadrature rules available from standard textbooks.
12.2 Gauss Quadrature
The workhorse of numerical integration is the Guass-Legendre quadrature rule. The th Gauss-Legendre rule approximates integrals of the form . The points of this quadrature rule are the zeros of the th Legendre polynomial, and the weights are the corresponding Newton-Cotes weights. That is, if is the th Lagrange polynomial such that , where are the Gauss-Legendre nodes, then: Gauss-Legendre rules are available for essentially any . Historically, the Golub-Welsch algorithm was used for rules up to moderate . More recently, different researchers have used asymptotics to estimate the node positions and have subsequently "polished" these rules by using them as a warm start for the moment-matching equations. These approaches have allowed for the computation of Gauss-Legendre rules with in the millions.
The th Gauss-Legendre rule has the remarkable property that it integrates all polynomials up to degree . That is, it integrates a polynomial with as many coefficients as the rule has degrees of freedom ( nodes and weights). Since the Gauss-Legendre rule can be easily shifted to any other interval under an affine transformation without altering its accuracy, it can be used as the basis of a composite rule in which an interval is adaptively refined to integrate a function that is not smooth enough to integrate with a high-order Gauss-Legendre rule directly.
12.3 Tensor-product Rules
Several one-dimensional quadrature rules can be combined to form tensor product rules in the natural way. For example, in two dimensions, if we have a quadrature rule with nodes and weights which approximates integrals of the form and another rule with nodes and weights which can be used to numerically integrate , then the tensor product rule formed from their combination looks like: Tensor-product Gauss-Legendre rules are a natural choice when it comes to integrating smooth functions on rectangles in multiple dimensions.
Although one-dimensional Gauss-Legendre rules are arguably the best choice of quadrature rule on a closed and bounded interval, this is not the case in two or more dimensions. There are numerous examples of quadrature rules which have even fewer nodes but which have essentially the same accuracy for certain special geometries, such as the triangle, the disk, the sphere, and so on. Nevertheless, because of their efficient tensor-product structure, they can be competitive with these rules despite their suboptimal node count.
12.4 Mapped Quadrature Rules
Building quadrature rules for complicated domains which have no known rules available analytically is a task of crucial importance in simulation. In this section and the next we will examine two strategies which often work in tandem. Useful in their own right, these "mapped composite rules" will also be fed as input to some of our later quadrature algorithms.
The first strategy works for domains which are the image of a simple reference domain under a nonlinear mapping. In this case, if we have a rule for a reference domain with points and weights , and we would like to integrate a function over a domain which is the image of under a nonlinear mapping , then a simple application of integration by substitution gives: The resulting "mapped quadrature rule" consists of the nodes and weights , where .
It is crucially important to note that a mapped rule which involves a nonlinear mapping will invariably be less accurate than the original rule. For example, if we start with a two-dimensional tensor-product Gauss-Legendre rule with nodes in each direction, then this rule integrates all monomials of the form where and . If the mapping is affine (i.e., a degree 1 polynomial in each component), then the same set of monomials is integrated over as well as itself. On the other hand, say is quadratic in each of its components. Then the subsequent mapped rule will be able to integrate monomials with and such that and at most. Since the quadrature error depends on how well we can approximate the integrand with a polynomial in the space being integrated, this is a serious reduction in accuracy.
12.5 Composite Quadrature Rules
A second strategy for constructing a quadrature rule for a complicated domain is to build a mesh of that domain and then to concatenate simple rules for each mesh element together into a composite quadrature rule. These simple rules will invariably be mapped rules.
As an example, assume that is a complicated domain and that we have a mesh of triangles which conforms exactly to the boundary of . Let be the th such triangle and let be an affine map that sends the standard triangle to . If we have a quadrature rule with nodes and weights () then we can approximate the integral of a function on as: Consequently, we can define a composite rule on with nodes and weights .
12.6 Matching Moments
Assume that we can expand a function over a domain in a series: where is a multi-index. If we have a quadrature rule with nodes and weights that integrate all monomials of total degree at most (i.e., monomials where ), then the quadrature error for is exactly:
In the context of finite elements, where we might have a smooth and an which satisfies , we can reasonably expect the quadrature error to behave like .
Working with the moment-matching equations assumes we have the ability to set up the right-hand side, consisting of the integrals . This almost certainly requires another quadrature rule. In IGA, where three-dimensional domains are primarily of interest, integrals of simple functions like monomials can be pushed to the boundary using the divergence theorem and approximated using a surface quadrature rule.
From this rough analysis, we can focus our efforts on finding quadrature rules which integrate certain polynomials or spaces of functions. This is a useful perspective to adopt when constructing quadrature rules on unstructured domains. Say our goal is to find a quadrature rule which integrates all bivariate polynomials of total degree at most . All such polynomials can be written as a linear combination of monomials of the form where . As a result, we can focus our efforts on finding a rule which integrates each monomial basis function individually. If we have a rule with nodes and weights , the set of equations which enforces this condition is: where is the total number of monomials to be integrated. These nonlinear equations in unknowns are known as the moment-matching equations. Note that the moment-matching equations are not particular to the monomial basis or even vector spaces of polynomials.
If the moment-matching equations consist of fewer unknowns than equations, it is not possible to solve them in general. Consequently, we can anticipate the requirement that . This inequality is achieved by Gauss-Legendre rules in 1D. We will define the efficiency of a quadrature rule to be: Since , we want to keep in mind the task of finding quadrature rules with an efficiency as close to 1 as possible. In dimensions two and higher, very little is known about efficient quadratures. It may not always be possible to achieve the minimum value of . There are some bizarre cases where the minimal quadrature for a given domain and polynomial space can only be achieved by placing nodes outside the domain! Quadrature rules which do achieve the minimum are called Gaussian or Gauss rules. Consequently, and somewhat confusingly, tensor-product Gauss-Legendre rules are not Gauss rules.
There are two approaches to using the moment-matching equations as a numerical tool for constructing quadrature rules. In the first, we solve the full, nonlinear system of equations using Gauss-Newton. Define the residual function: Then, Gauss-Newton solves the nonlinear least squares problem: If we let denote the vector of weights and nodes at the step, then the Gauss-Newton iteration is for . If is sufficiently close to a root of the moment-matching equations, then convergence will be quadratic.
The approach of solving the moment-matching equations as a nonlinear system is one in which both the nodes and weights are allowed to vary, and presupposes an inexact quadrature rule which is a good approximation of an exact quadrature (i.e., a root of the system). If the initial quadrature is far from a zero of the moment-matching equations, then the sensitivity of Newton’s method means that we will need to employ a more sophisticated algorithm for nonlinear optimization, such as sequential quadratic programming coupled with a line search. These algorithms are more expensive and unpredictable in terms of their runtime performance. For this reason, we do not solve the full moment-matching equations unless we can be confident that our initial quadrature is reasonably accurate.
Next, for the second approach, we fix the node positions and allow only the weights to vary. The result is a linear system of equations to solve. If there are more nodes than basis functions (), then the linear moment-matching equations will be underdetermined. Basic linear algebra tells us that an underdetermined linear system has an infinite number of solutions. Since our goal is to compute a quadrature rule with the minimum possible number of nodes, we could attempt to compute a sparse solution to the linear moment matching equations, in which case not all are nonzero. We will now expand on this idea a bit.
Recall that for univariate interpolation, the Vandermonde matrix is nonsingular if and only if the interpolation nodes are distinct. The same intuition applies to a multidimensional Vandermonde matrix with more columns than rows, although it is captured by the concept of "unisolvence". Roughly speaking, we can expect columns corresponding to the nodes that are in "general position" to be linearly independent. This suggests that a reasonable approach might be to greedily select a subset of quadrature nodes which are "good"; furthermore, we can hypothesize that we should not need more than nodes, which would lead to a quadrature rule with . As an example, a simple approach would be to run rank-revealing QR and terminate once we drive the residual sufficiently close to zero, although this might result in a rule with negative weights.
12.7 Nonnegative Least Squares Quadrature
A commonly desired feature of a quadrature rule is that for each —
For an interpolatory quadrature rule, . This inequality is attained exactly when the quadrature weights are positive. So, in some sense, the best stability that can be achieved (as measured by this metric!) is attainable only by positive quadrature rules. Nevertheless, finding a quadrature rule with both positive and negative weights for which where , say, is the common case rather than a rarity. Nevertheless, in this section, we consider how to explicitly compute quadrature rules with positive weights by solving the linear moment-matching equations.
We will assume that we have a fixed set of quadrature nodes and that we would like to integrate a set of linearly independent basis functions , where . Define the system matrix (a type of generalized, multivariate Vandermonde matrix) and right-hand side for the moment-matching equations by: If we write the vector of quadrature weights as , then the linear moment-matching equations are just: If , this linear system is underdetermined and has an infinite number of solutions.
Bearing in mind the fact that the moment-matching equations are likely to be underdetermined as well as our goal of producing a quadrature rule with positive weights, it’s natural to pose this problem in the form of a constrained minimization problem: Let’s recall some basic facts about the geometry of this optimization problem. First, for the unconstrained problem (i.e., find such that is minimized), since the system is underdetermined, if has a continuum of solutions if it has any solution. Hence, the set of global optima for the unconstrained minimization problem consists of a linear subspace of ; that is, a hyperplane. At the same time, the set of such that is just the nonnegative orthant in . The solution set for the constrained problem, then, is the intersection of the nonnegative orthant and the set of solutions of . So, we can see that in the general case, the constrained problem will also fail to have a unique solution: it will have a continuum of solutions unless the unconstrained solution set and the nonnegative orthant intersect exactly at the origin, which will not happen for geometry with nonzero measure.
There are lots of possible algorithms that can be used to solve equation . Since the solution to equation is nonunique, the choice of algorithm can easily bias the solution to have certain properties. For example, a simple coordinate relaxation method (gradient descent applied to one coordinate at a time), will converge reasonably quickly but will produce a solution with all weights nonzero. Indeed, equation is a positive definite quadratic program (QP) with a simple nonnegativity constraint, and can be solved with off-the-shelf QP solvers. However, a classical algorithm for solving equation due to Lawson and Hanson is particularly well-suited to nonnegative linear least squares problems of exactly this form. If it is given an initial guess of , it will relax one component of at a time and make rapid progress to a quadrature rule with a minimal number of nonzero components. Additionally, if implemented correctly, it is extraordinarily stable and accurate. This is the algorithm used internally to solve problems of this form.
Now that we know how to solve equation , all that remains is to choose an appropriate node distribution. One basic approach is to lay down a uniform grid of nodes on ’s bounding box, remove all nodes which do not lie inside , and solve the resulting optimization problem. This optimization problem may not have a solution which exactly satisfies . However, for a quadrature optimization problem of this form posed in 1D, a uniform grid with for some constant is enough to guarantee the existence of a subset of nodes which correspond to a quadrature rule that exactly integrates the desired basis. More generally, the existence of a quadrature rule with nodes is guaranteed by Tchakaloff’s theorem. Naively following the prescription to nodes can lead to an unnecessarily dense set of nodes. Instead, a more efficient approach starts with a coarse grid, attempts to solve equation to a given precision; and, on failure, doubles the number of nodes along each axis, proceeding this way until the grid is fine enough.
The approach outlined in the previous paragraph comprises the "grid-doubling" strategy for quadrature optimization. A word of caution is in order: the accuracy with which equation can be solved ultimately hinges on the accuracy of the geometry, our ability to accurately check whether a point lies inside , and how accurately we compute the right-hand side using a surface quadrature (or possibly by some other means). Each of these factors affects the geometry of the minimization problem and, ultimately, our ability to solve it accurately. It is important to get each piece right.
12.8 Pruned Quadrature Rules
In the previous section, we set up the right-hand side of equation by some means and determined a set of quadrature nodes to use for the rows of separately. Another approach takes a given quadrature rule and "prunes" it using the Lawson-Hanson algorithm.
Let and be an order quadrature rule and let be an approximation of for such that . Then, we can pose equation with the given set of nodes and the right-hand side . Assume further that our starting quadrature rule is sufficiently accurate so that we can use to compute our right-hand side; that is, just apply the starting quadrature rule to each basis function. In this case, if we run Lawson-Hanson to completion, we will recompute the weights for our quadrature rule with at most of them being nonzero. Lawson-Hanson will naturally terminate with a set of nonzero weights achieving the same error as the original quadrature rule. The fact that this is guaranteed to happen is a result due to Steinitz, which says that if a positive quadrature with nodes exists, then a positive quadrature with nodes exist. The proof of this result is essentially of a linear algebraic nature and has no immediate connection with quadrature.
To summarize the point of this section: if we have, e.g., a composite rule with positive weights which achieves a reasonable accuracy, then it is straightforward to run Lawson-Hanson to find a subset of nodes and a corresponding set of positive weights, thereby "pruning" the original quadrature rule.