Fast Simulation of Hyperplane-Truncated Multivariate Normal Distributions

1 Introduction

We investigate the problem of simulation from a multivariate normal (MVN) distribution whose samples are restricted to the intersection of a set of hyperplanes, which is shown to be inherently related to the simulation of a conditional distribution of a MVN distribution. A naive approach, which linearly transforms a random variable drawn from the conditional distribution of a related MVN distribution, requires a large number of intermediate variables that are often computationally expensive to instantiate. To address this issue, we propose a fast and exact simulation algorithm that directly projects a MVN random variable onto the intersection of a set of hyperplanes. We further show that sampling from a MVN distribution, whose covariance (precision) matrix can be decomposed as the sum (difference) of a positive-definite matrix, whose inversion is known or easy to compute, and a low-rank symmetric matrix, may also be made significantly fast by exploiting this newly proposed stimulation algorithm for hyperplane-truncated MVN distributions, avoiding the need of Cholesky decomposition that has a computational complexity of $O(k^{3})$ (Golub and Van Loan, 2012), where $k$ is the dimension of the MVN random variable.

2 Hyperplane-truncated and conditional MVNs

For the problem under study, we express a $k$-dimensional MVN distribution truncated on the intersection of $k_{2}<k$ hyperplanes as

and $\mbox{Rank}({\bf G})=k_{2}$. The probability density function can be expressed as

where $Z$ is a constant ensuring $\int p(\boldsymbol{x}\,|\,\boldsymbol{\mu},\boldsymbol{\Sigma},{\bf G},\boldsymbol{r})d\boldsymbol{x}=1$, and $\delta(x)=1$ if the condition $x$ is satisfied and $\delta(x)=0$ otherwise. Let us partition ${\bf G}$, $\boldsymbol{x}$, $\boldsymbol{\mu}$, $\boldsymbol{\Sigma}$, and $\boldsymbol{\Lambda}=\boldsymbol{\Sigma}^{-1}$ as

respectively, where $k=k_{1}+k_{2}$, $\boldsymbol{\Sigma}_{21}=\boldsymbol{\Sigma}_{12}^{T}$, and $\boldsymbol{\Lambda}_{21}=\boldsymbol{\Lambda}_{12}^{T}$.

A special case that frequently arises in real applications is when ${\bf G}_{1}=\mathbf{0}_{k_{2}\times k_{1}}$ and ${\bf G}_{2}={\bf I}_{k_{2}}$, which means $(\mathbf{0}_{k_{2}\times k_{1}},{\bf I}_{k_{2}})\boldsymbol{x}=\boldsymbol{x}_{2}=\boldsymbol{r}$ and the need is to simulate $\boldsymbol{x}_{1}$ given $\boldsymbol{x}_{2}=\boldsymbol{r}$. For a MVN random variable $\boldsymbol{x}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, it is well known, $e.g.$, in Tong (2012), that the conditional distribution of $\boldsymbol{x}_{1}$ given $\boldsymbol{x}_{2}=\boldsymbol{r}$, $i.e.$, the distribution of $\boldsymbol{x}$ restricted to $\mathcal{S}=\{\boldsymbol{x}:(\mathbf{0}_{k_{2}\times k_{1}},{\bf I}_{k_{2}})\boldsymbol{x}=\boldsymbol{r}\}$, can be expressed as

Alternatively, applying the Woodbury matrix identity to relate the entries of the covariance matrix $\boldsymbol{\Sigma}$ to those of the precision matrix $\boldsymbol{\Lambda}$, one may obtain the following equivalent expression as

More specifically, denoting $\boldsymbol{\Lambda}=[{\bf H}^{-1}\boldsymbol{\Sigma}({\bf H}^{-1})^{T}]^{-1}={\bf H}^{T}\boldsymbol{\Sigma}^{-1}{\bf H}$ as the precision matrix for $\boldsymbol{z}$, we have

and hence $\boldsymbol{x}$ truncated on $\mathcal{S}$ can be naively generated using the following algorithm, whose computational complexity is described in Table 1 of the Appendix.

