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_maxinterface for controlling the fits now use list generated by dedicated functions
The regression coefficients are now denoted by B, not Theta, such as B = t(Theta). We keep on sending back Theta as a field of myPLN$model_par$Theta, but this will soon be deprecated
PLNLDA() 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.