Sparse estimation of a precision matrix by the graphical Lasso, that is the
minimization over positive definite matrices \(\Theta\) of
$$-\log\det(\Theta) + \mathrm{tr}(S\Theta) + \|\rho \circ \Theta\|_1.$$
This is the solver used internally by PLNnetwork() and ZIPLNnetwork().
graphical_lasso(
S,
rho,
thr = 1e-04,
maxit = 10000L,
w_init = NULL,
wi_init = NULL,
trace = FALSE,
stall_patience = 1000L
)a list with components
w: the estimated covariance matrix,
wi: the estimated precision matrix (symmetric),
niter: the number of outer sweeps performed,
converged: TRUE if the convergence criterion was met,
status: how the solve ended, one of "converged", "stalled",
"max_iter", "inner_failure" or "degenerate" (see Details),
delta: the best value reached by the convergence criterion, relative to
the threshold it had to cross. delta <= 1 means convergence; a stalled
solve typically sits between 1 and 3, that is, just short of it,
dw_trace: the per-sweep criterion when trace = TRUE.
a symmetric p x p (empirical) covariance matrix.
the penalty: either a non-negative scalar, applied to all entries (the diagonal included), or a symmetric p x p matrix of non-negative per-entry penalties (e.g. with a zero diagonal to leave it unpenalized).
convergence threshold, relative to the average absolute
off-diagonal entry of S. Default is 1e-4, as in glassoFast.
maximal number of outer sweeps. Default is 10000, as in
glassoFast.
optional warm start: the w and wi of a previous
solve, typically at a nearby penalty along a regularization path. Both must
be given, with the dimensions of S, to be used. Note that a warm start
stops closer to the starting point than a cold one at the same thr, so it
does not reproduce a cold solve exactly.
if TRUE, also return the per-sweep convergence criterion, to
diagnose a solve that does not converge. Default is FALSE.
number of consecutive sweeps without progress after
which the algorithm concludes that it is cycling and stops (see Details).
Inf disables the detection. Default is 1000.
The algorithm is the block coordinate descent of Friedman, Hastie and Tibshirani (2008), in the implementation of Sustik and Calderhead (2012): the code is a C++ port of the Fortran routine of the glassoFast package, and returns the same result on ordinary input. It departs from it on degenerate input only:
it always terminates: non-finite input, or a coordinate with
\(S_{ii} + \rho_{ii} \leq 0\), is rejected (the result is filled with
NA and converged is FALSE), and the inner coordinate descent is
bounded, where glassoFast can loop forever on a nearly collapsed
covariance matrix;
failure to converge is reported through converged rather than silently;
it can be interrupted from R;
when S has no off-diagonal mass, the (diagonal) solution
\(1 / (S_{ii} + \rho_{ii})\) is returned, where glassoFast
returns \(1 / \max(\rho_{ii}, \epsilon)\);
it detects when it is cycling rather than converging, and stops.
That last point matters on an ill-conditioned S -- typically a
rank-deficient covariance, as arises when the number of variables approaches
the number of samples. The sweeps then settle into a small limit cycle: the
convergence criterion stops decreasing and oscillates just above its
threshold forever, while the solution itself no longer moves. glassoFast
spends its whole sweep budget on these and reports success regardless; here
the cycle is detected after stall_patience sweeps without progress, the
solve stops, and status reports "stalled".
Stopping early costs nothing, because sweeping on does not buy accuracy: the
iterates wander inside the cycle rather than settle. On a oaks-derived
covariance the solve stops after ~1100 sweeps instead of 10000, and both that
answer and the 10000-sweep one sit within 1e-3 of a 50000-sweep grind, with
the same number of edges and a couple of borderline ones differing --
the amplitude of the cycle, which no budget removes. The useful consequence
of "stalled" is thus not a warning about the solution, but the information
that the problem is in that regime: usually too weak a penalty for a
covariance that is (nearly) rank-deficient.
The remaining statuses are "max_iter" (maxit reached while still
progressing), "inner_failure" and "degenerate" (numerical trouble, the
result may contain NA).
J. Friedman, T. Hastie and R. Tibshirani (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3), 432--441.
M. A. Sustik and B. Calderhead (2012). GLASSOFAST: An efficient GLASSO implementation. UTCS Technical Report TR-12-29, The University of Texas at Austin.
data(trichoptera)
S <- cov(log1p(as.matrix(trichoptera$Abundance)))
fit <- graphical_lasso(S, rho = 0.1)
fit$converged
sum(fit$wi[upper.tri(fit$wi)] != 0) # number of edges
## penalty weights, leaving the diagonal unpenalized
W <- matrix(1, ncol(S), ncol(S)); diag(W) <- 0
fit <- graphical_lasso(S, rho = 0.1 * W)
Run the code above in your browser using DataLab