For illustration, we consider a simple 2-dimensional example with $\boldsymbol{\mu}=(1,1.2)^{T}$, $\boldsymbol{\Sigma}=[(1,0.3)^{T},(0.3,1)^{T}]$, ${\bf G}=(1,1)$, and $\boldsymbol{r}=1$. If we choose ${\bf H}_{1}=(-0.7071,0.7071)^{T}$ and ${\bf H}_{2}=(1.3,1.3)^{T}$, then we have $z_{2}=({\bf G}{\bf H}_{2})^{-1}\boldsymbol{r}=(2.6)^{-1}=0.3846$, $\boldsymbol{\Lambda}_{11}=1.4285$, and $\boldsymbol{\Lambda}_{12}=0$; as shown in Figure 1, we may generate $\boldsymbol{x}$ using

where $z_{1}\sim\mathcal{N}(0.1414,0.7)$ and $z_{2}=0.3846$.

For high dimensional problems, however, Algorithm 1 in general requires a large number of intermediate variables that could be computationally expensive to compute. In the following discussion, we will show how to completely avoid instantiating these intermediate variables.

3 Fast and exact simulation of MVN distributions

Instead of using Algorithm 1, we first provide a theorem to show how to efficiently and exactly simulate from a hyperplane-truncated MVN distribution. In the Appendix, we provide two different proofs. The first proof facilitates the derivations by employing an existing algorithm of Hoffman and Ribak (1991) and Doucet (2010), which describes how to simulate from the conditional distribution of a MVN distribution shown in (4) without computing $\boldsymbol{\Sigma}_{11}-\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\boldsymbol{\Sigma}_{21}$ and its Cholesky decomposition. Note it is straightforward to verify that the algorithm in Hoffman and Ribak (1991) and Doucet (2010), as shown in the Appendix, can be considered as a special case of the proposed algorithm with ${\bf G}=[\boldsymbol{0},{\bf I}]$.

The above algorithm and theorem, whose computational complexity is described in Table 2 of the Appendix, show that one may draw $\boldsymbol{y}$ from the unconstrained MVN as $\boldsymbol{y}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ and directly map it to a vector $\boldsymbol{x}$ on the intersection of hyperplanes using $\boldsymbol{x}=\boldsymbol{\Sigma}{\bf G}^{T}({\bf G}\boldsymbol{\Sigma}{\bf G}^{T})^{-1}\boldsymbol{r}+\left[{\bf I}-\boldsymbol{\Sigma}{\bf G}^{T}({\bf G}\boldsymbol{\Sigma}{\bf G}^{T})^{-1}{\bf G}\right]\boldsymbol{y}.$ For illustration, with the same $\boldsymbol{\mu}$, $\boldsymbol{\Sigma}$, ${\bf G}$, and $\boldsymbol{r}$ as those in Figure 1, we show in Figure 2 a simple two dimensional example, where the unrestricted Gaussian distribution $\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ is represented with a set of ellipses, and the constrained sample space $\mathcal{S}$ is represented as a straight line in the two-dimensional setting. With $\boldsymbol{\Sigma}{\bf G}^{T}({\bf G}\boldsymbol{\Sigma}{\bf G}^{T})^{-1}\boldsymbol{r}=(0.5,0.5)^{T}$, $\left[{\bf I}-\boldsymbol{\Sigma}{\bf G}^{T}({\bf G}\boldsymbol{\Sigma}{\bf G}^{T})^{-1}{\bf G}\right]=[(0.5,-0.5)^{T},(-0.5,0.5)^{T}]$, one may directly maps a sample $\boldsymbol{y}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ to a vector on the constrained space. For example, if $\boldsymbol{y}=(1,2)^{T}$, then it would be mapped to $\boldsymbol{x}=(0,1)^{T}$ on the straight line.

For fast simulation of MVN distributions with structured covariance or precision matrices, our idea is to relate them to higher-dimensional hyperplane-truncated MVN distributions, with block-diagonal covariance matrices, that can be efficiently simulated with Algorithm 2. We first introduce an efficient algorithm for the simulation of a MVN distribution, whose covariance matrix is a positive-definite matrix subtracted by a low-rank symmetric matrix. Such kind of covariance matrices commonly arise in the conditional distributions of MVN distributions, as shown in (4). We then further extend this algorithm to the simulation of a MVN distribution whose precision (inverse covariance) matrix is the sum of a positive-definite matrix and a low-rank symmetric matrix. Such kind of precision matrices commonly arise in the conditional posterior distributions of the regression coefficients in both linear regression and generalized linear models.

The probability density function (PDF) of the MVN distribution

