Group-sparse and sparse PCA of numerical data

Introduction

The function groupSparsePCA() in the sparsePCA package implements group-sparse principal component analysis (PCA) using either a block optimization algorithm or a deflation-based algorithm. Sparse PCA is recovered as the special case where each variable forms its own group (the default when groups = NULL). Group sparsity means that the loadings associated with predefined groups of variables are simultaneously set to zero. The variance of the potentially nonorthogonal principal components is computed using the optimal projected variance.

library(sparsePCA)

A small simulated data example

The function simuPCA() implements the simulation procedure described in Shen and Huang (2008) to generate data from an underlying PCA model. A data matrix \({\bf X}\) of dimension \(n \times p\) is randomly drawn from a centered multivariate normal distribution whose covariance matrix is constructed from specified eigenvalues and leading eigenvectors.

Here, we randomly generate a data matrix \({\bf X}\) with:

• \(p = 10\) eigenvalues \((200, 100, 50, 50, 6, 5, 4, 3, 2, 1)\);

• \(m = 2\) leading eigenvectors:

  • \(v_1 = (1, 1, 1, 1, 0, 0, 0, 0, 0.9, 0.9)\);
  • \(v_2 = (0, 0, 0, 0, 1, 1, 1, 1, -0.3, 0.3)\).

These two loading vectors are sparse, with their zero entries corresponding to three predefined groups of variables of sizes 4, 4, and 2.

# True loading vectors
v1 <- c(1, 1, 1, 1, 0, 0, 0, 0, 0.9, 0.9)
v2 <- c(0, 0, 0, 0, 1, 1, 1, 1, -0.3, 0.3)
norms <- sqrt(colSums(cbind(v1, v2)^2))
Vtrue <- sweep(cbind(v1, v2), 2, norms, "/") # true normalized loadings

# True eigenvalues
eigvals <- c(200, 100, 50, 50, 6, 5, 4, 3, 2, 1)

Let us now generate a data matrix \({\bf X}\) with \(n = 50\) observations from this PCA model.

# Simulation of a numerical data matrix

set.seed(1)
X <- simuPCA(n = 100, 
             Vtrue = cbind(v1, v2), 
             eigvals = eigvals)

Standard PCA (based on the covariance matrix) is first performed. Here, the loadings are computed by singular value decomposition of the centered data matrix \({\bf X}\).

# Center the data
Z <- scale(X, center = TRUE, scale = FALSE)

# PCA loadings
V_pca <- svd(Z)$v[, 1:2]  

The PCA loadings are close to the true loading vectors defined above, but the zero loadings are not identified. To promote sparsity in the loadings, sparse PCA is applied with the following arguments:

• \(m = 2\) components;
• \(\lambda = (0.1, 0.1)\), corresponding to a common sparsity parameter of \(0.1\) for both components;
• scale = FALSE, meaning that sparse PCA is performed on the covariance matrix (i.e., the variables are centered but not scaled).

# Sparse PCA loadings
out <- groupSparsePCA(X = X, 
                 m = 2,
                 lambda = c(0.1, 0.1), 
                 scale = FALSE)
V_spca <- out$V

# Selected variables
out$sel 

## $dim1
## [1] "X1"  "X2"  "X3"  "X4"  "X6"  "X7"  "X9"  "X10"
## 
## $dim2
## [1] "X1"  "X2"  "X5"  "X6"  "X7"  "X8"  "X9"  "X10"

The sparse PCA loadings recover some of the true zero loadings, but not all of them. To promote group sparsity in the loadings, group-sparse PCA is finally applied with the following arguments:

• \(m = 2\) components;
• \(\lambda = (0.3, 0.3)\);
• scale = FALSE.

In addition, the following arguments are specified:

