Orthant Normal Distribution

• Orthant normal distribution is defined as the probability that all components of a multivariate Gaussian are nonnegative, with applications in Gaussian process theory and polyhedral geometry.
• Specialized techniques like Steck’s reduction for equicorrelated structures and asymptotic bounds provide efficient estimation methods for these multidimensional probabilities.
• The holonomic gradient method reformulates the problem as a system of ODEs, offering enhanced numerical stability and competitive accuracy for high-dimensional integration challenges.

An orthant normal distribution is a classical but highly nontrivial object in the study of multivariate Gaussian measures: it refers to the computation of the probability that a multivariate normal vector falls within a specified orthant, most often the non-negative orthant, i.e., the set $\{x\in\mathbb{R}^n:x_1\ge0,\dots,x_n\ge0\}$. The study of such probabilities is essential in multivariate probability, Gaussian process theory, polyhedral geometry, and has nontrivial applications in random polynomials on the simplex. For equicorrelated covariance structures, this problem admits specialized reductions and asymptotic analysis, while modern computation exploits connections to holonomic systems.

1. Formal Definition and Integral Representation

Let $X=(X_1,\dots,X_n)\sim N_n(0,\Sigma)$, where every $X_i$ is real-valued and $\Sigma$ is a symmetric positive-definite covariance matrix. The orthant probability for the non-negative orthant is

$$f(n,\Sigma) = P\{X_1\ge0,\ldots,X_n\ge0\} = \int_{[0,\infty)^n} \phi_n(x;0,\Sigma)\,dx$$

with multivariate density

$$\phi_n(x;0,\Sigma) = (2\pi)^{-n/2}\det(\Sigma)^{-1/2}\exp\left(-\frac{1}{2}x^T\Sigma^{-1}x\right).$$

For the equicorrelated case $\Sigma_{ij} = 1$ for $i = j$ and $\Sigma_{ij} = \rho\in(0,1)$ for $i\ne j$, the orthant probability admits a one-dimensional reduction due to Steck (1962): $f(n,\rho) = \mathbb{E}[\Phi^n(Z\sqrt{s})]$ where $s=\rho/(1-\rho)$, $Z\sim N(0,1)$, and $\Phi$ is the univariate normal CDF (Pinasco et al., 2020). More generally,

$$\Phi_n(\mu,\Sigma)=P(X_1\ge0,\dots,X_n\ge0) = \int_{[0,\infty)^n} \phi_n(x;\mu,\Sigma)\,dx$$

with extensions allowing arbitrary mean and covariance (Koyama et al., 2012).

2. Analytical Bounds for the Equicorrelated Orthant Probability

For the equicorrelated normal, $f(n, \rho)$ displays distinct asymptotic regimes depending on $\rho$:

• For $\rho > 1/2$ and $n\ge2$: There exists a universal constant $c_1>0$ such that

$$\frac{c_1}{\sqrt{2+(1/\rho-1)\log n}} \le \frac{f(n,\rho)}{n^{1-1/\rho}} \le \frac{1}{2^{2-1/\rho}\sqrt{\frac{1-\rho}{\rho}(1+B(2,1/\rho-1))}}$$

with $B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b)$ the Beta function.

• For $\rho < 1/2$ and suitably large $n$:

$$\frac{1}{2^{2-1/\rho}\sqrt{\frac{1-\rho}{\rho}(\Gamma(1/\rho-1)-1)}} \le \frac{f(n,\rho)}{n^{1-1/\rho}}$$

and for all $n\ge n_0(\rho)$,

$$\frac{f(n,\rho)}{n^{1-1/\rho}} \le \sqrt{\frac{1-\rho}{\rho}[(1/\rho-1)\log^2 n]^{1/\rho-2}}$$

For all fixed $\rho\in(0,1)$, $f(n,\rho) \asymp n^{1-1/\rho}$ up to polylogarithmic factors in $n$ (Pinasco et al., 2020).

Steck’s reduction, Markov inequalities, and Mill’s ratio-type bounds underpin these estimates, with critical refinements available for the logarithmic correction regimes in $n$.

3. Holonomic Gradient Method and Pfaffian Systems

Direct computation of orthant probabilities for arbitrary covariance matrices is numerically challenging. The holonomic gradient method (HGM) establishes a system of linear PDEs—the holonomic system—annihilating the multivariate integral