is the same as the PDF of the marginal distribution of $\boldsymbol{x}_{1}=(x_{1},\ldots,x_{k_{1}})^{T}$ in $\boldsymbol{x}=(\boldsymbol{x}_{1}^{T},x_{k_{1}+1},\ldots,x_{k})^{T}$, whose PDF is expressed as

is the same as the PDF of the marginal distribution of $\boldsymbol{x}_{1}$ in $\boldsymbol{x}=(\boldsymbol{x}_{1}^{T},x_{k_{1}+1},\ldots,x_{k})^{T}$, whose PDF is expressed as

where $\boldsymbol{x}_{2}=(x_{k_{1}+1},\ldots,x_{k})^{T}$, $Z$ is a normalization constant, and

Further applying Theorem 1 to Corollary 3, as described in detail in the Appendix, a MVN random variable $\boldsymbol{x}$ with a structured covariance matrix can be generated as in Algorithm 3, where there is no need to compute $\boldsymbol{\Sigma}_{11}-\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\boldsymbol{\Sigma}_{21}$ and its Cholesky decomposition. Suppose the covariance matrix $\boldsymbol{\Sigma}_{11}$ admits some special structure that makes it easy to invert and computationally efficient to simulate from $\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}_{11})$, then Algorithm 3 could lead to a significant saving in computation if $k_{2}\ll k_{1}$. On the other hand, when $k_{2}\gg k_{1}$ and $\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}$ admits no special structures, Algorithm 3 may not bring any computational advantage and hence one may resort to the naive Cholesky decomposition based procedure. Detailed computational complexity analyses for both methods are provided in Tables 3 and 4 of the Appendix, respectively.

A random variable simulated with Algorithm 3 is distributed as $\boldsymbol{x}_{1}\sim\mathcal{N}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11}-\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\boldsymbol{\Sigma}_{21}).$

The random variable obtained with Algorithm 4 is distributed as $\boldsymbol{\beta}\sim\mathcal{N}(\boldsymbol{\mu}_{\beta},\boldsymbol{\Sigma}_{\beta})$, where $\boldsymbol{\Sigma}_{\beta}=({\bf A}+\boldsymbol{\Phi}^{T}\boldsymbol{\Omega}\boldsymbol{\Phi})^{-1}$.

4 Illustrations

Below we provide several examples to illustrate Theorem 1, which shows how to efficiently simulate from a hyperplane-truncated MVN distribution, and Corollary 4 (Corollary 5), which shows how to efficiently simulate from a MVN distribution with a structured covariance (precision) matrix. We run all our experiments on a 2.9 GHz computer.

We first compare Algorithms 1 and 2, whose generated random samples follow the same distribution, as suggested by Theorem 1, to highlight the advantages of Algorithm 2 over Algorithm 1. We then employ Algorithm 2 for a real application whose data dimension is high and sample size is large.

First, to verify Theorem 1, we conduct an experiment with $k=5000$ data dimension, $k_{2}=20$ hyperplanes, and a diagonal $\boldsymbol{\Sigma}$. Contour plots of two randomly selected dimensions of the 10,000 random samples simulated with Algorithms 1 and 2 are shown in the top and bottom rows of Figure 3, respectively. The clear matches between the contour plots of these two different algorithms suggest the correctness of Theorem 1.

To demonstrate the efficiency of Algorithm 2, we first carry out a series of experiments with the number of hyperplane constraints fixed at $k_{2}=20$ and the data dimension increased from $k=50$ to $k=5000$. The computation time of simulating 10,000 samples averaged over five random trials is shown in Figure 4 for non-diagonal $\boldsymbol{\Sigma}$’s and in Figure 4 for diagonal ones. It is clear that, when the data dimension $k$ is high, Algorithm 2 has a clear advantage over Algorithm 1 by avoiding computing unnecessary intermediate variables, which is especially evident when $\boldsymbol{\Sigma}$ is diagonal. We then carry out a series of experiments where we vary not only $k$, but also $k_{2}$ from $0.1k$ to $0.9k$ for each $k$. As shown in Figure 4, it is evident that Algorithm 2 dominates Algorithm 1 in all scenarios, which can be explained by the fact that Algorithm 2 needs to compute much fewer intermediate variables. Also observed is that a larger $k_{2}$ leads to slower simulation for both algorithms, but to a much lesser extent for Algorithm 2. Moreover, the curvatures of those curves indicate that Algorithm 2 is more practical in a high dimensional setting. Note that since Algorithm 2 can naturally exploit the structure of the covariance matrix $\boldsymbol{\Sigma}$ for fast simulation, it is clearly more capable of benefiting from having a diagonal or block-diagonal $\boldsymbol{\Sigma}$, demonstrated by comparing Figures 4 and 4 with Figures 4 and 4. All these observations agree with our computational complexity analyses for Algorithms 1 and 2, as shown in Table 1 and 2 of the Appendix, respectively.

