spca.cavi.Laplace {VBsparsePCA}

Function for the PX-CAVI algorithm using the Laplace slab

Description

This function employs the PX-CAVI algorithm proposed in Ning (2020). The g in the slab density of the spike and slab prior is chosen to be the Laplace density, i.e., N(0, \sigma^2/\lambda_1 I_r). Details of the model and the prior can be found in the Details section in the description of the 'VBsparsePCA()' function. This function is not capable of handling the case when r > 1. In that case, we recommend to use the multivariate distribution instead.

Usage

spca.cavi.Laplace(
  x,
  r = 1,
  lambda = 1,
  max.iter = 100,
  eps = 0.001,
  sig2.true = NA,
  threshold = 0.5,
  theta.int = NA,
  theta.var.int = NA,
  kappa.para1 = NA,
  kappa.para2 = NA,
  sigma.a = NA,
  sigma.b = NA
)

Arguments

x

Data an n*p matrix.

r

Rank.

lambda

Tuning parameter for the density g.

max.iter

The maximum number of iterations for running the algorithm.

eps

The convergence threshold; the default is 10^{-4}.

sig2.true

The default is false, \sigma^2 will be estimated; if sig2 is known and its value is given, then \sigma^2 will not be estimated.

threshold

The threshold to determine whether \gamma_j is 0 or 1; the default value is 0.5.

theta.int

The initial value of theta mean; if not provided, the algorithm will estimate it using PCA.

theta.var.int

The initial value of theta.var; if not provided, the algorithm will set it to be 1e-3*diag(r).

kappa.para1

The value of \alpha_1 of \pi(\kappa); default is 1.

kappa.para2

The value of \alpha_2 of \pi(\kappa); default is p+1.

sigma.a

The value of \sigma_a of \pi(\sigma^2); default is 1.

sigma.b

The value of \sigma_b of \pi(\sigma^2); default is 2.

Value

iter

The number of iterations to reach convergence.

selection

A vector (if r = 1 or with the jointly row-sparsity assumption) or a matrix (if otherwise) containing the estimated value for \boldsymbol \gamma.

theta.mean

The loadings matrix.

theta.var

The covariance of each non-zero rows in the loadings matrix.

sig2

Variance of the noise.

obj.fn

A vector contains the value of the objective function of each iteration. It can be used to check whether the algorithm converges

Examples

#In this example, the first 20 rows in the loadings matrix are nonzero, the rank is 1
set.seed(2021)
library(MASS)
library(pracma)
n <- 200
p <- 1000
s <- 20
r <- 1
sig2 <- 0.1
# generate eigenvectors
U.s <- randortho(s, type = c("orthonormal"))
U <- rep(0, p)
U[1:s] <- as.vector(U.s[, 1:r])
s.star <- rep(0, p)
s.star[1:s] <- 1
eigenvalue <- seq(20, 10, length.out = r)
# generate Sigma
theta.true <- U * sqrt(eigenvalue)
Sigma <- tcrossprod(theta.true) + sig2*diag(p)
# generate n*p dataset
X <- t(mvrnorm(n, mu = rep(0, p), Sigma = Sigma))
result <- spca.cavi.Laplace(x = X, r = 1)
loadings <- result$theta.mean