src/graphical_lasso.h is now the very same file as in PLNmodels, which
ported the same glassoFast Fortran: the two copies had drifted apart. It is
shared verbatim (namespace graphical_lasso), and brings to normalblockr:
ensure_pd() clamped its
eigenvalues and made it dense;graphical_lasso_fit() also returns status, delta and shift.The package now fits two complementary families. The original Normal-Block
model clusters variables by how they covary; the new mean-block family
clusters them by how they respond to the covariates, constraining the mean
as mu_i = C B' X_i with one regression profile per cluster.
normal_block(..., model = "mean") fits it, for a known clustering
(NormalBlockMeanKnownClusters) or an unknown one inferred by variational
EM (NormalBlockMeanUnknownClusters), over a range of cluster counts
(NormalBlockMeanCollectionClusters), a penalty path
(NormalBlockMeanCollectionSparsity) or both
(NormalBlockMeanCollectionClustersSparsity).ZINormalBlockMeanKnownClusters,
ZINormalBlockMeanUnknownClusters), and collections of those over a range
of cluster counts.noise_covariance gains a "full" shape for this family, whose Sigma is
the p x p residual covariance rather than a diagonal noise term.
"diagonal" is the default: a full Sigma spends p(p+1)/2 parameters that
drown the mean structure the criteria are weighing, and selects the number
of clusters markedly worse as p approaches n. Asking for sparsity > 0
implies "full", since a diagonal precision matrix has nothing for the
graphical lasso to penalize.normal_block_sequential() chains the two families: a mean-block fit, then
a variance-block fit on its residuals. A heuristic two-stage estimator, not
a joint model.NormalBlockData$new() gains a zeros argument, to carry an explicit
zero-inflation mask when the matrix handed to a model is no longer the one
carrying the zeros (the residuals of a first stage, say).This also fixes the segfault CRAN reported on
r-devel-linux-x86_64-fedora-gcc when re-building
breast-cancer-proteomics.Rmd (address 0x580, cause 'memory not mapped'),
which our own CI had been hitting intermittently on Linux.
glassoFast is no longer a dependency. The C++ (V)EM used to call back into
R once per M-step to reach it -- thousands of round trips from inside a
single .Call for one sparsity path, and a recurring source of intermittent
crashes. The solver now lives in src/graphical_lasso.h.glassoFast are fixed on the way. When the empirical covariance
has no off-diagonal mass the problem separates and the solution is diagonal;
it returned 1 / max(rho_ii, eps) there, dropping the variance term, so with
an unpenalized diagonal it returned about 9e15 instead of 1 / S_ii -- and
every q = 1 problem takes that branch. Separately, its inner loop could not
terminate on non-finite input.X0 with two or
more columns) failed with "non-conformable arguments": B0 was built
transposed. A second bug behind it made kappa non-finite when the response
had no zeros at all. Both were unreachable before, so no working fit changes.static).First CRAN submission
print(), summary(), plot(), logLik() and BIC() for fitted models (any NormalBlockVarBase subclass), and print()/summary()/logLik()/BIC() for collections of models; accessing $loglik on a collection now raises an informative error instead of silently returning NULL.@examples, cross-referenced $plot_network()/$latent_network() in each other's documentation, and replaced a few inefficient matrix operations (solve() on symmetric positive-definite matrices, M %*% t(C)) with chol2inv(chol())/tcrossprod().ZINormalBlockVarKnownClusters/ZINormalBlockVarUnknownClusters) for data with an excess of exact zeros.NormalBlockVarCollectionSparsity, sparsity = TRUE in normal_block()).ward2, kmeans, spectral, sbm, selectable via NB_control(clustering_init = )), and best_of_inits() to try several and keep the best-ELBO fit.refine() on NormalBlockVarCollectionClusters: post-hoc split/merge search seeded from neighboring cluster counts, to escape mediocre local optima left by independent per-q cold starts.brca_rppa (breast cancer proteomics), onema (French stream fish biomass, zero-inflated), university (WebKB text data), each with a dedicated vignette.NormalBlock* to NormalBlockVar*) for consistency; NormalBlockData is unchanged.NormalBlockData rescales columns of Y by default; fitted values and regression coefficients are reported back on the original scale.ClustOfVar dependency (the kmeansvar clustering heuristic was dropped after benchmarking showed it was both the worst-ranked and least reliable of the available heuristics).