A practical application of Algorithm 2

More specifically, we focus on a big data setting in which the globally shared simplex-constrained model parameters could be linked to some latent counts via the multinomial likelihood. When there are tens of thousands or millions of observations in the dataset, scalable Bayesian inference for the simplex-constrained globally shared model parameters is highly desired, for example, for inferring the topics’ distributions over words in latent Dirichlet allocation (Blei et al., 2003; Hoffman et al., 2010) and Poisson factor analysis (Zhou et al., 2012, 2016).

However, this batch-learning inference procedure would become inefficient and even impractical when the dataset size $N$ grows to a level that makes it too time consuming to finish even a single iteration of updating all local variables $n_{vj\kappa}$. To address this issue, we consider constructing a mini-batch based Bayesian inference procedure that could make substantial progress in posterior simulation while the batch-learning one may still be waiting to finish a single iteration.

Instead of waiting for all $\boldsymbol{n}_{j}$ to be updated before performing a single update of $\boldsymbol{\phi}$, we develop a mini-batch based Bayesian inference algorithm under a general framework for constructing stochastic gradient Markov chain Monte Carlo (SG-MCMC) (Ma et al., 2015), allowing $\boldsymbol{\phi}$ to be updated every time a mini-batch of $\boldsymbol{n}_{j}$ are processed. For the sake of completeness, we concisely describe the derivation for a SG-MCMC algorithm, as outlined below, for simplex-constrained globally shared model parameters. We refer the readers to Cong et al. (2017) for more details on the derivation and its application to scalable inference for topic modeling.

It is clear that (16) corresponds to simulation of a $V-1$ dimensional truncated MVN distribution with $V$ inequality constraints. Since the number of constraints is larger than the dimension, previously proposed iterative simulation methods such as the one in Botev (2016) are often inappropriate. Note that, by omitting the non-negative constraints, the update in (17) corresponds to simulation of a hyperplane-truncated MVN simulation with a diagonal covariance matrix, which can be efficiently sampled as described in the following example.

Example 1: Simulation of a hyperplane-truncated MVN distribution as

Return $\boldsymbol{x}=\boldsymbol{y}+(1-\mathbf{1}^{T}\boldsymbol{y})\boldsymbol{\phi}$.

Example 2: Simulation from (16) can be approximately but rapidly realized as

Calculate $\boldsymbol{z}=\boldsymbol{y}+(1-\mathbf{1}^{T}\boldsymbol{y})\boldsymbol{\phi}_{t}$;

For comparison, we choose the same SG-MCMC inference procedure but consider simulating (16), as performed every time a mini-batch of data samples are provided, either as in Example 2 or with the Gibbs sampler of Rodriguez-Yam et al. (2004). Simulating (16) with the Gibbs sampler of Rodriguez-Yam et al. (2004) is realized by updating all the $V$ dimensions, one dimension at a time, in each Gibbs sampling iteration. We set the total number of Gibbs sampling iterations for (16) in each mini-batch based update as 1, 5, or 10. Note that in practice, the $\boldsymbol{n}_{j}$ belonging to the current mini-batch are often latent and are updated conditioning on the data samples in the mini-batch and $\boldsymbol{\phi}$. For simplicity, all $\boldsymbol{n}_{j}$ here are simulated once and then fixed.