• groups = c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3), indicating the group membership of the 10 variables, with three groups: G1 = {X1, X2, X3, X4}, G2 = {X5, X6, X7, X8}, and G3 = {X9, X10}.
• groupnames = c("G1", "G2", "G3"), specifying the names of the three groups.
• groupsize = FALSE, indicating that group sizes are not taken into account in the thresholding operator.

# Group-sparse PCA loadings
out <- groupSparsePCA(X, 
                 m = 2, 
                 lambda = c(0.3,0.3), 
                 groups = rep(c(1,2,3), c(4,4,2)),
                 scale=F,
                 groupnames = c("G1", "G2", "G3"),
                 groupsize = FALSE)
V_gspca <- out$V

# Selected groups
out$sel 

## $dim1
## [1] "G1" "G3"
## 
## $dim2
## [1] "G2"

The groups of variables with zero loadings are now well recovered, except for group G3, which is not selected for the second component.

Block or deflation approach

The groupSparsePCA() function implements group-sparse and sparse PCA using either a block optimization algorithm or a deflation-based algorithm.

The block algorithm uses a numeric vector of strictly positive weights \(\mu_k\), \(k = 1, \ldots, m\), usually chosen to be either strictly decreasing (\(\mu_k = 1/k\), \(k = 1, \ldots, m\)) or constant (\(\mu_k = 1\) for all \(k\)). Strictly decreasing weights recover the standard PCA solution when all sparsity parameters are set to zero. With constant weights, the standard PCA solution is not necessarily recovered; instead, an orthogonal rotation of the PCA solution may be obtained.

In practice, three algorithms can be used to compute group-sparse or sparse PCA:

• If block = FALSE, a deflation-based algorithm is used to recursively compute the loading vectors.
• If block = TRUE and mu = 1/(1:m), a block algorithm with decreasing weights is used. It recovers the standard PCA solution when lambda = c(0, ..., 0). These are the default parameters of the function.
• If block = TRUE and mu = c(1, ..., 1), a block algorithm with equal weights is used. It does not necessarily recover the standard PCA solution when lambda = c(0, ..., 0), but rather an orthogonal rotation of it. However, this algorithm directly maximizes the penalized optimal projected variance, rather than a weighted version of it.

out <- groupSparsePCA(X = X, 
                 m = 2,
                 lambda = c(0.3, 0.3), 
                 groups = rep(c(1,2,3), c(4,4,2)),
                 block = FALSE)
out$V

##           dim1      dim2
## X1  -0.4431646 0.0000000
## X2  -0.4322402 0.0000000
## X3  -0.4712233 0.0000000
## X4  -0.4539533 0.0000000
## X5   0.0000000 0.4633705
## X6   0.0000000 0.5144252
## X7   0.0000000 0.5325856
## X8   0.0000000 0.4868338
## X9  -0.3003683 0.0000000
## X10 -0.3137314 0.0000000

out <- groupSparsePCA(X = X, 
                 m = 2,
                 lambda = c(0.3, 0.3), 
                 groups = rep(c(1,2,3), c(4,4,2)),
                 block = TRUE,
                 mu = c(1,1))
out$V

##           dim1      dim2
## X1  -0.4441751 0.0000000
## X2  -0.4332577 0.0000000
## X3  -0.4715423 0.0000000
## X4  -0.4555815 0.0000000
## X5   0.0000000 0.4625772
## X6   0.0000000 0.5113720
## X7   0.0000000 0.5294358
## X8   0.0000000 0.4941848
## X9  -0.2934401 0.0000000
## X10 -0.3146148 0.0000000

out <- groupSparsePCA(X = X, 
                 m = 2,
                 lambda = c(0.3, 0.3), 
                 groups = rep(c(1,2,3), c(4,4,2)),
                 block = TRUE,
                 mu = c(1, 1/2))
out$V 

