Optim plnpca - explore block Newton with preconditionning and trust-region#171
Merged
Conversation
Profiles out the variational (M, S) block with a concave per-observation
Newton VE-step (the envelope theorem gives the reduced gradient of
g(B,C) = max_{M,S} ELBO for free), then optimises the loadings (B, C) with a
saddle-aware trust-region Newton: analytic Schur-complement Hessian-vector
products, a Jacobi preconditioner (analytic diagonal of L_theta_theta), and a
preconditioned Steihaug-CG with a P-norm trust region.
Unlike block-coordinate ascent, it escapes the saddle points of the
non-convex low-rank landscape, reaching a higher variational bound than the
nlopt CCSAQ backend on every tested dataset, while being faster on moderate p
(barents ~1.5x) and at parity on large p (oaks). Opt-in via
PLNPCA_param(backend = "trnewton"); the nlopt default is unchanged. Covariate
normalisation is skipped for this backend, whose trust region is not
scale-invariant. The 85 existing PLNPCA tests pass.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…olves Fixes failures found on large data (microcosm n=880, p=259), where trnewton crashed (solve(): solution not found) and never converged: * Scale-invariant stopping: the gradient norm scales with n and the counts, so the absolute gtol was never reached on large problems. Use a relative rule ||g|| <= gtol * max(||g0||, 1). * Robust linear solves for large counts: the per-observation M-Newton solve now uses the non-throwing arma::solve (dropping no_approx, falling back to a zero step), and the inner 2q x 2q block inverse uses a relative ridge with a pinv fallback instead of an absolute ridge that was negligible against large A. * Defaults retuned: maxit_out 60 -> 150 (the cap was cutting convergence short), gtol 1e-4 -> 1e-3 relative (sweet spot: best loglik at ~parity speed on small/moderate p; 1e-2 underconverges, 1e-4 is needlessly slow). trnewton now converges on microcosm (1.76x faster than nlopt, better bound in aggregate). The 85 PLNPCA tests pass; the nlopt default is unchanged. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Replace the previous PLNPCA "builtin" optimizer (spectral Barzilai-Borwein) with the profiled trust-region Newton, and expose it under the name "builtin" for consistency with the other PLN variants (where "builtin" is the home-made optimizer). No user-facing "trnewton" backend anymore; nlopt stays the default. * Delete the old spectral builtin rank optimizer (builtin_optimize_rank and builtin_optimize_vestep_rank). * Rename the trust-region Newton implementation to the "builtin" name: trnewton_plnpca.h -> builtin_plnpca.h (namespace trnewton -> builtin), wrappers_trnewton_optim_plnpca.cpp -> wrappers_builtin_optim_plnpca.cpp, trnewton_optimize_rank -> builtin_optimize_rank, backend tag "trnewton" -> "builtin". * R wiring: drop the "trnewton" backend choice; PLNPCA_param(backend = "builtin") now carries the trust-region config (cg_maxit, maxit_out, ftol_out, gtol, delta0); the builtin prediction VE-step falls back to the nlopt VE-step (the Newton backend has no dedicated VE-step); update the backend documentation. Rationale: across datasets the trust-region Newton reaches a higher variational bound than both the old builtin and nlopt, and it is not slower than the old builtin (which was in fact the slowest backend on large data). The 88 PLNPCA / backend tests pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Improves the outer trust-region convergence of the builtin (trust-region Newton) backend on large data, where it was crawling. A trace on microcosm showed the iterations were wasted on (i) an oversized initial radius (3 rejections before it settled) and (ii) a sawtooth radius with many rejected steps, each of which needlessly recomputed the gradient, the per-observation inner-block inverses and the preconditioner even though the iterate had not moved. * On a rejected step the iterate is unchanged, so the gradient / inner-block Hessian / preconditioner are now cached and refreshed only on acceptance; a rejection costs just one extra Steihaug-CG solve at the reduced radius. * Seed the initial radius 4x smaller and shrink to 0.25*||s|| (a fraction of the failed step) instead of 0.25*Delta. On microcosm (n=880, p=259) ranks that previously hit the iteration cap now converge (status 3), faster and at a higher variational bound. Small/moderate data is unchanged or slightly better (oaks +241 vs nlopt at ~parity speed). The 85 PLNPCA tests pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Adds, under trace > 1, the CG iteration count, exit reason (converged / negative curvature / trust-region boundary / iteration cap) and residual ratio for each outer trust-region step. This is what pinned down the large-n convergence behaviour: on microcosm the reduced Hessian is indefinite and the Steihaug-CG never fully converges (it hits the trust-region boundary or the iteration cap), so the outer loop is limited by the trust-region radius rather than by an under-solved Newton system. No change to the optimizer's behaviour. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… + DEVLOG * Note in the backend documentation that, on large data, the builtin (trust-region Newton) backend is bounded by maxit_out rather than by the gradient tolerance — the higher ranks keep improving the bound and trade quality for time; users wanting the best bound should raise maxit_out. * Add DEVLOG_2026-07-02-03.md summarising the session: SQUAREM evaluation (negative), the profiled trust-region Newton for PLNPCA, its promotion to the "builtin" backend, scale robustness, trust-region tuning, and the (understood, fundamental) large-n convergence limit. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Covers the previously untested builtin PLNPCA optimizer (the profiled trust-region Newton) and its backend-specific R path, fixing the codecov/patch failure: a family fit with backend = "builtin" checks it reaches a variational bound at least as good as nlopt, behaves like a normal PLNPCAfit, and honours the trust-region tuning keys (cg_maxit, maxit_out). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
A Quarto article (website-only, not built by R CMD check) documenting, for each model, which optimizer is used and its default, and detailing the mathematical principles of the home-made "builtin" backend: the joint (mean, variance) Newton VE-step shared by PLN / PLNnetwork / ZIPLN, the closed-form and graphical-lasso covariance M-steps, and the profiled saddle-aware trust-region Newton for PLNPCA (envelope-theorem gradient, matrix-free Schur reduced Hessian, Jacobi-preconditioned Steihaug-CG). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The singular "article/" (holding only the shared PLNreferences.bib) was easily confused with the plural "articles/" (the pkgdown convention for website-only articles, introduced for optimization.qmd). Rename to something unambiguous and update the bibliography: path in the 8 vignettes that reference it. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.