PLNnetwork() and ZIPLNnetwork() now use an
internal graphical Lasso, a C++ port of the GLASSOFAST algorithm (Sustik
and Calderhead, 2012) shared with the normalblockr
package, instead of calling glassoFast::glassoFast().
glassoFast is no longer a dependency (it moves to
Suggests, for tests only).setTimeLimit()/R.utils::withTimeout(), which
then raise their usual error.PLNnetwork() and ZIPLNnetwork() fits with the
default "builtin" backend are unchanged (log-likelihoods
within 1e-11). Fits with backend = "nlopt" can differ
slightly along the path: that backend is sensitive to perturbations at
the level of the last floating-point digit. Speed is the same or
slightly better.graphical_lasso(S, rho, thr, maxit, w_init, wi_init),
returning w, wi, niter,
converged, status and delta, with
the same defaults as glassoFast and an optional warm start.PLNnetwork() produces when p ~ n) the sweeps settle into a
small limit cycle: the convergence criterion stops decreasing and
oscillates just above its threshold forever. glassoFast spends its whole
10000-sweep budget on these and reports success regardless. The cycle is
now detected after 1000 sweeps without progress (tunable through
stall_patience), and reported as
status = "stalled".$optim_par$glasso_nonconverged and
$optim_par$glasso_stalled, counted along the alternating
optimization). A warning is now raised only for the numerical failures
("degenerate", "inner_failure",
"max_iter"), not for a stalled solve, which says something
about the problem (too weak a penalty for a nearly rank-deficient
covariance) rather than about the solution.1 / (S_ii + rho_ii), where glassoFast returned
1 / max(rho_ii, 1.1e-16).PLNnetworkfit and
ZIPLNfit_sparse is now the one of Foygel and Drton
(2010), BIC - 2 gamma |E| log(p), which was
designed for the graphical Lasso, instead of the Stirling approximation
of Chen and Chen (2008)’s original criterion that was used so far. The
tuning parameter is exposed as $ebic_gamma (default
0.5, 0 gives back the BIC), on a single fit or
on a whole collection. EBIC values change, and so may the model selected
by getBestModel("EBIC").$density of a network is now
|E| / (p (p - 1) / 2): it divided the edge count by
p^2 rather than by the number of possible edges, and was
thus understated by a factor (p - 1) / p.igraph output of plot() on a
network fit, the opacity of an edge is now proportional to the strength
of its partial correlation, as its width already was. Dense networks are
no longer a solid blob of colour. The opacity of the weakest edge is the
new edge.alpha argument (default 0.2, set it
to 1 to restore uniformly opaque edges).PLNPCA()’s variational bound carried a
spurious p / 2 per observation, the entropy
constant of the full-covariance models: a rank-q model’s
variational distribution is over the q-dimensional scores,
and its q / 2 constant was already in the Kullback-Leibler
term. loglik, BIC and ICL of
every PLNPCAfit therefore decrease by
n * p / 2. Rank selection is unchanged (the term does not
depend on the rank), but a PLNPCA() fit is now comparable
with a PLN() one, or with any other model. Reported by
Nguyen Quang Huy (Actuarial Science Laboratory, National Economics
University, Vietnam).
PLNnetwork() and ZIPLNnetwork()
no longer fail outright when one model along the penalty path cannot be
fitted: the collection is truncated with a warning instead, and
stability_selection() treats a model with no estimated
network as fully unstable. A second line of defense on top of the solver
above (thanks @rfriedman22, #176).
ZIPLNnetwork() fitted its inception with the
regularization path’s truncated optimizer which is appropriate for
models that are warm-started from one another, but the inception starts
from scratch and was given the same budget. The inception is now fitted
with an untruncated optimizer, as PLNnetworkfamily.New built-in Newton optimizer
(backend = "builtin") for PLN, ZIPLN and PLNnetwork:
envelope-theorem Newton steps with strong Wolfe line search, no
dependency on NLOPT. Substantially faster and more accurate than nlopt
on large datasets with full covariance.
PLNPCA’s "builtin" backend rewritten as a
profiled trust-region Newton: the variational block
(M, S) is profiled out with a per-observation Newton
VE-step, and the loadings (B, C) are optimised with a
saddle-aware trust-region Newton on the resulting objective (analytic
Schur-complement Hessian-vector products, Jacobi-preconditioned
Steihaug-CG) — replacing the previous spectral projected-gradient
"builtin". It reliably reaches a higher variational bound
than "nlopt" at comparable-to-better speed; tuning keys
cg_maxit, maxit_out, ftol_out,
gtol, delta0 in config_optim (see
?PLNPCA_param). "nlopt" remains the default
for PLNPCA.
Backend defaults revisited package-wide, based
on extensive benchmarking: PLN and PLNPCA keep "nlopt" (PLN
now consistently faster thanks to profiled = TRUE, see
below); PLNnetwork and ZIPLNnetwork now default to
"builtin", which finds a better optimum at a modest speed
cost; ZIPLN keeps its "builtin" default. All backends
remain configurable via the backend argument; see the
corresponding *_param() documentation for the trade-offs.
The torch backend is now clearly marked
experimental everywhere.
Quality and speed improvements:
config_optim$profiled = TRUE is now the default for
full-covariance nlopt fits (faster, slightly better loglik); ZIPLN’s
variational step now optimises (M, ψ, R) jointly via Newton
instead of sequentially; PLNPCA shares a single SVD initialisation
across ranks and can warm-start from a pre-fitted PLNfit
for large ranks (inception/init_method, see
?PLNPCA_param); ZIPLN’s starting point no longer relies on
pscl::zeroinfl (now an internal LM + binomial GLM routine),
which is both much faster and a better starting point —
pscl is no longer a dependency.
Fixed a critical nlopt convergence bug affecting PLN/PLNPCA: ill-conditioned covariate scaling could trigger the XTOL stopping criterion after very few iterations, well before convergence. The built-in backend was never affected; nlopt is now also fixed via better parameter scaling.
CovTraitsBase in covariance_pln.h) instead of
being duplicated per structure. PLNPCA and ZIPLN’s variational step now
reuse the same machinery instead of separate hand-rolled
implementations, removing a substantial amount of duplicated code and
fixing minor inefficiencies along the way (e.g. ZIPLN’s VE-step used to
treat the precision matrix as dense even for diagonal/spherical
covariance, at unnecessary O(np^2) cost).{backend}_optimize_{structure} convention.Uniform covariate normalization: a
normalize_covariates() helper (zero mean, unit variance per
column) is now applied consistently in all optimize()
methods (PLN, PLNPCA, PLNnetwork, ZIPLN). This makes the nlopt XTOL
criterion scale-invariant and stabilises the torch backend.
Parallelism backend:
future.apply::future_lapply is replaced by
parallel::mclapply throughout (stability selection for
PLNnetwork / ZIPLNnetwork). Use options(mc.cores = N) to
set the number of cores.
Bug fixes: PLNnetwork/ZIPLNnetwork’s inception
(warm-start) model didn’t inherit
ftol_em/maxit_em from the user’s
config_optim, silently falling back to defaults and
producing a wrong penalty grid; the PLNPCA rank-model objective used
A − Y where it should use A − Y ⊙ Z; various
ZIPLN prediction/initialization fixes (#146, #149, #150, #152).
microcosm data now included (#153, #154); AIC added for PLN and ZIPLN classes (#151); other fixes (#155).
CRAN check fixes (no user-visible effect):
removed an unused -fopenmp compilation flag in
src/Makevars that was inadvertently turning on Armadillo’s
internal OpenMP parallelisation and inflating the CPU/elapsed time ratio
of several examples; fixed a GCC -Wmismatched-new-delete
false positive in src/packing.cpp’s internal test helper by
restructuring the code (no diagnostic-suppressing pragma
involved).
variance_jackknife() and
variance_bootstrap() to prevent estimation recycling,
results from those functions are now comparable to doing jackknife /
bootstrap “by hand”.predict() function for PLNfit model
to (i) return fitted values if newdata is missing or (ii) perform one VE
step to improve fit if responses are provided (fix issue #114)scale argument compute_offset() to force the
offsets (RLE, CSS, GMPR, Wrench) to be on the same scale as the counts,
like TSS.clusters) is not of the form 1:K_maxPLNLDA() and
changing extract_model() to conform with
model.frame()$VEStep() for PLN-PCA, dealing
with low rank matrices$project() for PLN-PCA, used to
project newdata into PCA space$latent_pos() which is equivalent to
active binding $latentNEWS.md file to track changes to the
package.