##           dim1      dim2
## X1  -0.4435191 0.0000000
## X2  -0.4325931 0.0000000
## X3  -0.4713464 0.0000000
## X4  -0.4545482 0.0000000
## X5   0.0000000 0.4621565
## X6   0.0000000 0.5098144
## X7   0.0000000 0.5278418
## X8   0.0000000 0.4978793
## X9  -0.2978345 0.0000000
## X10 -0.3141133 0.0000000

A real data example

The protein dataset is described in Hand et al., A Handbook of Small Data Sets (1994, p. 297). It contains protein consumption from nine food groups in 25 European countries.

data(protein, package = "PCAmixdata")
?PCAmixdata::protein

Let us first apply standard PCA based on the correlation matrix. Here, PCA is performed using the singular value decomposition of the standardized data matrix \({\bf Z}\). The squared singular values (properly normalized) are the eigenvalues of the correlation matrix and correspond to the variances of the principal components.

n <- nrow(protein)
Z <- scale(protein)*sqrt(n/(n-1))
eig <- svd(Z)$d^2/n

The barplot of the eigenvalues suggest to choose \(m=4\) components.

barplot(eig, names = 1:9, xlab = "dim", ylab = "pev (in %)")

# Proportion of variance explained by the principal components
pev <- eig/sum(eig) 

The first four principal components explain 85.8% of the total variance of the standardized data. The proportions of variance explained by each component are as follows:

The PCA loading vectors below show that some variables are not essential for building the principal components, as their loadings are close to zero. For instance, the variable Red.Meat is almost irrelevant for the construction of the second principal component, with a loading of -0.056.

# PCA loading matrix

Vpca <- svd(Z)$v[,1:4] 

In order to promote sparsity, sparse PCA (using the block algorithm with decreasing weights) is applied to the protein dataset with \(m = 4\) components.

Once the number of components \(m\) has been selected, a common sparsity parameter \(\lambda\) must be chosen. For this purpose, we jointly examine the regularization paths of the proportion of variance explained and the number of selected variables.

A one-dimensional grid of common values of the sparsity parameter is first defined: \[0 \leq \lambda = \lambda_1 = \ldots = \lambda_4 \leq 1,\] by letting \(\lambda\) vary from 0 to 1 in steps of 0.01.

grid <- seq(0, 1, by = 0.01) # grid of common lambda values

The loading matrices are computed for each value in this one-dimensional grid, together with the corresponding proportions of variance explained (computed using the optimal projected variance) and the number of selected variables for each component.

m <- 4
pev <- matrix(NA, length(grid), m) # proportion of explained variance
nsel <- matrix(NA, length(grid), m) # number of selected variables

for (k in 1:length(grid))
{
  out <- groupSparsePCA(protein, m = m, lambda = rep(grid[k], m)) 
  pev[k,] <- out$pev
  nsel[k,] <- out$nsel
}

colnames(pev) <- colnames(nsel) <- paste0("dim", 1:m)

The figure below shows the proportion of variance explained by the principal components (pev) as a function of the common sparsity parameter \(\lambda\).

kmax <- 80
matplot(grid[1:kmax], 
        pev[1:kmax,], 
        xlab = "lambda", 
        ylab = "pev",
        type = "l", 
        col = 1:m, 
        lty = 1) 

legend("topright",
       legend=colnames(pev),
       lty = 1,
       col = 1:m,
       bty = "n")

abline(v = 0.6, lty = 2)

The next figure shows the number of selected variables (nsel) as a function of the common sparsity parameter \(\lambda\).

matplot(grid[1:kmax],
        nsel[1:kmax,], 
        xlab = "lambda", 
        ylab = "nsel",
        type = "l",
        col = 1:m,
        lty = 1) 

legend("topright",
       legend=colnames(nsel),
       lty = 1,
       col = 1:m,
       bty = "n")
abline(v = 0.6, lty = 2)

These two regularization paths suggest using \(\lambda=0.6\) (vertical dotted line). The sparse loadings obtained with this value of \(\lambda\) are given below.