Using $\boldsymbol{\phi}_{post}^{*}=(\sum_{j=1}^{N}\boldsymbol{n}_{j}+\eta)/(\sum_{j=1}^{N}n_{\boldsymbol{\cdot}j}+\eta V)$, the posterior mean of $\boldsymbol{\phi}$ in a batch-learning setting, as the reference, we show in Figure 5 how the residual errors for the estimated $\boldsymbol{\phi}^{*}$, defined as $\left\|\boldsymbol{\phi}^{*}-\boldsymbol{\phi}\right\|_{2}$, change both as a function of the number of processed mini-batches and as a function of computation time under various settings of the mini-batch based SG-MCMC algorithm. The curves shown in Figure 5 suggest that for each mini-batch, to simulate (16) with the Gibbs sampler of Rodriguez-Yam et al. (2004), it is necessary to have more than one Gibbs sampling iteration to achieve satisfactory results. It is clear from Figure 5 that the Gibbs sampler with 5 or 10 iterations for each mini-batch, even though each mini-batch has only 10 data samples, provides residual errors that quickly approach that of the batch posterior mean with a tiny gap, indicating the effectiveness of the SG-MCMC updating in (16). While simulating (16) with Gibbs sampling could in theory lead to unbiased samples if the number of Gibbs sampling iterations is large enough, it is much more efficient to simulate (16) with the procedure described in Example 2, which provides a performance that is undistinguishable from those of the Gibbs sampler with as many as 5 or 10 iterations for each mini-batch, but at the expense of a tiny fraction of a single Gibbs sampling iteration.

4.2 Simulation of MVNs with structured covariance matrices

To illustrate Corollary 4, we mimic the truncated MVN simulation in (16) and present the following simulation example with a structured covariance matrix.

Example 3: Simulation of a MVN distribution as

Return $\boldsymbol{x}_{1}=\boldsymbol{\mu}_{1}+\boldsymbol{y}_{1}-(\mathbf{1}^{T}\boldsymbol{y}_{1}+a\boldsymbol{y}_{2})\boldsymbol{\phi}_{1}$.

Denoting $\boldsymbol{x}=(\boldsymbol{x}_{1}^{T},x_{k})^{T}$, $\boldsymbol{\phi}=(\boldsymbol{\phi}_{1}^{T},\phi_{k})^{T}$, $\boldsymbol{\mu}=(\boldsymbol{\mu}_{1}^{T},\mu_{k})^{T}$, and $\mu_{k}=1-\mathbf{1}^{T}\boldsymbol{\mu}_{1}$, the above sampling steps can also be equivalently expressed as follows.

Return $\boldsymbol{x}_{1}=\boldsymbol{y}_{1}+(1-\mathbf{1}^{T}\boldsymbol{y})\boldsymbol{\phi}_{1}$.

As shown in Figure 6, for the proposed Algorithm 3, the average time of simulating 10,000 random samples increases linearly in the dimension $k$. By contrast, for the naive Cholesky decomposition based simulation algorithm, whose computational complexity is $O(k^{3})$ (Golub and Van Loan, 2012), the average simulation time increases at a significantly faster rate as the dimension $k$ increases.

For explicit verification, with the 10,000 simulated $k=10^{4}$ dimensional random samples in a random trial, we randomly choose two dimensions and display their joint distribution using a contour plot. As in Figure 7, shown in the first row are the contour plots of five different random trials for the naive Cholesky implementation, whereas shown in the second row are the corresponding ones for the proposed Algorithm 3. As expected, the contour lines of the two figures in the same column closely match each other.

To further examine when to apply Algorithm 3 instead of the naive Cholesky decomposition based implementation in a general setting, we present the computational complexity analyses in Tables 3 and 4 of the Appendix for the naive approach and Algorithm 3, respectively. In addition, we mimic the settings in Section 0.4.1 to conduct a set of experiments with randomly generated $\boldsymbol{\Sigma}_{12}$, diagonal $\boldsymbol{\Sigma}_{11}$, and diagonal $\boldsymbol{\Sigma}_{22}$. We fix $k_{1}=4000$ and vary $k_{2}$ from 1 to 8000. The computation time for one sample averaged over 50 random trials is presented in Figure 8. It is clear from Tables 3 and 4 and Figure 8 that, as a general guideline, one may choose Algorithm 3 when $k_{2}$ is smaller than $k_{1}$ and $\boldsymbol{\Sigma}_{11}$ admits some special structure that makes it easy to invert and computationally efficient to simulate from $\mathcal{N}(\bf{0},\boldsymbol{\Sigma}_{11})$.

4.3 Simulation of MVNs with structured precision matrices

To examine when to apply Algorithm 4 instead of the naive Choleskey decomposition based procedure, we first consider a series of random simulations in which the sample size $n$ is fixed while the data dimension $p$ is varying. We then show that Algorithm 4 can be applied for high-dimensional regression whose $p$ is often much larger than $n$.

