Spike and slab Bayesian sparse principal component analysis

Abstract

Sparse principal component analysis (SPCA) is a popular tool for dimensionality reduction in high-dimensional data. However, there is still a lack of theoretically justified Bayesian SPCA methods that can scale well computationally. One of the major challenges in Bayesian SPCA is selecting an appropriate prior for the loadings matrix, considering that principal components are mutually orthogonal. We propose a novel parameter-expanded coordinate ascent variational inference (PX-CAVI) algorithm. This algorithm utilizes a spike and slab prior, which incorporates parameter expansion to cope with the orthogonality constraint. Besides comparing to two popular SPCA approaches, we introduce the PX-EM algorithm as an EM analogue to the PX-CAVI algorithm for comparison. Through extensive numerical simulations, we demonstrate that the PX-CAVI algorithm outperforms these SPCA approaches, showcasing its superiority in terms of performance. We study the posterior contraction rate of the variational posterior, providing a novel contribution to the existing literature. The PX-CAVI algorithm is then applied to study a lung cancer gene expression dataset. The $\textsf{R}$ package $\textsf{VBsparsePCA}$ with an implementation of the algorithm is available on the Comprehensive R Archive Network (CRAN).

Appendix A: Derivation of (12)–(18)

First, we need the following result:

$$\begin{aligned}&\mathbb {E}_{w|\Theta ^{(t)}} \left[ \frac{1}{2\sigma ^2} \sum _{i=1}^n (X_{ij} - {{\widetilde{\beta }}}_j w_i)^2 \right] \nonumber \\&= \frac{1}{2\sigma ^2} \sum _{i=1}^n \left( X_{ij}^2 - 2X_{ij} {{\widetilde{\beta }}}_j {{\widetilde{\omega }}}_i + {{\widetilde{\beta }}}_j H_i {{\widetilde{\beta }}}_j' \right) , \end{aligned}$$

(36)

where $H_i = {{\widetilde{\omega }}}_i{{\widetilde{\omega }}}_i' + {\widetilde{V}}_w$ and the expressions of ${\widetilde{w}}_i$ and ${\widetilde{V}}_w$ are given in (10).

Since the ELBO is a summation of p terms, we solve $u_j$ and $M_j$ for each j. As the posterior conditional on $\gamma _j =0$ is singular to the Dirac measure, we only need to consider the case $\gamma _j = 1$. This leads to minimize the function

$$\begin{aligned}&\mathbb {E}_{{\widetilde{u}}_j, {\widetilde{M}}_j, z_j|\gamma _j = 1}\Bigg [ \frac{1}{2\sigma ^2} \sum _{i=1}^n \left( - 2X_{ij} {{\widetilde{\beta }}}_j {{\widetilde{\omega }}}_i + {{\widetilde{\beta }}}_j H_i {{\widetilde{\beta }}}_j' \right) \\&\qquad \qquad \qquad + \log \frac{N({\widetilde{u}}_j, \sigma ^2 {\widetilde{M}}_j)}{\kappa _j^\circ g( {{\widetilde{\beta }}}_j|\lambda _1)} \Bigg ]\\ {}&\quad = C - \frac{1}{\sigma ^2} \sum _{i=1}^n X_{ij} {\widetilde{u}}_j {{\widetilde{\omega }}}_i \\&\hspace{1cm}+ \frac{1}{2\sigma ^2} \sum _{i=1}^n \left( {\widetilde{u}}_j H_i {\widetilde{u}}_j' + {{\,\text {Tr}\,}}\left( \sigma ^2{\widetilde{M}}_jH_i\right) \right) \\ {}&\hspace{1cm} + \lambda _1 \sum _{k=1}^r f({\widetilde{u}}_{jk}, {\widetilde{M}}_{j,kk}), \end{aligned}$$

where $\kappa _j^\circ = \int \pi (\gamma _j|\kappa ) d\Pi (\kappa )$. Then we take the derivative of ${\widetilde{u}}_j$ and ${\widetilde{M}}_j$ to obtain (12) and (13). The solutions in (14) are obtained by changing $\lambda _1 \sum _{k=1}^r f({\widetilde{u}}_{jk}, {\widetilde{M}}_{j,kk})$ in the last display with $\frac{\lambda _1}{2\sigma ^2} \left( {\widetilde{u}}_j {\widetilde{u}}_j' + \sigma ^2 {{\,\textrm{Tr}\,}}(M_j)\right)$.

To derive (15), we have

$$\begin{aligned}&\mathbb {E}_P \left( \mathbb {E}_{w|\Theta ^{(t)}} \pi ({{\widetilde{\beta }}}_j, w, X) - \log q({{\widetilde{\beta }}}_j) \right) \nonumber \\&= C + \mathbb {E}_{{{\widetilde{\mu }}}_j, {\widetilde{M}}_j, z_j} \Bigg [ \frac{1}{2\sigma ^2} \sum _{i=1}^n \left( -2X_{ij} {{\widetilde{\beta }}}_j {{\widetilde{\omega }}}_i + {{\widetilde{\beta }}}_j H_i {{\widetilde{\beta }}}_j' \right) \nonumber \\&\quad + \mathbb {1}_{\{\gamma _j = 0\}}\log \frac{1 - z_j}{1-\kappa _j^\circ } \nonumber \\&\quad + \mathbb {1}_{\{\gamma _j = 1\}}\log \frac{z_j N({{\widetilde{\mu }}}_j, \sigma ^2 M_j)}{\kappa _j^\circ g({{\widetilde{\beta }}}_j|\lambda _1)} \Bigg ]\nonumber \\&= C+ (1-z_j) \log \frac{1-z_j}{1-\kappa _j^\circ }\nonumber \\&\quad +z_j \Bigg \{ \frac{1}{2\sigma ^2} \sum _{i=1}^n \left( {{\widetilde{\mu }}}_j H_i {{\widetilde{\mu }}}_j' +\sigma ^2 {{\,\text {Tr}\,}}({\widetilde{M}}_j H_i) - 2X_{ij} {{\widetilde{\mu }}}_j {{\widetilde{\omega }}}_i \right) \nonumber \\&\quad + r\log \left( \frac{\sqrt{2}}{\sqrt{\pi } \sigma \lambda _1}\right) - \frac{1}{2}\log \det ({\widetilde{M}}_j) - \frac{1}{2}\nonumber \\&\quad + \lambda _1 \sum _{k=1}^r f({{\widetilde{\mu }}}_{jk},\sigma ^2 {\widetilde{M}}_{j,kk}) + \log \frac{z_j}{\kappa _j^\circ } \Bigg \}. \end{aligned}$$

(37)

The solution of ${\widehat{h}}_j$ can be obtained by minimizing $z_j$ from the last line of the above display. Similarly, (16) is obtained by minimizing $z_j$ from the following expression

$$\begin{aligned}&C+z_j\Bigg \{ \frac{1}{2\sigma ^2} \sum _{i=1}^n \left( {{\widetilde{\mu }}}_j H_i {{\widetilde{\mu }}}_j' +\sigma ^2 {{\,\text {Tr}\,}}({\widetilde{M}}_j H_i) - 2X_{ij} {{\widetilde{\mu }}}_j {{\widetilde{\omega }}}_i \right) \nonumber \\&- \frac{r\log \lambda _1 + 1}{2}- \frac{1}{2}\log \det ({\widetilde{M}}_j) + \frac{\lambda _1}{2\sigma ^2} \left( {\widetilde{u}}_j {\widetilde{u}}_j' + {{\,\text {Tr}\,}}(\sigma ^2 {\widetilde{M}}_j) \right) \nonumber \\&+ \log \frac{z_j}{\kappa _j^\circ }\Bigg \} + (1-z_j) \log \frac{1-z_j}{1-\kappa _j^\circ }. \end{aligned}$$

(38)

Last, to obtain (17), we first sum the expressions in (37) for all $j =1, \dots , p$. Next, we write down the explicit expression of C which involves $\sigma ^2$, i.e.,

$$\begin{aligned} pC_{\sigma ^2} = \frac{(np + 2\sigma _a +2)\log \sigma ^2}{2} + \frac{{{\,\textrm{Tr}\,}}(X'X) + 2\sigma _b}{2\sigma ^2}. \end{aligned}$$

Last, we plugging the above expression and solve $\sigma ^2$. The solution (18) can be obtained similarly using (38).