out <- groupSparsePCA(protein, lambda=rep(0.6, m), m = m)
out$V # sparse PCA loadings

##                     dim1 dim2       dim3 dim4
## Red.Meat       0.0000000    0  0.0000000   -1
## White.Meat     0.0000000    0  0.5230873    0
## Eggs          -0.5267939    0  0.0000000    0
## Milk          -0.3581631    0  0.0000000    0
## Fish           0.0000000    0 -0.8522791    0
## Cereals        0.6295215    0  0.0000000    0
## Starchy.Foods  0.0000000    0  0.0000000    0
## Nuts           0.4448708    0  0.0000000    0
## Fruite.veg.    0.0000000    1  0.0000000    0

out$pev # proportion of explained variable

##       dim1       dim2       dim3       dim4 
## 0.30138587 0.10632326 0.13267270 0.08862838

out$sel # selected variables

## $dim1
## [1] "Eggs"    "Milk"    "Cereals" "Nuts"   
## 
## $dim2
## [1] "Fruite.veg."
## 
## $dim3
## [1] "White.Meat" "Fish"      
## 
## $dim4
## [1] "Red.Meat"

When \(\lambda = 0.6\), four variables are selected to build the first principal component, one variable is selected for the second principal component, two variables are selected for the third principal component, and one variable is selected for the fourth principal component. The four principal components obtained from these sparse loadings explain 62.9% of the variance, compared with 85.8% for standard PCA.

The principal components \({\bf y}_k = {\bf Z} v_k\) are new synthetic variables computed using only the variables with nonzero loadings, with coefficients derived from the corresponding sparse loading vectors:

\[{\bf y}_1=-0.527 \; Eggs -0.358 \; Milk + 0.630 \; Cereals + 0.445 \; Nuts\] \[{\bf y}_2=Fruit.veg\] \[{\bf y}_3=0.523 \; White.Meat - 0.852 \; Fish\]

\[{\bf y}_4=-Red.Meat\]

head(out$scores) # principal component scores

##              dim1        dim2       dim3        dim4
## Alban  2.73803997 -1.37825141  0.1052588 -0.08294065
## Aust  -1.45349534  0.09278868  1.4404957  0.28297397
## Belg  -1.13284250 -0.07694671  0.1476879 -1.11969872
## Bulg   2.67116619  0.03621022  0.5144084  0.61839570
## Czech -0.02808723 -0.07694671  1.0903050  0.03903089
## Den   -1.92153870 -0.98220215 -1.0160337 -0.23540507

plot(out) # plot of the 25 countries of the two first dimensions.

The predict method uses the coefficients of the principal components described above to predict the scores of new observations.

Remark: the function optvar computes the optimal projected variance of a matrix \({\bf Y}\) containing not necessarily orthogonal principal components. The proportion of explained variance (PEV) is then defined as the optimal projected variance of the components divided by the total variance of the data (here, the variance of the standardized data matrix \({\bf Z}\)).

# Total variance of Z is the number p = 9 of variables
tot <- sum((apply(Z, 2, var) * (n-1)/n)) 

# The principal components are correlated
cor(out$scores)

##             dim1        dim2        dim3        dim4
## dim1  1.00000000  0.21964918 -0.05405997  0.57483240
## dim2  0.21964918  1.00000000 -0.23549355  0.07422123
## dim3 -0.05405997 -0.23549355  1.00000000 -0.02554245
## dim4  0.57483240  0.07422123 -0.02554245  1.00000000

# Optimal variance of the principal components
optvar(out$scores)

##      dim1      dim2      dim3      dim4 
## 2.7124729 0.9569093 1.1940543 0.7976554

optvar(out$scores)/tot # pev

##       dim1       dim2       dim3       dim4 
## 0.30138587 0.10632326 0.13267270 0.08862838

out$pev

##       dim1       dim2       dim3       dim4 
## 0.30138587 0.10632326 0.13267270 0.08862838