$$g(x,y) = \int_{[0,\infty)^n} \exp\left(\frac{1}{2}\sum_{i,j=1}^n x_{ij} t_i t_j + \sum_{i=1}^n y_i t_i\right) dt_1\cdots dt_n,$$

where parameters are $x=-\Sigma^{-1}$ and $y=\Sigma^{-1}\mu$ (Koyama et al., 2012). The corresponding orthant probability is then

$$\Phi_n(\mu,\Sigma) = (2\pi)^{-n/2}|\Sigma|^{-1/2}\exp\left(\frac{1}{2}\mu^T\Sigma^{-1}\mu\right)g(x,y).$$

The HGM reformulates the problem as a first-order Pfaffian ODE for the vector $G(x,y) = \{g_J(x,y): J\subset\{1,\ldots,n\}\}$ of $2^n$ components, with recurrence relations for derivatives in $x_{ij}$ and $y_i$. This reduction ensures that the numerical procedure is well-conditioned as long as the covariance remains positive definite.

The algorithm iteratively integrates from a diagonal start (zero mean, diagonal covariance) to the desired $(x_{\mathrm{target}}, y_{\mathrm{target}})$ along a straight path, solving the ODE with a method such as fourth-order Runge–Kutta. The principal computational cost stems from the exponential state scaling, $O(n^2 2^n N_\text{steps})$, but practical calculations are reported up to $n=50$ under moderate accuracy (Koyama et al., 2012).

4. Plackett Recurrence and Comparative Numerical Strategies

Classical approaches, such as Plackett's 1954 recurrence, relate derivatives of orthant probabilities to lower-dimensional analogs through explicit recurrence relations. However, the recursion involves denominators $\sqrt{1-\rho_{ij}^2}$, introducing numerical instability for near-singular correlation matrices or as $|\rho_{ij}| \uparrow 1$.

The holonomic gradient method, by contrast, eliminates such problematic factors and operates on rational recurrences over an exponentially-sized state space, yielding improved stability on high-dimensional or ill-conditioned problems. Comparative studies indicate HGM’s competitive or superior computational timing for large $n$ and improved accuracy benchmarks up to $10^{-9}$ (Koyama et al., 2012).

5. Asymptotics and Applications in Random Polynomials on the Simplex

A substantial application concerns random homogeneous polynomials evaluated on the standard $n$-simplex $\Delta$. Specifically, for a random Bombieri-normalized $k$-homogeneous polynomial $P$, the event that $P$ attains a (relative) maximum at a vertex $a_i$ is equivalent to the event that an associated equicorrelated multivariate normal $X(a_i)$ has all positive entries. The relevant correlation parameter for vertices is

$$\rho_n = \frac{nk + (k-1)}{n(k+1) + (k-1)} \to \frac{k}{k+1} \quad \text{as } n\to\infty.$$

For $k>4$, the events $\{P \text{ has max at } a_i\}$ become almost independent in total variation as $n\to\infty$, and the probability that $P$ has a maximum at some vertex converges to $1$ (Pinasco et al., 2020).

6. Numerical Implementation and Performance

The HGM algorithm operates through an initial parameter transformation to $(x, y)$, initialization on a diagonal covariance, ODE setup via block-sparse recurrences, and numerical integration along a chosen path. Key implementation items include:

• Sparse matrix optimization given that each recurrence couples subsets differing by at most two elements.
• Precomputation of recurrence sparsity patterns for efficiency.
• Monitoring of covariance positivity to stay within the holonomic domain.

Numerical experiments reported by Koyama & Takemura show high-precision agreement with exact results (up to $10^{-9}$), competitive runtimes for small $n$ compared to Miwa–Hayter–Kuriki recursion, and clear acceleration for large $n$ (e.g., HGM is faster for $n=20,50$). Limitations can occur when $\|\mu\|$ is large and $\Sigma$ nearly singular, causing ODE stiffness (Koyama et al., 2012).

7. Connections to Broader Contexts

Orthant normal probabilities underpin various problems in statistical inference, Gaussian process excursion sets, and the study of random fields, as well as in polyhedral geometry and stochastic optimization. Their computation, both by analytic reduction in symmetric (equicorrelated) cases and by holonomic systems in the general setting, remains an active area with substantial algorithmic and theoretical developments (Pinasco et al., 2020, Koyama et al., 2012).