An Introduction to amanpg

Algorithm Description

The A-ManPG algorithm can be applied to solve the general manifold optimization problem

where H(A, B) is a smooth function of A, B with a Lipschitz continuous gradient, f(⋅) and g(⋅) are lower semicontinuous (possibly nonsmooth) convex functions, and ℳ₁, ℳ₂ are two embedded submanifolds in the Euclidean space.

For sparse PCA, the following function definitions are used:

• H(A, B) = Tr(B^TX^TXB) − 2Tr(A^TX^TXB)
• f(A) ≡ 0
• $g(B)=\lambda_2\sum_{j=1}^k\|B_j\|_2^2+\sum_{j=1}^k\lambda_{1,j}\|B_j\|_1$
• ℳ₁ = St(p, k)
• ℳ₂ = ℝ^p × k

where X is the n × p data matrix or n × n covariance matrix, A is the scores and B is the loadings, k is the rank of the matrices (in other words, how many principal components are desired), λ₁ is the L1 norm penalty and λ₂ is the L2 norm penalty. Both the L1 and L2 norm are used as elastic net regularization to impose sparseness within the loadings. Note that a different L1 norm penalty is used for every principal component, and the algorithm operates differently when the L2 norm penalty is set to a large constant (np.inf or Inf).

The A-ManPG algorithm uses the following subproblems with an alternating updating scheme to solve sparse PCA, computed in a Gauss-Seidel manner for faster convergence.

where A_k + 1 is obtained via a retraction operation (in this case, polar decomposition), t₁ ≤ L_A, and t₂ ≤ L_B. L_A and L_B are the least upper bounds of the Lipschitz constants for ∇_AH(A, B) and ∇_BH(A, B), respectively. The subproblems are solved using an adaptive semismooth Newton method.

Convergence

Let ϵ represent a tolerance level to detect convergence. An ϵ-stationary point is defined as a point (A, B) with corresponding D^A and D^B that satisfy the following:

The algorithm reaches an ϵ-stationary point in at most

$$\frac{2(F(A_0,B_0)-F^*)}{ ((\gamma\bar{\alpha}_1t_1+\gamma\bar{\alpha}_2t_2)\epsilon^2)}$$

iterations, where:

• (A₀, B₀) are the initial values
• F^* is the lower bound of F, from the general manifold optimization problem
• ᾱ₁ and ᾱ₂ are positive constants

Pseudocode

The following describes the algorithm used for solving the general manifold optimization problem using A-ManPG. For solving sparse PCA, the algorithm is implemented with the aforementioned definitions.

Input initial point (A0,B0) and necessary parameters for the required problem 

for i=0,1,... do
  Solve the first subproblem for Da
  Set alpha = 1
  
  while F(Retr(alpha * Da),B) > F(A,B) - alpha * norm(Da)^2 / (2 * t1) do
    alpha = gamma * alpha
  end while
  
  Set A = Retr(alpha * Da)
  
  Solve the second subproblem for Db
  Set alpha = 1
  
  while F(A,Retr(alpha * Db)) > F(A,B) - alpha * norm(Db)^2 / (2 * t2) do
    alpha = gamma * alpha
  end while
  
  Set B = Retr(alpha * Db)
end for

Return A as the scores and B as the sparse loadings

R Usage

spca.amanpg(z, lambda1, lambda2, 
            f_palm = 1e5, x0 = NULL, y0 = NULL, k = 0, type = 0, 
            gamma = 0.5, maxiter = 1e4, tol = 1e-5, 
            normalize = TRUE, verbose = FALSE)

Python Usage

spca(z, lambda1, lambda2, 
     x0=None, y0=None, k=0, gamma=0.5, type=0, 
     maxiter=1e4, tol=1e-5, f_palm=1e5,
       normalize=True, verbose=False):

Arguments

Values

Python returns a dictionary with the following key-value pairs, while R returns a list with the following elements:

R Example

As with other libraries, begin by loading amanpg in R.

library(amanpg)

Before proceeding, it is helpful to determine a few parameters. Let the rank of the sparse loadings matrix be k = 4 (returning four principal components), the input data matrix be n × p where n = 1000 and p = 500, λ₁ be a 4 × 1 “matrix” where λ_i, 1 = 0.1, and λ₂ = 1.

# parameter initialization
k <- 4
n <- 1000
p <- 500
lambda1 <- matrix(data=0.1, nrow=k, ncol=1)
lambda2 <- 1

For this example, the data matrix z is randomly generated from the normal distribution. Although it should be centered to mean 0 and normalized to Euclidean length 1, the function will automatically preprocess the input matrix when normalize=TRUE.

# data matrix generation
set.seed(10)
z <- matrix(rnorm(n * p), n, p)

# only show a subset of the data matrix for brevity
knitr::kable(as.data.frame(z)[1:10,1:4])