We fix $n=4000$, vary $p$ from 1 to 8000, and mimic the settings in Section 0.4.1 to randomly generate $\boldsymbol{\Phi}$, diagonal ${\bf A}$, and diagonal $\boldsymbol{\Omega}$. As a function of dimensions $p$, the computation time for one sample averaged over 50 random trials is shown in Figure 9. It is evident that, identical to the complexity analysis in Tables 5 and 6, Algorithm 4 has a linear complexity with respect to $p$ under these settings, which will bring significant acceleration in a high-dimensional setting with $p\gg n$. If the sample size $n$ is large enough that $n>p$, then one may directly apply the naive Cholesky decomposition based implementation.

Example 4: Simulation of the MVN distribution

Sample $\boldsymbol{y}_{1}\sim\mathcal{N}(\mathbf{0},{\bf A}^{-1})$ and $\boldsymbol{y}_{2}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Omega}^{-1})$ ;

Return $\boldsymbol{\beta}=\boldsymbol{y}_{1}+{\bf A}^{-1}\boldsymbol{\Phi}^{T}(\boldsymbol{\Omega}^{-1}+\boldsymbol{\Phi}{\bf A}^{-1}\boldsymbol{\Phi}^{T})^{-1}\left(\boldsymbol{t}-\boldsymbol{\Phi}\boldsymbol{y}_{1}-\boldsymbol{y}_{2}\right)$, which can be realized using

Solve $\boldsymbol{\alpha}$ such that $(\boldsymbol{\Omega}^{-1}+\boldsymbol{\Phi}{\bf A}^{-1}\boldsymbol{\Phi}^{T})\boldsymbol{\alpha}=\boldsymbol{t}-\boldsymbol{\Phi}\boldsymbol{y}_{1}-\boldsymbol{y}_{2}$;

Return $\boldsymbol{\beta}=\boldsymbol{y}_{1}+{\bf A}^{-1}\boldsymbol{\Phi}^{T}\boldsymbol{\alpha}$.

Note that if $\boldsymbol{\Omega}={\bf I}_{n}$, then the simulation algorithm in Example 4 reduces to the one in Proposition 2.1 of Bhattacharya et al. (2016), which is shown there to be significantly more efficient than that of Rue (2001) for high-dimensional regression if $p\gg n$.

5 Conclusions

A fast and exact simulation algorithm is developed for a multivariate normal (MVN) distribution whose sample space is constrained on the intersection of a set of hyperplanes, which is shown to be inherently related to the conditional distribution of a unconstrained MVN distribution. The proposed simulation algorithm is further generalized to efficiently simulate from a MVN distribution, whose covariance (precision) matrix can be decomposed as the sum (difference) of a positive-definite matrix and a low-rank symmetric matrix, using a higher dimensional hyperplane-truncated MVN distribution whose covariance matrix is block-diagonal.

References

Proofs

Let us denote $\boldsymbol{\Lambda}={\bf H}^{T}\boldsymbol{\Sigma}^{-1}{\bf H}$ as a precision matrix that can be partitioned as in (7). For Algorithm 1, instead of directly simulating $\boldsymbol{z}_{1}$ given $\boldsymbol{z}_{2}$ using the conditional distribution of the MVN, we apply Algorithm 5 (Hoffman and Ribak, 1991; Doucet, 2010) to modify its sampling steps as follows.

Therefore, we can equivalently generate $\boldsymbol{x}$ as follows.

Sample $\boldsymbol{y}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$.

Return $\boldsymbol{x}=({\bf H}_{1},{\bf H}_{1}\boldsymbol{\Lambda}_{11}^{-1}\boldsymbol{\Lambda}_{12}){\bf H}^{-1}\boldsymbol{y}+({\bf H}_{2}-{\bf H}_{1}\boldsymbol{\Lambda}_{11}^{-1}\boldsymbol{\Lambda}_{12})({\bf G}{\bf H}_{2})^{-1}\boldsymbol{r}$.

The computation can be significantly simplified if $\boldsymbol{\Lambda}_{12}=\mathbf{0}$, which means

Since ${\bf H}_{1}^{T}{\bf G}^{T}=0$ by definition, to make $\boldsymbol{\Lambda}_{12}=\mathbf{0}$, if and only if we have ${\bf H}_{2}$ as

Sample $\boldsymbol{y}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$.

