O(N) Hierarchical algorithm for computing the expectations of truncated multi-variate normal distributions in N dimensions
Jingfang Huang, Fuhui Fang, George Turkiyyah, Jian Cao, Marc G. Genton, David E. Keyes
TL;DR
This work tackles the high-dimensional TMVN expectation $\phi(\mathbf{a}, \mathbf{b}; A)$ and develops an asymptotically optimal $O(N)$ algorithm exploiting hierarchical low-rank structure in $A$ and a low-rank form for $H(\mathbf{x})$.A divide-and-conquer strategy on a hierarchical tree decouples variables through Fourier transforms, performing a downward pass to relate parent and child effective variables and an upward pass to assemble Fourier series representations of low-dimensional functions at each node.Two concrete matrix models are analyzed: (i) SPD tridiagonal systems and (ii) exponential covariance-based matrices, with potential-theory justifications via Green's functions and dimension-reduction techniques to maintain at most two effective variables per node.Numerical results demonstrate accuracy and near-linear scaling with dimension for both cases, while outlining generalizations and limitations, including the need for careful basis choices and handling of large condition numbers.
Abstract
In this paper, we study the $N$-dimensional integral $φ(a,b; A) = \int_{a}^{b} H(x) f(x | A) \text{d} x$ representing the expectation of a function $H(X)$ where $f(x | A)$ is the truncated multi-variate normal (TMVN) distribution with zero mean, $x$ is the vector of integration variables for the $N$-dimensional random vector $X$, $A$ is the inverse of the covariance matrix $Σ$, and $a$ and $b$ are constant vectors. We present a new hierarchical algorithm which can evaluate $φ(a,b; A)$ using asymptotically optimal $O(N)$ operations when $A$ has "low-rank" blocks with "low-dimensional" features and $H(x)$ is "low-rank". We demonstrate the divide-and-conquer idea when $A$ is a symmetric positive definite tridiagonal matrix, and present the necessary building blocks and rigorous potential theory based algorithm analysis when $A$ is given by the exponential covariance model. Numerical results are presented to demonstrate the algorithm accuracy and efficiency for these two cases. We also briefly discuss how the algorithm can be generalized to a wider class of covariance models and its limitations.
