Merge pull request #171 from PLN-team/optim-plnpca

Optim plnpca - explore block Newton with preconditionning and trust-region

Author: jchiquet

DESCRIPTION

@@ -54,6 +54,7 @@ Suggests:
     pheatmap,
     rmarkdown,
     spelling,
+    quarto,
     testthat
 LinkingTo:
     nloptr,

DEVLOG_2026-07-02-03.md

@@ -0,0 +1,140 @@
+# Journal de développement — 2-3 juillet 2026
+
+Branche : `optim-plnpca`
+
+## Contexte
+
+Exploration de méthodes d'accélération pour les procédures d'optimisation de PLNmodels,
+puis développement d'un nouvel optimiseur pour PLNPCA. Deux volets :
+1. Évaluation de SQUAREM (accélération de point-fixe EM) sur ZIPLN et PLNmixture.
+2. Conception, implémentation et intégration d'un **Newton profilé à région de confiance**
+   pour PLNPCA, qui remplace l'ancien backend `builtin`.
+
+---
+
+## 1. Évaluation SQUAREM (ZIPLN, PLNmixture) — négatif
+
+SQUAREM (Varadhan & Roland 2008) accélère une carte de point-fixe EM lente. Prototypes
+non-invasifs (scratchpad) reconstruisant les cartes EM de ZIPLN (`ZIPLNfit$optimize`) et
+PLNmixture (`PLNmixturefit$optimize`).
+
+**Verdict** : pas de gain de vitesse.
+- ZIPLN / PLNmixture : SQUAREM coupe ~1.5–2× les *cycles externes*, mais chaque cycle coûte
+  ~2–3 évaluations de la carte (VE-step / fits de composantes), et la garde de monotonie
+  ELBO rejette l'extrapolation → net plus lent (25–40 %). Léger mieux d'ELBO (échappe des
+  plateaux) mais pas de vitesse.
+- PLN / PLNPCA de base : hors-scope (solveur conjoint en C++, pas de carte EM séparable lente).
+
+Conclusion : SQUAREM = éventuelle option de robustesse, jamais un défaut. Le vrai levier pour
+PLNPCA est ailleurs (cf. §2).
+
+---
+
+## 2. PLNPCA — Newton profilé à région de confiance (nouveau backend)
+
+### 2.1 Diagnostic de départ
+
+À q fixé, l'optim PLNPCA (nlopt CCSAQ ou l'ancien builtin spectral BB) faisait des **milliers**
+d'itérations premier-ordre. L'ascension par blocs (VE-step Newton pour (M,S), puis (B,C))
+stagne sur des **points-selles** du paysage non-convexe des loadings ; c'est pourquoi les
+méthodes conjointes premier-ordre payent tant d'itérations (prix de la robustesse aux selles).
+
+### 2.2 L'idée
+
+Profiler les paramètres variationnels : `g(B,C) = max_{M,S} ELBO`.
+- Le bloc (M,S) est **concave** → VE-step Newton par observation (solve q×q), ~9 passes.
+- Le **théorème de l'enveloppe** donne `∇g` gratuitement (= ∂ELBO/∂(B,C) au (M,S) optimal).
+- Optimiser `g(B,C)` par un Newton **conscient des selles** (région de confiance + courbure
+  négative), en utilisant le **complément de Schur** `H_θθ − H_θφ H_φφ⁻¹ H_φθ` où `H_φφ` est
+  bloc-diagonal par observation (2q×2q) → produits Hessien-vecteur matrix-free bon marché.
+
+Validé en R d'abord : un pas le long de la courbure négative à la selle escalade vers un
+meilleur optimum que le package. Le Newton conscient des selles bat nlopt PARTOUT en ~10–40
+itérations externes (vs 463–5406 CCSAQ).
+
+### 2.3 Implémentation C++
+
+- `src/builtin_plnpca.h` : `ve_solve` (VE Newton), `inner_blocks_inv` (bloc interne analytique
+  2q×2q), `hess_dir` (Hessienne appliquée à une direction — dérivée analytiquement, `dA=A⊙dη`
+  + règle du produit, cohérente terme à terme avec `rank_obj_grad`), `precond_diag`
+  (préconditionneur de Jacobi = diagonale analytique de `L_θθ`, strictement positive).
+- `src/wrappers_builtin_optim_plnpca.cpp` : `builtin_optimize_rank` — gradient-enveloppe,
+  produit Hessien-vecteur de Schur (2 appels `hess_dir` + solve du bloc interne), **Steihaug-CG
+  préconditionné** (région de confiance en norme-P), boucle TR (test de ratio ρ, mise à jour
+  du rayon Δ).
+
+**Pièges rencontrés** :
+- `PLNPCAfit$optimize` normalise les covariables ; la région de confiance / les produits
+  Hessien-vecteur ne sont PAS invariants d'échelle → la normalisation DÉGRADE (oaks −1553 ll).
+  Fix : sauter la normalisation pour ce backend (échelle naturelle).
+- Compile : `C2` doit être déclaré AVANT les lambdas `[&]`.
+- Les produits Hessien-vecteur coûtent *intrinsèquement* 2 évaluations directionnelles de
+  gradient (une par direction θ/φ) — l'analytique gagne en propreté (pas d'ε FD), pas en vitesse.
+
+### 2.4 Préconditionnement
+
+Le Jacobi (diagonale de `L_θθ`) absorbe la disparité d'échelle B-vs-C → le Steihaug-CG converge
+en bien moins d'itérations (`cg_maxit` abaissé). Sur oaks (grand p), sans préconditionnement
+cg=6 donnait −2388 ll (inutilisable) ; avec, cg=8 → parité vitesse + meilleure qualité.
+
+### 2.5 builtin ← trnewton (renommage)
+
+Décision : garder torch (experimental), **supprimer l'ancien builtin spectral BB de PLNPCA** et
+exposer le Newton région-de-confiance sous le nom `builtin` (cohérent avec les autres variantes
+PLN où builtin = optimiseur maison). nlopt reste le défaut. Renommage complet
+(fichiers/namespace/fonction/tag), plus de backend "trnewton" user-facing.
+
+### 2.6 Comparaison des 4 backends
+
+Qualité (Δloglik vs meilleur de ligne, ranks 1:5) : **trnewton(builtin) > builtin-ancien > nlopt
+> torch** sur les 15 cas (trichoptera/barents/oaks). torch = pire partout ET le plus lent.
+L'ancien builtin spectral était même LE PLUS LENT à l'échelle (microcosm 351s vs nlopt 273s vs
+nouveau builtin 155s), sans niche → suppression justifiée.
+
+### 2.7 Robustesse à l'échelle
+
+Bug trouvé sur microcosm (n=880, p=259) : crash (`solve(): solution not found`) et non-convergence
+(gtol absolu jamais atteint car ‖g‖ scale en millions). Corrigés :
+- convergence **relative** `‖g‖ ≤ gtol·max(‖g₀‖,1)` ;
+- solves robustes (`arma::solve` non-levant + fallback ; ridge relatif + `pinv` sur le bloc) ;
+- défauts : `maxit_out 60→150`, `gtol 1e-4→1e-3` relatif (point idéal), `cg_maxit=8`.
+
+### 2.8 Gestion de la région de confiance
+
+Diagnostic (trace TR) : rayon initial surdimensionné (3 rejets) + dents de scie (30 rejets/120),
+chacun recalculant inutilement gradient/`inner_blocks_inv`/préconditionneur alors que le point ne
+bouge pas. Fix : **rejets bon marché** (quantités mises en cache, recalculées seulement à
+l'acceptation), seed Δ ÷4, rétrécissement `Δ=0.25·‖s‖`. Résultat : q1/q3 convergent maintenant
+(status 3), + rapides + meilleure qualité ; petit/moyen non régressé.
+
+### 2.9 Convergence lente à grand n — cause fondamentale (levier A épuisé)
+
+Instrumentation du Steihaug-CG (exit reason + résidu, dans le trace) : sur microcosm q3, le CG
+**ne converge JAMAIS** (0/120), il bute sur la frontière TR (83) ou le cap avec **résidu qui
+GROSSIT** (26). **Cause racine : le Hessien réduit est INDÉFINI à l'échelle** (paysage riche en
+selles) ET `g(B,C)` est très **non-quadratique aux grands pas** (le profilage de (M,S) est une
+carte non-linéaire violente) → pas limités par le rayon TR. Mais les pas restent **précis (ρ≈1)**
+et réduisent bien f : la rampe (~100 itér × ~2000 d'amélioration) est **fondamentale**.
+- Garde-fou "stopper CG si résidu grossit" → EMPIRE (les itérations errantes produisent quand
+  même de bons pas). Retiré.
+- Préconditionneur diagonal impuissant face à l'indéfini / directions plates non-alignées.
+- Seul levier théorique restant = solveur TR indéfini exact (GLTR/Lanczos), ~150 lignes, ROI
+  incertain (pas déjà exact, balle TR petite par non-quadraticité) → non retenu.
+
+### 2.10 Leviers pragmatiques — pas de gain gratuit
+
+- **Warm-start des rangs** (sequential=TRUE) : NET NÉGATIF pour le TR Newton (oaks 0.55× + pire).
+  Contrairement à un 1er-ordre, le coût par rang est dans la rampe (paysage indéfini), pas dans la
+  distance à l'init ; partir du rang q−1 biaise le bassin. Écarté.
+- **maxit_out** : pur cadran temps↔qualité (les hauts rangs rampent productivement : q5
+  +5.8k→+175k de maxit 150 à 350). Défaut 150 = compromis ; documenté que l'utilisateur peut
+  relever `maxit_out` pour la qualité à l'échelle (ajout doc `@param backend`).
+
+### 2.11 État final
+
+builtin (TR Newton) = meilleure qualité petit/moyen partout ; à l'échelle plus rapide (1.8×) et
+meilleur en agrégat, dominant sur tous les rangs si `maxit_out` relevé. La lenteur de convergence
+à grand n est comprise et fondamentale. 85 tests PLNPCA passent, nlopt reste défaut.
+
+**Commits** : `0880c681` (backend trnewton) · `98c6c35d` (robustesse échelle) · `a791a163`
+(trnewton→builtin) · `c9f38bec` (gestion TR) · `e641fe13` (diagnostics CG). Point d'arrêt naturel.

NEWS.md

@@ -2,9 +2,11 @@
 
 ## New backends and optimizers
 
-* **New built-in Newton optimizer** (`backend = "builtin"`) for PLN, ZIPLN and PLNPCA: 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.
+* **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.
 
-* **Backend defaults revisited package-wide**, based on extensive benchmarking: PLN keeps `"nlopt"` (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 four remain configurable via the `backend` argument; see the corresponding `*_param()` documentation for the trade-offs. The `torch` backend is now clearly marked **experimental** everywhere.
+* **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.
 

R/PLNPCA.R

@@ -59,11 +59,21 @@ PLNPCA <- function(formula, data, subset, weights, ranks = 1:5, control = PLNPCA
 #'
 #' @param backend optimization backend, either `"nlopt"` (default, NLOPT/CCSAQ, recommended
 #'   for PLNPCA: conservative per-variable steps reliably find the global basin even when
-#'   the singular-value ratio d1/sqrt(n) is large), `"builtin"` (joint L-BFGS with strong
-#'   Wolfe line search on all parameters simultaneously — faster per iteration but may
-#'   converge to inferior local optima on ill-conditioned datasets),
-#'   or `"torch"` (**experimental**: automatic differentiation via the torch package; tends to
-#'   find lower loglik than `"nlopt"` and `"builtin"` on most datasets — not recommended).
+#'   the singular-value ratio d1/sqrt(n) is large), `"builtin"` (the home-made optimizer, a
+#'   profiled trust-region Newton), or `"torch"` (**experimental**: automatic differentiation
+#'   via the torch package; tends to find lower loglik than `"nlopt"` and `"builtin"` on most
+#'   datasets — not recommended).
+#'   The `"builtin"` backend profiles out the variational parameters `(M, S)` with a
+#'   per-observation Newton VE-step and optimises the loadings `(B, C)` with a saddle-aware
+#'   trust-region Newton on the resulting objective (analytic Schur Hessian-vector products,
+#'   Jacobi-preconditioned Steihaug-CG). It reliably reaches a higher variational bound than
+#'   `"nlopt"` on small/moderate data at comparable speed, and is faster on large data;
+#'   tuning keys in `config_optim`: `cg_maxit`, `maxit_out`, `ftol_out`, `gtol`, `delta0`.
+#'   On large data (many samples and variables) the reduced-Hessian landscape is indefinite
+#'   and the outer trust-region iteration converges slowly but keeps improving the bound, so
+#'   it is capped by `maxit_out` (default 150) rather than by the gradient tolerance: the
+#'   higher ranks then trade quality for time. For the best variational bound on large data,
+#'   raise it, e.g. `config_optim = list(maxit_out = 300)`.
 #' @inheritParams PLN_param trace config_optim config_post
 #' @param init_method character: strategy used to compute the starting point for the shared SVD.
 #'   - `"LM"` (default): fast multivariate `lm.fit` on log-transformed counts. Good for
@@ -122,7 +132,11 @@ PLNPCA_param <- function(
   if (backend == "torch")
     message("torch backend is experimental: may converge to suboptimal solutions or fail on some datasets.")
   config_opt <- make_config_optim(backend, config_optim, trace,
-                                  builtin_default = config_default_plnpca)
+                                  builtin_default = config_default_plnpca,
+                                  extra = if (backend == "builtin")
+                                    list(maxit_out = 150L, ftol_out = 1e-8, cg_maxit = 8L,
+                                         gtol = 1e-3, delta0 = 1.0)
+                                  else list())
   config_opt$sequential <- sequential
 
   structure(list(

R/PLNPCAfit-class.R

@@ -253,7 +253,9 @@ PLNPCAfit <- R6Class(
         } else {
           private$optimizer$main <- nlopt_optimize_rank
         }
-        private$optimizer$vestep <- if (control$backend == "builtin") builtin_optimize_vestep_rank else nlopt_optimize_vestep_rank
+        ## The builtin (trust-region Newton) backend has no dedicated VE-step;
+        ## prediction reuses the nlopt VE-step for every backend.
+        private$optimizer$vestep <- nlopt_optimize_vestep_rank
         if (!is.null(control$svdM)) {
           svdM <- control$svdM
         } else {
@@ -303,7 +305,12 @@ PLNPCAfit <- R6Class(
       ## Optimization ----------------------
       #' @description Call to the C++ optimizer and update of the relevant fields
       optimize = function(responses, covariates, offsets, weights, config) {
-        nrm  <- normalize_covariates(covariates)
+        ## The builtin (trust-region Newton) backend is naturally well-conditioned
+        ## (analytic Hessian) and its trust region is not scale-invariant, so covariate
+        ## normalization degrades it: keep the natural scale for that backend.
+        nrm  <- if (identical(private$optimizer$main, builtin_optimize_rank))
+                  list(X_sc = covariates, scales = rep(1, ncol(covariates)))
+                else normalize_covariates(covariates)
         args <- list(data   = list(Y = responses, X = nrm$X_sc, O = offsets, w = weights),
                      params = list(B = sweep(private$B, 1, nrm$scales, "*"),
                                    C = private$C, M = private$M, S2 = private$S2),

R/RcppExports.R

@@ -81,10 +81,6 @@ builtin_optimize_rank <- function(data, params, config) {
     .Call('_PLNmodels_builtin_optimize_rank', PACKAGE = 'PLNmodels', data, params, config)
 }
 
-builtin_optimize_vestep_rank <- function(data, params, config) {
-    .Call('_PLNmodels_builtin_optimize_vestep_rank', PACKAGE = 'PLNmodels', data, params, config)
-}
-
 builtin_optimize_vestep_zipln_full <- function(data, params, config) {
     .Call('_PLNmodels_builtin_optimize_vestep_zipln_full', PACKAGE = 'PLNmodels', data, params, config)
 }

_pkgdown.yml

@@ -5,11 +5,11 @@ template:
 
 navbar:
   structure:
-    left:  [home, reference, slideshow, articles, news]
+    left:  [home, reference, vignettes, articles, news]
     right: [github]
   components:
     home: ~
-    articles:
+    vignettes:
       text: "Vignettes"
       menu:
       - text: "The trichoptera data set"
@@ -29,6 +29,13 @@ navbar:
       - text: "ZIPLN: Zero-inflated models for counts with excess zeros"
         href: articles/ZIPLN.html
 
+    articles:
+      text: "Articles"
+      menu:
+      - text: "Optimization in PLNmodels"
+        href: articles/optimization.html
+
+
 reference:
 - title: 'Top-level fitting functions'
   desc: >