Return $\boldsymbol{x}=({\bf H}_{1},\mathbf{0}_{k\times k_{2}}){\bf H}^{-1}\boldsymbol{y}+\boldsymbol{\Sigma}{\bf G}^{T}({\bf G}\boldsymbol{\Sigma}{\bf G}^{T})^{-1}\boldsymbol{r}$, or return

Let us denote $(\mathbf{0}_{k\times k_{1}},\boldsymbol{\Sigma}{\bf G}^{T}{\bf M}){\bf H}^{-1}={\bf C}$. We have

and hence ${\bf C}{\bf H}_{1}=\mathbf{0}$ and ${\bf C}\boldsymbol{\Sigma}{\bf G}^{T}{\bf M}=\boldsymbol{\Sigma}{\bf G}^{T}{\bf M}$. Since ${\bf G}{\bf H}_{1}=0$, we have ${\bf C}=\boldsymbol{\Sigma}{\bf G}^{T}({\bf G}\boldsymbol{\Sigma}{\bf G}^{T})^{-1}{\bf G}$. The proof is completed by substituting $(\mathbf{0}_{k\times k_{1}},\boldsymbol{\Sigma}{\bf G}^{T}{\bf M}){\bf H}^{-1}$ in (20) with $\boldsymbol{\Sigma}{\bf G}^{T}({\bf G}\boldsymbol{\Sigma}{\bf G}^{T})^{-1}{\bf G}$. ∎

To solve the problem in (3), one may solve an equivalent problem in (6) by defining an invertible transformation matrix ${\bf H}$ that satisfies ${\bf G}{\bf H}_{1}=\mathbf{0}_{k_{2}\times k_{1}}$. Let us denote $\boldsymbol{\Lambda}={\bf H}^{T}\boldsymbol{\Sigma}^{-1}{\bf H}$ as a precision matrix that can be partitioned as in (7). To simply the problem in (6), we choose the transformation matrix ${\bf H}$ to make $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$ be independent to each other. Since $\boldsymbol{z}$ follows a MVN distribution, $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$ are independent to each other if and only if

Since ${\bf H}_{1}^{T}{\bf G}^{T}=\mathbf{0}$ by definition, to make $\boldsymbol{\Lambda}_{12}=\mathbf{0}$, if and only if we have ${\bf H}_{2}$ as

Thus with ${\bf H}$ satisfying ${\bf G}{\bf H}_{1}=\mathbf{0}_{k_{2}\times k_{1}}$ and ${\bf H}_{2}=\boldsymbol{\Sigma}{\bf G}^{T}{\bf M}$, one can transform the original problem in (3) to that in (6), where $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$ are independent and the restrictions ${\bf G}\boldsymbol{x}=\boldsymbol{r}$ and $\boldsymbol{z}_{2}=({\bf G}{\bf H}_{2})^{-1}\boldsymbol{r}$ imply each other. Following the naive approach shown in Algorithm 1, one can generate $\boldsymbol{x}$ from (3) as follows

Find ${\bf H}=({\bf H}_{1},{\bf H}_{2})$ with ${\bf H}_{2}=\boldsymbol{\Sigma}{\bf G}^{T}{\bf M}$ and with ${\bf H}_{1}$ satisfying ${\bf G}{\bf H}_{1}=\mathbf{0}_{k_{2}\times k_{1}}$;

Sample $\boldsymbol{z}_{1}\sim\mathcal{N}[({\bf I}_{k_{1}},\mathbf{0}_{k_{1}\times k_{2}}){\bf H}^{-1}\boldsymbol{\mu},({\bf H}^{T}_{1}\boldsymbol{\Sigma}^{-1}{\bf H}_{1})^{-1}]$;

Return \boldsymbol{x}={\bf H}\left[\begin{array}[]{c}\boldsymbol{z}_{1}\\ {\bf M}^{-1}({\bf G}\boldsymbol{\Sigma}{\bf G}^{T})^{-1}\boldsymbol{r}\end{array}\right].

However, this naive approach contains intermediate variables that could be computationally expensive to compute. Below we present a method to bypass these intermediate variables. Since the last step could be reexpressed as

which means $\boldsymbol{z}\sim\mathcal{N}[{\bf H}^{-1}\boldsymbol{\mu},{\bf H}^{-1}\boldsymbol{\Sigma}({\bf H}^{-1})^{T}]$, and choose ${\bf M}$ to make