Alternatively, the data can be normalized beforehand and the parameter is set to FALSE in the function call. However, this example won’t do so, but the output of normalize is displayed below.

# see the effects of normalize()
knitr::kable(as.data.frame(normalize(z))[1:10,1:4])

Now the function is called, passing through matrix a, lambda1, lambda2, and our desired rank k. The output is stored as a list in fin_sprout. Note that if a different initial point is desired, x0 and y0 should be modified, but the default value as the first k right singular vectors is sufficient for this example.

If further printout is desired, set verbose=TRUE for progress updates (time, difference for convergence, F value) per iteration.

# function call
fin_sprout <- spca.amanpg(z, lambda1, lambda2, k=4)
print(paste(fin_sprout$iter, "iterations,", fin_sprout$sparsity, "sparsity,", fin_sprout$time))

## [1] "280 iterations, 0.491 sparsity, 0.43851900100708"

The loadings can be viewed from fin_sprout$loadings. Note that many entries are set to zero as a result of the induced sparsity.

# View loadings. Only first 10 rows for brevity
knitr::kable(as.data.frame(fin_sprout$loadings)[1:10,])

The resulting scree plot (Figure 1) looks like this.

pr.var <- (apply(fin_sprout$x, 2, sd))^2
pve <- pr.var / sum(pr.var)

par(mfrow=c(1,2))
plot(pve, 
     xlab="Sparse PC", 
     ylab="Proportion of Variance Explained", 
     ylim=c(0,1), 
     type="b")
plot(cumsum(pve), 
     xlab="Sparse PC", 
     ylab="Cumulative Proportion of Variance Explained", 
     ylim=c(0,1),
     type="b")

Scree plots for the finite lambda case.

The resulting biplot (Figure 2), with the zero loadings filtered out, can be obtained using the following:

y_sub = apply(fin_sprout$loadings, 1, function(row) all(row != 0))
loadings = fin_sprout$loadings[y_sub, ]

par(mfrow=c(1,1))
biplot(fin_sprout$x, loadings, xlab="PC 1", ylab="PC 2")

Biplot for the finite lambda case.

Now consider an alternative situation where we set λ₂ to a large constant Inf. The algorithm changes by directly retracting B without using a while loop to determine an appropriate retraction step size and only iterating A.

# infinite lambda2
inf_sprout <- spca.amanpg(z, lambda1, lambda2=Inf, k=4) 
print(paste(inf_sprout$iter, "iterations,", inf_sprout$sparsity, "sparsity,", inf_sprout$time))

## [1] "344 iterations, 0.253 sparsity, 0.371917009353638"

# extract loadings. Only first 10 rows for brevity
knitr::kable(as.data.frame(inf_sprout$loadings)[1:10,])

We obtain the scree plot (Figure 3) and the biplot (Figure 4) using the same method.

pr.var <- (apply(inf_sprout$x, 2, sd))^2
pve <- pr.var / sum(pr.var)

par(mfrow=c(1,2))
plot(pve, 
     xlab="Sparse PC", 
     ylab="Proportion of Variance Explained", 
     ylim=c(0,1), 
     type="b")
plot(cumsum(pve), 
     xlab="Sparse PC", 
     ylab="Cumulative Proportion of Variance Explained", 
     ylim=c(0,1),
     type="b")

Scree plot for lambda=inf case.

For data with high dimensionality, the biplot is still much harder to read with to a lower sparsity value.

y_sub = apply(inf_sprout$loadings, 1, function(row) all(row != 0))
loadings = inf_sprout$loadings[y_sub, ]

par(mfrow=c(1,1))
biplot(inf_sprout$x, loadings, xlab="PC 1", ylab="PC 2")

Biplot for lambda=inf case. Observe that with lower sparsity in the loadings and high-dimensional data, the biplot becomes less readable.

Python Example

Note that the Python package depends on numpy.

The following example accomplishes the same situation (down to the same randomly-generated data) but in Python.

import numpy as np
from sparsepca import spca

k = 4  # rank
p = 500  # dimensions
n = 1000  # sample size
lambda1 = 0.1 * np.ones((k, 1))
lambda2 = 1

np.random.seed(10)
z = np.random.normal(0, 1, size=(n, p))  # generate random normal 1000x500 matrix

fin_sprout = spca(z, lambda1, lambda2, k=k)
print(f"Finite: {fin_sprout['iter']} iterations with final value 
        {fin_sprout['f_manpg']}, sparsity {fin_sprout['sparsity']}, 
        timediff {fin_sprout['time']}.")

fin_sprout['loadings']

inf_sprout = spca(z, lambda1, np.inf, k=k)
print(f"Infinite: {inf_sprout['iter']} iterations with final value 
        {inf_sprout['f_manpg']}, sparsity {inf_sprout['sparsity']}, 
        timediff {inf_sprout['time']}.")

inf_sprout['loadings']