Contents ↑

• The Gram matrix trick
  • Why should we have expected the Gram matrix to appear?
• Gaussian orthants
  • Opening moves
  • Genz’s approach
  • Wrapping up

Volumes of high-dimensional sphere sectors

The following problem came up recently while thinking about “semantic-geometry” approaches to solving Semantle, following a question by my former colleague Victoria Ritvo.

Problem: Let $v_1, \cdots, v_n \in \mathbb{R}^d$, where $n \le d$. What is the volume of the region $$\left\{x \in S^{d-1} | \langle x, v_i \rangle \ge 0\text{ for all }1 \le i \le n\right\}?$$

Figuring out how to solve this problem was a very fun little rabbit hole that taught me several new tricks and reminded me of some I’d forgotten. Let’s dive in!

As a warm-up, we’ll note two easy facts:

• The only interesting regime is when $n \le d$. For $n > d$, a generic choice of $v_i$ produces an empty spherical sector, and any choice of $n > d$ vectors $v_i$ that does produce a nonempty sector must have redundant constraints.
  
• For reasonably large $n$, naive Monte Carlo integration will be very inefficient because almost all points on $S^{d-1}$ will not lie in the desired region.

The Gram matrix trick

The first observation we can make is that the problem is projectively symmetric: there’s no real dependence on $\|x\|$. So we could express the volume fraction of the desired region as a probability $\Pr\left[ \langle x, v_i \rangle \ge 0 | x \sim \mathcal{D}\right]$ for any spherically-symmetric distribution $\mathcal{D}$. The natural choice of spherically-symmetric distribution is the multivariate Gaussian $N(\mathbf{0}, \mathbf{I}_d)$.

From here, we arrive at a remarkable fact: if $x \sim N(\mathbf{0}, \mathbf{I}_d)$, then the random variates $Z_i = \langle x, v_i\rangle$ are themselves distributed as a multivariate Gaussian. The easiest way to argue this is to use the Cramér-Wold theorem that a (Borel) probability measure is characterized by its one-dimensional projections. All of the one-dimensional projections of a multivariate Gaussian are themselves univariate Gaussians, so if we have a distribution all of whose one-dimensional projections are univariate Gaussians it must be a multivariate Gaussian. And indeed, if we consider some random variate $C = \sum_i c_iZ_i$, we find that, writing $x = (X_k)$, $$C = \sum_i c_i \langle x, v_i \rangle = \sum_i c_i \sum_{k=1}^d (v_i)_k X_k.$$ Since the $X_k$ are independent Gaussians, $C$ is itself a Gaussian, and hence $Z_i$ are jointly distributed as a multivariate Gaussian.

We can calculate the mean and covariance of the vector $Z = (Z_i)$ fairly easily. We have $\mathbb{E}[Z_i] = \sum_{j=1}^d (v_i)_j \mathbb{E}[X_j] = 0$ and $$\mathbb{E}[Z_iZ_j] = \sum_{1 \le k, l \le d}(v_i)_k (v_j)_l \mathbb{E}[X_kX_l] = \sum_{k=1}^d (v_i)_k (v_j)_k = \langle v_i, v_j \rangle.$$ In other words, letting $\mathbf{G}$ be the Gram matrix of the $v_i$, we have $$\Pr\left[ \langle x, v_i \rangle \ge 0 | x \sim N(\mathbf{0}, \mathbf{I}_d)\right] = \Pr\left[ Z_i \ge 0 | Z \sim N(\mathbf{0}, \mathbf{G})\right].$$

All this algebra is neat, but the point here is really valuable: we’ve transformed this somewhat-bespoke high-dimensional geometry problem into a more-straightforward (and, thankfully, well-studied) question in probability: what’s the measure of an orthant under a mean-zero Gaussian distribution?

Why should we have expected the Gram matrix to appear?

This can seem pretty magical, so let’s stop and think about why we could have expected the Gram matrix to pop up. The key point is that the Gram matrix is a complete invariant of sets of vectors up to orthogonal transformations. In more invariant-theoretic language, if we consider the simultaneous $O(d)$-action on $(\mathbb{R}^d)^n$, the Gram matrix construction separates $O(d)$-orbits. The spherical volumes we’re interested in are clearly unchanged under orthogonal transformations, so it makes sense to consider whether the geometry of this problem is easily expressed in terms of the Gram matrix. It’s still a little bit surprising to me that the answer is so neat, but hopefully this makes it less mysterious why the Gram matrix pops up.

Gaussian orthants

At this point, we’re left with the problem of computing probabilities of Gaussian orthants. Note that this Gaussian-orthant framing is no more susceptible to naive Monte Carlo integration than our original framing: we still have exponential-in-$n$ decay of the volume fractions. So we have to do a little more work, but thankfully this is a well-studied problem following seminal work of Genz.

Opening moves

Let’s start by generalizing the problem a little bit. Above we consider only mean-zero distributions with domain of integration the positive orthant. We’ll generalize to allow arbitrary means and integration domains arbitrary parallelepipeds, i.e. regions of the form $\mathbf{l} \le \mathbf{Ax} \le \mathbf{u}$ for arbitrary vectors $\mathbf{l}$, $\mathbf{u}$ and an arbitrary matrix $\mathbf{A} \in \mathbb{R}^{n \times n}$ of full rank: that is, integrals of the form $$\int_{\mathbf{l} \le \mathbf{Ax} \le \mathbf{u}} \exp\left[-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu}) \right] d\mathbf{x}.$$

We can immediately convert this more-general case into the special case of mean-zero Gaussians and axis-aligned parallelepipeds (i.e., those with $\mathbf{A} = \mathbf{I}$) via the following two steps:

• A shift $\mathbf{y} = \mathbf{x} - \boldsymbol{\mu}$ centers the Gaussian, and preserves both general and axis-aligned parallelepipeds.

• A linear change of variables $\mathbf{z} = \mathbf{A}\mathbf{y}$ aligns the region of integration along the axes, and changes the covariance matrix
  Recall that the PDF of a multivariate Gaussian uses the precision matrix $\boldsymbol{\Sigma}^{-1}$, not the covariance matrix.
  $\boldsymbol{\Sigma}$ to $\mathbf{A}\boldsymbol{\Sigma}\mathbf{A}^T.$

Genz’s approach

Genz takes as starting point a Gaussian integral over an axis-aligned parallelepiped, i.e. some expression of the form $$P = \int_{\mathbf{a} \le \mathbf{x} \le \mathbf{b}} \exp\left[-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu}) \right] d\mathbf{x}.$$ Having just worked to axis-align our region of integration, the first step in Genz’s approach will be to undo this alignment. Why on earth would we do this? There are two upshots:

• We’ll get to fully standardize the Gaussian integrand.

• Our new region of integration will come with a convenient variable ordering where the limits for $y_i$ depend only on $y_{i-1}, \cdots, y_1$.

The key idea, which is common in multivariate statistics, is simply to take the Cholesky decomposition $\boldsymbol{\Sigma} = \mathbf{L}\mathbf{L}^T$ of the covariance matrix and perform the change of variables $\mathbf{x} = \mathbf{L}\mathbf{y}$ so that the covariance matrix becomes the identity.

The integration constraints $a_i \le (\mathbf{L}\mathbf{y})_i \le b_i$ then reduce to (writing $L_{ij}$ for the entries of $\mathbf{L}$):

$$\frac{a_i - \sum_{j=1}^{i-1} L_{ij} y_j}{L_{ii}} \le y_i \le \frac{b_i - \sum_{j=1}^{i-1} L_{ij} y_j}{L_{ii}}.$$

Note the index constraints on the sums above: because $\mathbf{L}$ is lower-triangular, the integration limits for $y_i$ depend only on the earlier variables $y_1, \ldots, y_{i-1}$. Then, writing $a'_i(y_1, \cdots, y_{i-1})$ and $b'_i(y_1, \cdots, y_{i-1})$ for the left- and right-hand sides of the constraints, our integral becomes $$P = \int_{a'_1}^{b'_1} \frac{1}{\sqrt{2\pi}}\exp(-y_1^2/2) \int_{a'_2(y_1)}^{b'_2(y_1)} \frac{1}{\sqrt{2\pi}}\exp(-y_2^2/2) \cdots \int_{a'_n(y_1, \cdots, y_{n-1})}^{b'_n(y_1, \cdots, y_{n-1})} \frac{1}{\sqrt{2\pi}}\exp(-y_n^2/2) \;d\mathbf{y}.$$

From here, the rest of Genz’s approach is simply some straightforward reparametrizations to make the integral more suitable for off-the-shelf integration methods. We begin by making the substitution $y_i = \Phi^{-1}(z_i)$, where $\Phi$ is the standard normal CDF. Writing $d_i = \Phi(a'_i)$, $e_i = \Phi(b'_i)$, our integral becomes

$$P = \int_{d_1}^{e_1} \int_{d_2(z_1)}^{e_2(z_1)} \cdots \int_{d_n(z_1,\ldots,z_{n-1})}^{e_n(z_1,\ldots,z_{n-1})} dz_n \cdots dz_1.$$

Now the integrand is really simple–it’s literally 1–but the region of integration is complicated. One more reparametrization fixes that: writing $z_i = d_i + w_i(e_i-d_i)$ allows us to integrate in the $w_i$ over the unit cube as

$$\int_{[0,1]^n} (e_1 - d_1)(e_2 - d_2) \cdots (e_n - d_n) \; d\mathbf{w}.$$

Carrying all the substitutions through (though leaving the $y_i$ as intermediate variables), we have

$$y_{i-1} = \Phi^{-1}(d_{i-1} + w_{i-1}(e_{i-1} - d_{i-1})),$$ $$d_i = \Phi\!\left(\frac{a_i - \sum_{j=1}^{i-1} L_{ij} y_j}{L_{ii}}\right), e_i = \Phi\!\left(\frac{b_i - \sum_{j=1}^{i-1} L_{ij} y_j}{L_{ii}}\right).$$

If we squint, each $e_i-d_i$ factor in the integrand is a conditional probability: explicitly, it’s the probability that a standard normal falls in $[a'_i, b'_i]$ given the values of $y_1, \ldots, y_{i-1}$, so that the full product is a telescoping product of conditional probabilities coming from the Cholesky decomposition and the integral computes the joint probability we wanted.

Wrapping up

At this point we can just turn to standard Monte Carlo techniques and get a pretty good estimator in moderate dimensions. I believe that the sequential nature of the integrand (i.e., as a product of conditional probabilities) lends itself well to more-complicated Monte Carlo algorithms that save us relatively-expensive evaluations of Gaussian CDFs and inverse CDFs, but I’m not familiar enough with this literature to say anything more specific.

There’s also a subsidiary question of what variable ordering we should choose; this is essentially the only input to our Cholesky decomposition. Genz suggests having the variables with the widest integration intervals $[a_i, b_i]$ innermost, and notes that this can make a substantial speedup; I think this question is fleshed out a bit more in follow-up literature.

To actually compute these orthant probabilities, we can use the R package mvtnorm written by Genz himself. Analogous implementations exist in SciPy, MATLAB, and Julia, and possibly elsewhere.