which means ${\bf G}\boldsymbol{\Sigma}{\bf G}^{T}({\bf M}\boldsymbol{z}_{2})={\bf G}{\bf H}\boldsymbol{z}$. Thus we have

In addition, since if $\boldsymbol{z}\sim\mathcal{N}[{\bf H}^{-1}\boldsymbol{\mu},{\bf H}^{-1}\boldsymbol{\Sigma}({\bf H}^{-1})^{T}]$, then $\boldsymbol{y}={\bf H}\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$. Therefore, without the need to compute any intermediate variables, one may use Algorithm 2 to generate $\boldsymbol{x}$ from the hyperplane truncated MVN distribution.

Using the matrix inversion lemma on (8), we have

Using (11), we have ${\bf G}\boldsymbol{\mu}={\bf G}_{1}\boldsymbol{\mu}_{1}+{\bf G}_{2}\boldsymbol{\mu}_{2}=\boldsymbol{r}$. Since ${\bf G}\boldsymbol{x}=\boldsymbol{r}$, we further have ${\bf G}(\boldsymbol{x}-\boldsymbol{\mu})=0$ and hence ${\bf G}_{1}(\boldsymbol{x}_{1}-\boldsymbol{\mu}_{1})=-{\bf G}_{2}(\boldsymbol{x}_{2}-\boldsymbol{\mu}_{2})$. Therefore, given the construction of $\boldsymbol{\mu}_{2}$ as in (11), we can replace the equality constraint ${\bf G}\boldsymbol{x}=\boldsymbol{r}$ on $\boldsymbol{x}$ by requiring $(\boldsymbol{x}_{2}-\boldsymbol{\mu}_{2})=-{\bf G}_{2}^{-1}{\bf G}_{1}(\boldsymbol{x}_{1}-\boldsymbol{\mu}_{1})$. Using this equivelent constraint together with (3), we have

Applying Theorem 1 to Corollary 3, we can generate $\boldsymbol{x}$ with

Thus we can let $\boldsymbol{x}_{1}=\boldsymbol{\mu}_{1}+\boldsymbol{y}^{\prime}_{1}-\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\left(\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{y}^{\prime}_{1}+\boldsymbol{y}_{2}\right)$, where $\boldsymbol{y}^{\prime}_{1}=\boldsymbol{y}_{1}-\boldsymbol{\mu}_{1}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma}_{11})$. ∎

Using the matrix inversion lemma, we have

The proof is completed by using Corollary 4 with $\boldsymbol{\mu}_{1}=\boldsymbol{\mu}_{\beta}$, $\boldsymbol{\Sigma}_{11}={\bf A}^{-1}$, $\boldsymbol{\Sigma}_{12}={\bf A}^{-1}\boldsymbol{\Phi}^{T}$, and $\boldsymbol{\Sigma}_{22}=\boldsymbol{\Omega}^{-1}+\boldsymbol{\Phi}{\bf A}^{-1}\boldsymbol{\Phi}^{T}$. ∎

Computational Complexity

We present the computational complexities of all proposed algorithms in the following tables, where we highlight with bold the lowest complexity in each row.

Brief derivation of SG-MCMC for a simplex-constrained vector

Based on a comprehensive framework for constructing SG-MCMC algorithms in Ma et al. (2015), we have the mini-batch update rule for a global variable $\boldsymbol{z}$ as

For simplicity, we adopt the same specifications that lead to the stochastic gradient Riemannian Langevin dynamics (SGRLD) inference algorithm for simplex-constrained model parameters (Patterson and Teh, 2013; Ma et al., 2015), namely ${\bf D}\left(\boldsymbol{z}\right)={\bf G}{\left(\boldsymbol{z}\right)^{-1}}$, ${\bf Q}\left(\boldsymbol{z}\right)=\boldsymbol{0}$, and ${\hat{}{\bf B}_{t}}=\boldsymbol{0}$, where ${\bf G}\left(\boldsymbol{z}\right)$ denotes the Fisher information matrix (FIM) (Girolami and Calderhead, 2011).

Substituting both (Brief derivation of SG-MCMC for a simplex-constrained vector) and (26) into (Brief derivation of SG-MCMC for a simplex-constrained vector), we have (16) as

Next we prove that equation (16) can be equivalently represented as (17), namely