Skip to contents

Cluster variables using L2 fusion penalty and simultaneously estimates a gaussian graphical model structure with the addition of L1 sparsity penalty.

Usage

mglasso(
  x,
  lambda1 = 0,
  fuse_thresh = 0.001,
  maxit = NULL,
  distance = c("euclidean", "relative"),
  lambda2_start = 1e-04,
  lambda2_factor = 1.5,
  precision = 0.01,
  weights_ = NULL,
  type = c("initial"),
  compact = TRUE,
  verbose = FALSE
)

Arguments

x

Numeric matrix (\(n x p\)). Multivariate normal sample with \(n\) independent observations.

lambda1

Positive numeric scalar. Lasso penalty.

fuse_thresh

Positive numeric scalar. Threshold for clusters fusion.

maxit

Integer scalar. Maximum number of iterations.

distance

Character. Distance between regression vectors with permutation on symmetric coefficients.

lambda2_start

Numeric scalar. Starting value for fused-group Lasso penalty (clustering penalty).

lambda2_factor

Numeric scalar. Step used to update fused-group Lasso penalty in a multiplicative way..

precision

Tolerance for the stopping criterion (duality gap).

weights_

Matrix of weights.

type

If "initial" use classical version of MGLasso without weights.

compact

Logical scalar. If TRUE, only save results when previous clusters are different from current.

verbose

Logical scalar. Print trace. Default value is FALSE.

Value

A list-like object of class mglasso is returned.

out

List of lists. Each element of the list corresponds to a clustering level. An element named levelk contains the regression matrix beta and clusters vector clusters for a clustering in k clusters. When compact = TRUE out has as many elements as the number of unique partitions. When set to FALSE, the function returns as many items as the the range of values taken by lambda2.

l1

the sparsity penalty lambda1 used in the problem solving.

Details

Estimates a gaussian graphical model structure while hierarchically grouping variables by optimizing a pseudo-likelihood function combining Lasso and fuse-group Lasso penalties. The problem is solved via the COntinuation with NEsterov smoothing in a Shrinkage-Thresholding Algorithm (Hadj-Selem et al. 2018). Varying the fusion penalty \(\lambda_2\) in a multiplicative fashion allow to uncover a seemingly hierarchical clustering structure. For \(\lambda_2 = 0\), the approach is theoretically equivalent to the Meinshausen-Bühlmann (2006) neighborhood selection as resuming to the optimization of pseudo-likelihood function with \(\ell_1\) penalty (Rocha et al., 2008). The algorithm stops when all the variables have merged into one cluster. The criterion used to merge clusters is the \(\ell_2\)-norm distance between regression vectors.

For each iteration of the algorithm, the following function is optimized : $$1/2 \sum_{i=1}^p || X^i - X^{\ i} \beta^i ||_2 ^2 + \lambda_1 \sum_{i = 1}^p || \beta^i ||_1 + \lambda_2 \sum_{i < j} || \beta^i - \tau_{ij}(\beta^j) ||_2.$$

where \(\beta^i\) is the vector of coefficients obtained after regression \(X^i\) on the others and \(\tau_{ij}\) is a permutation exchanging \(\beta_j^i\) with \(\beta_i^j\).

See also

conesta for the problem solver, plot_mglasso for plotting the output of mglasso.

Examples

# \donttest{
install_conesta()
#> mglasso requires the r-reticulate virtual environment. Attempting to create...
#> virtualenv: r-reticulate
#> Using virtual environment 'r-reticulate' ...
#> + '/home/runner/.virtualenvs/r-reticulate/bin/python' -m pip install --upgrade 'scipy == 1.7.1' 'scikit-learn'
#> Installing pylearn-parsimony
#> pylearn-parsimony is installed.
n = 50
K = 3
p = 9
rho = 0.85
blocs <- list()
for (j in 1:K) {
  bloc <- matrix(rho, nrow = p/K, ncol = p/K)
  for(i in 1:(p/K)) { bloc[i,i] <- 1 }
  blocs[[j]] <- bloc
}

mat.covariance <- Matrix::bdiag(blocs)
mat.covariance
#> 9 x 9 sparse Matrix of class "dgCMatrix"
#>                                                   
#>  [1,] 1.00 0.85 0.85 .    .    .    .    .    .   
#>  [2,] 0.85 1.00 0.85 .    .    .    .    .    .   
#>  [3,] 0.85 0.85 1.00 .    .    .    .    .    .   
#>  [4,] .    .    .    1.00 0.85 0.85 .    .    .   
#>  [5,] .    .    .    0.85 1.00 0.85 .    .    .   
#>  [6,] .    .    .    0.85 0.85 1.00 .    .    .   
#>  [7,] .    .    .    .    .    .    1.00 0.85 0.85
#>  [8,] .    .    .    .    .    .    0.85 1.00 0.85
#>  [9,] .    .    .    .    .    .    0.85 0.85 1.00

set.seed(11)
X <- mvtnorm::rmvnorm(n, mean = rep(0,p), sigma = as.matrix(mat.covariance))
X <- scale(X)

res <- mglasso(X, 0.1, lambda2_start = 0.1)
res$out[[1]]$clusters
#> [1] 1 2 3 4 5 6 7 8 9
res$out[[1]]$beta
#>              [,1]        [,2]         [,3]        [,4]         [,5]        [,6]
#>  [1,]  0.00000000  0.50942304  0.420856249  0.08059251  0.198555576 -0.15744283
#>  [2,]  0.57897509  0.00000000  0.375480194 -0.11129404 -0.024338821  0.03622193
#>  [3,]  0.50877541  0.39402762  0.000000000  0.09138873 -0.244767942  0.13197185
#>  [4,]  0.09644133 -0.11232346  0.096386587  0.00000000  0.380559101  0.53949430
#>  [5,]  0.34031463 -0.03189842 -0.367833554  0.55046923  0.000000000  0.28792248
#>  [6,] -0.23822144  0.05604748  0.156337318  0.67569926  0.247844203  0.00000000
#>  [7,]  0.07147603  0.00000000 -0.066314827  0.19093102  0.005503117 -0.21086078
#>  [8,] -0.17121934  0.13577447  0.006274291  0.04037628  0.017037923  0.01709042
#>  [9,]  0.00000000 -0.03662807  0.030397589 -0.27537816  0.024525025  0.21891289
#>               [,7]         [,8]        [,9]
#>  [1,]  0.088104950 -0.136441557  0.00000000
#>  [2,] -0.005820201  0.112868548 -0.04462052
#>  [3,] -0.078583570  0.004854681  0.03876315
#>  [4,]  0.274032064  0.031338363 -0.32834722
#>  [5,]  0.011236774  0.027884703  0.06051871
#>  [6,] -0.374671576  0.008336123  0.32264027
#>  [7,]  0.000000000  0.389473752  0.56690235
#>  [8,]  0.592191132  0.000000000  0.30687389
#>  [9,]  0.688368824  0.245147593  0.00000000
# }