lrpsadmm: Fit the the low-rank plus sparse estimator using the...

In benjaminfrot/lrpsadmm: Low-Rank Plus Sparse Estimation via Alternating Direction Method of Multipliers (ADMM)

Description

Given a n x p data matrix X and its empirical correlation matrix Σ, an alternating direction method of multipliers (ADMM) algorithm is used to obtain the solutions S, L to

(S, L) = argmin_{A, B} -logdet(A - B) + Tr((A-B)Σ) + λ_1 ||A||_1 + λ_2Tr(B),

subject to A - B > 0 and B ≥ 0.

Usage

lrpsadmm(Sigma, Lambda1, Lambda2, init = NULL, maxiter = 1000,
  mu = 1, abs_tol = 1e-04, rel_tol = 0.01, print_progress = TRUE,
  print_every = 10, zeros = NULL, backend = "RcppEigen")

Arguments

Sigma	An estimate of the correlation matrix.
Lambda1	Penalty on the l1 norm of S
Lambda2	Penalty on the sum of the eigenvalues of L
init	The output of a previous run of the algorithm. For warm starts.
maxiter	Maximal number of iterations
mu	Stepsize of the ADMM algorithm.
abs_tol	Absolute tolerance required to stop the algorithm. Default 1e-04.
rel_tol	Relative tolerance required to stop the algorithm. The algorithm stops when both the change in parameters is below tolerance and the constraints are satisfied. Default 1e-02.
print_progress	Whether the algorithm should report on its progress.
print_every	How often should the algorithm report on its progress (in terms of #iterations).
zeros	A p x p matrix with entries set to 0 or 1. Whereever its entries are 0, the entries of the estimated S will be forced to 0.
backend	A character vector. It should be either 'RcppEigen' or 'R'. If 'R' then use the pure R implementation, if 'RcppEigen' then use the C++ implementation.

Details

Given a n x p data matrix X and its empirical correlation matrix Σ, an alternating direction method of multipliers (ADMM) algorithm is used to obtain the solutions S, L to

(S, L) = argmin_{A, B} -logdet(A - B) + Tr((A-B)Σ) + λ_1 ||A||_1 + λ_2Tr(B),

subject to A - B > 0 and B ≥ 0. This is the estimator suggested in Chandrasekaran et al.

The optimisation problem is decomposed as a three-block ADMM optimisation problem, as described in Ye et al. Because it is a so-called consensus problem, the ADMM is guaranteed to converge.

The tuning parameters λ_1 and λ_2 are typically reparametrised as λ_1 = λ γ and λ_2 = λ (1 - γ), for γ \in (0,1). Here, for a fixed γ, λ controls the overall shrinkage along the path defined by γ. γ controls the tradeoff on the penalties between sparse and low-rank components.

For numerical stability, a smaller value of γ is preferable when n is close to p. See examples.

Value

An S3 object of class lrpsadmm. It is essentially a list with keys:

S

A p x p matrix. The sparse estimate S

L

A p x p matrix. The low-rank estimate L

termcode

An integer. Its value determines whether the algorithm terminated normally or with an error. 0: Convergence reached. -1: Maxiter reached. -2: Shrinkage too strong.

termmsg

A character vector. The message corresponding to the value termcode.

history

A numerical dataframe with the objective function at each iterations, the norm and dual norm as well the primal and dual tolerance criteria for convergence. The algorithm exits when r_norm < eps_pri and s_norm < eps_dual.

U

A p x p matrix. Augmented Lagrangian multiplier. It is stored in order to allow warm starts.

References

Chandrasekaran, Venkat; Parrilo, Pablo A.; Willsky, Alan S. Latent variable graphical model selection via convex optimization. Ann. Statist. 40 (2012), no. 4, 1935–1967. doi:10.1214/11-AOS949. https://projecteuclid.org/euclid.aos/1351602527

Gui-Bo Ye, Yuanfeng Wang, Yifei Chen, Xiaohui Xie; Efficient Latent Variable Graphical Model Selection via Split Bregman Method. https://arxiv.org/abs/1110.3076

See Also

lrpsadmm.cv lrpsadmm.path

Examples

set.seed(1)
# Generate data for a well-powered dataset
sim.data <- generate.latent.ggm.data(n=2000, p=100, h=5, outlier.fraction = 0.0,
                                     sparsity = 0.02, sparsity.latent = 0.7)
X <- sim.data$obs.data; Sigma <- cor(X) # Sample correlation matrix

### Fit the estimator for some value of the tuning parameters
lambda <- 0.7; gamma <- 0.1 # The tuning parameters.
l1 <- lambda * gamma; l2 <- lambda * (1 - gamma)
fit <- lrpsadmm(Sigma = Sigma, Lambda1 = l1, Lambda2 = l2, 
                abs_tol=1e-06, rel_tol=1e-04)
plot(fit) # Use the S3 method plot
estS <- fit$S # Sparse estimate
image(estS!=0) # Visualise its non-zero pattern
estL <- fit$L # Low-rank estimate
plot(eigen(estL)$values) # Visualise the spectrum of the low-rank estimate
plot(fit$history$Objval, type='l') # The log-likelihood from iteration 15 onwards

### Fit for another value of the tuning parameters and compare cold/warm starts
lambda <- 0.4; gamma <- 0.1 # Change the tuning parameters
l1 <- lambda * gamma; l2 <- lambda * (1 - gamma)
# Reuse the previous fit as warm start:
warm.fit <- lrpsadmm(Sigma = Sigma, Lambda1 = l1, Lambda2 = l2,
                     init=fit, rel_tol = 1e-04, abs_tol=1e-06)
# Fit without warm start
cold.fit <- lrpsadmm(Sigma = Sigma, Lambda1 = l1, Lambda2 = l2, rel_tol=1e-04,
                     abs_tol=1e-06)
plot(cold.fit$history$Objval, type='l', col='red', ylab = "Log-Likelihood", xlab = "#Iteration")
lines(warm.fit$history$Objval, col='blue')
xleg = 0.5 * nrow(cold.fit$history)
yleg = 0.5 * (max(cold.fit$history$Objval) + min(cold.fit$history$Objval))
legend(x = xleg, y=yleg, legend=c("Cold Start", "Warm Start"), col=c("red", "blue"), lty=c(1,1))

### Force the sparsity pattern of the sparse component
zeros = 1 * (sim.data$precision.matrix != 0) # A mtrix of 0 and 1.
zeros = zeros[1:100,1:100] # Keep only the observed part
# Whereever zeros[i,j] = 0, the estimated S will be 0.
fit.zeros <- lrpsadmm(Sigma, l1, l2, abs_tol=1e-06, rel_tol=1e-04, zeros = zeros)
fit.no.zeros <- lrpsadmm(Sigma, l1, l2, abs_tol=1e-06, rel_tol=1e-04)
image(fit.zeros$S!=0) # Comparing the sparsity patterns
image(fit.no.zeros$S!=0)

### Fit the estimator when the problem is not so well-posed (n close to p)
set.seed(0)
# n = 80, p = 100 with 5 latent variables.
sim.data <- generate.latent.ggm.data(n=80, p=100, h=5, outlier.fraction = 0.0,
                                     sparsity = 0.02, sparsity.latent = 0.7)
X <- sim.data$obs.data; Sigma <- cor(X) # Sample correlation matrix

lambda <- 2; gamma <- 0.1 # Here gamma is fairly small
l1 <- lambda * gamma; l2 <- lambda * (1 - gamma)
fit.small.gamma <- lrpsadmm(Sigma = Sigma, Lambda1 = l1, Lambda2 = l2)
plot(eigen(fit.small.gamma$L)$value) # Spectrum of L
# This is too high for this ratio n/p
l1 <- lambda * gamma; l2 <- lambda * (1 - gamma)
fit.large.gamma <- lrpsadmm(Sigma = Sigma, Lambda1 = l1, Lambda2 = l2)
plot(eigen(fit.large.gamma$L)$value) # Spectrum of L
# Numerical stability and convergence are not guaranteed for such an ill-posed proble.
# Gamma is too high.
plot(fit.large.gamma$history$Objval, type = 'l', xlab= "#Iterations", ylab = "Log-Likelihood")

### Fit the estimator with a robust estimator of the correlation matrix
# Generate data with 5% of outliers
set.seed(0)
sim.data <- generate.latent.ggm.data(n=2000, p=100, h=5, outlier.fraction = 0.05,
                                     sparsity = 0.02, sparsity.latent = 0.7)
X <- sim.data$obs.data;
Sigma <- cor(X) # Sample correlation matrix
Sigma.Kendall <- Kendall.correlation.estimator(X) # The robust estimator

lambda <- 0.7; gamma <- 0.1 # The tuning parameters.
l1 <- lambda * gamma; l2 <- lambda * (1 - gamma)
# Outliers make the problem very ill-posed
fit <- lrpsadmm(Sigma = Sigma, Lambda1 = l1, Lambda2 = l2,
                abs_tol=1e-06, rel_tol=1e-04, print_every = 200) 
# Use the Kendall based estimator
Kendall.fit <- lrpsadmm(Sigma = Sigma.Kendall, Lambda1 = l1,
                        Lambda2 = l2, abs_tol=1e-06, rel_tol=1e-04) 
image(fit$S!=0)
image(Kendall.fit$S!=0)
plot(fit$history$Objval, xlab="#Iterations", ylab="Log-Likelihood")
plot(Kendall.fit$history$Objval, xlab="#Iterations", ylab="Log-Likelihood")

benjaminfrot/lrpsadmm documentation built on Oct. 19, 2019, 8:13 a.m.