Woodbury Optimizations

[aschneuw]

• Improves Woodbury Identity optimizations
• properly takes into account pivoting for all Cholesky Decomposition based operations
• updates scikit-sparse to the latest version 0.5.0 (had breaking API changes)
• Introduces Ztilde_t <- phi.rphi %*% Zt as a further optimization

[aschneuw]

this was exactly my thinking :)

[fabian-s]

Review Summary

Thanks for the work on this, Arno — the getVb() refactoring and the Woodbury identity path are well-structured and the core linear algebra is correct. Numerical equivalence testing (7 scenarios: Gaussian/binomial/Poisson with various smooth types, random intercepts/slopes) confirms machine-precision agreement between the R-only path and master, and between the Woodbury and direct paths.

However, there are several issues that need to be addressed before this can be merged. See inline comments for details.

Must Fix

1. python_cholmod hardcoded to TRUE — will crash any user without Python/scikit-sparse
2. Permutation direction in phi.rphi — likely uses perm where inv_perm is needed
3. root.phi potentially undefined — when no non-smooth random effects exist

Should Fix

4. Comment typo introduced ("--get need")
5. Stale/misleading comment on root.phi line
6. scikit-sparse==0.5.0 exact pin is fragile
7. Dead V computation
8. Dead if (TRUE) {} else {} branches
9. Indentation inconsistency

[fabian-s]

Stale comment: (was ind,ind - seems wrong) — but the code IS using [ind,ind], same as master. The parenthetical should be removed to avoid confusion.

[fabian-s]

Sorry for the chaos - briefly lost control of my coding agent ... ;)
I fully endorse the comments it made on my behalf above.

[fabian-s]

I can’t push to aschneuw/r-gam-gamm4 from here (permission error on the fork), so I’m putting the proposed fix set here directly, see also https://github.com/r-gam/gamm4/tree/woodbury-fixes which you could simply merge into your PR.

I rechecked this locally with sksparse enabled via reticulate, and I think the remaining changes needed on this branch are:

1. Restore the Python Woodbury fallback when factoring phi fails.
2. Fix the direct python_cholmod path: the current manual permutation handling is not correct for scikit-sparse 0.5.0, and the current V_py = r_to_py(V) also hits a read-only buffer failure.
3. Reset the accidental executable bits on the touched package files.
4. Keep the small cleanup items (--get need typo, dead V <- Diagonal(...) line).

The direct-path fix is to use cho_factor(..., lower=TRUE) plus solve(system = "A"), because in sksparse 0.5.0 that is the solve mode that applies CHOLMOD’s permutation correctly.

Suggested code changes in R/gamm4.r:

if (woodbury && nrow(Zt)) {
  if (python_cholmod) {
    phi_py = reticulate::r_to_py(phi)$copy()
    cho_phi_py <- try(
      cholmod$cho_factor(phi_py, supernodal_mode="auto", lower=TRUE),
      silent = TRUE
    )
    if (inherits(cho_phi_py, "try-error")) {
      woodbury <- FALSE
    } else {
      perm <- as.integer(reticulate::py_to_r(cho_phi_py$get_perm()) + 1)
      inv_perm <- perm
      inv_perm[perm] <- 1:length(perm)
      phi.rphi <- reticulate::py_to_r(scipy_sparse$csc_matrix(cho_phi_py$R))[,inv_perm]
    }
  } else {
    phi.inv <- try(chol2inv(chol(phi)), silent = TRUE)
    if (inherits(phi.inv,"try-error")) woodbury <- FALSE
  }
}

if (python_cholmod) {
  # Only system="A" applies CHOLMOD's internal permutation correctly.
  V_py = reticulate::r_to_py(V)$copy()
  R <- cholmod$cho_factor(V_py, supernodal_mode="auto", lower=TRUE)
  XFP_py <- reticulate::r_to_py(Xfp)$copy()
  XF_py <- reticulate::r_to_py(Xf)$copy()
  XVXS <- t(Xfp) %*% reticulate::py_to_r(scipy_sparse$csc_matrix(R$solve(XFP_py, system="A"))) + Sp^2/scale
  XVX <- t(Xf) %*% reticulate::py_to_r(scipy_sparse$csc_matrix(R$solve(XF_py, system="A")))
} else {
  R <- mgcv::mchol(V);piv <- attr(R,"pivot")
  WX <- as(solve(t(R),Xfp[piv,]),"matrix")
  XVX <- as(solve(t(R),Xf[piv,]),"matrix")
  XVX <- crossprod(XVX)
  XVXS <- crossprod(WX)+Sp^2/scale
}

gam.terms <- attr(gmf,"terms") # terms object for `gam' part of fit -- need this for prediction to work properly

and remove the dead line:

V <- Diagonal(n=length(v),x=v)

I also recommend resetting the file modes on:

• .gitignore
• ChangeLog
• DESCRIPTION
• MD5
• NAMESPACE
• R/gamm4.r
• man/gamm4.Rd

Locally I validated the patched version with:

• python_cholmod=TRUE vs FALSE model-level smoke fits (Gaussian smooth-only, Gaussian + random intercept, Poisson + random intercept): all matched to numerical precision
• randomized direct-path getVb(..., woodbury=FALSE) comparisons against the native-R branch: no mismatches after the direct-path fix above
• a singular / semidefinite phi case: the Woodbury Python branch now falls back cleanly instead of hard-failing during factorization

From my side, once the above is applied, this looks ready to merge.

[aschneuw]

I can’t push to aschneuw/r-gam-gamm4 from here (permission error on the fork), so I’m putting the proposed fix set here directly, see also https://github.com/r-gam/gamm4/tree/woodbury-fixes which you could simply merge into your PR.

I rechecked this locally with sksparse enabled via reticulate, and I think the remaining changes needed on this branch are:

1. Restore the Python Woodbury fallback when factoring phi fails.
2. Fix the direct python_cholmod path: the current manual permutation handling is not correct for scikit-sparse 0.5.0, and the current V_py = r_to_py(V) also hits a read-only buffer failure.
3. Reset the accidental executable bits on the touched package files.
4. Keep the small cleanup items (--get need typo, dead V <- Diagonal(...) line).

The direct-path fix is to use cho_factor(..., lower=TRUE) plus solve(system = "A"), because in sksparse 0.5.0 that is the solve mode that applies CHOLMOD’s permutation correctly.

Suggested code changes in R/gamm4.r:

if (woodbury && nrow(Zt)) {
  if (python_cholmod) {
    phi_py = reticulate::r_to_py(phi)$copy()
    cho_phi_py <- try(
      cholmod$cho_factor(phi_py, supernodal_mode="auto", lower=TRUE),
      silent = TRUE
    )
    if (inherits(cho_phi_py, "try-error")) {
      woodbury <- FALSE
    } else {
      perm <- as.integer(reticulate::py_to_r(cho_phi_py$get_perm()) + 1)
      inv_perm <- perm
      inv_perm[perm] <- 1:length(perm)
      phi.rphi <- reticulate::py_to_r(scipy_sparse$csc_matrix(cho_phi_py$R))[,inv_perm]
    }
  } else {
    phi.inv <- try(chol2inv(chol(phi)), silent = TRUE)
    if (inherits(phi.inv,"try-error")) woodbury <- FALSE
  }
}

if (python_cholmod) {
  # Only system="A" applies CHOLMOD's internal permutation correctly.
  V_py = reticulate::r_to_py(V)$copy()
  R <- cholmod$cho_factor(V_py, supernodal_mode="auto", lower=TRUE)
  XFP_py <- reticulate::r_to_py(Xfp)$copy()
  XF_py <- reticulate::r_to_py(Xf)$copy()
  XVXS <- t(Xfp) %*% reticulate::py_to_r(scipy_sparse$csc_matrix(R$solve(XFP_py, system="A"))) + Sp^2/scale
  XVX <- t(Xf) %*% reticulate::py_to_r(scipy_sparse$csc_matrix(R$solve(XF_py, system="A")))
} else {
  R <- mgcv::mchol(V);piv <- attr(R,"pivot")
  WX <- as(solve(t(R),Xfp[piv,]),"matrix")
  XVX <- as(solve(t(R),Xf[piv,]),"matrix")
  XVX <- crossprod(XVX)
  XVXS <- crossprod(WX)+Sp^2/scale
}

gam.terms <- attr(gmf,"terms") # terms object for `gam' part of fit -- need this for prediction to work properly

and remove the dead line:

V <- Diagonal(n=length(v),x=v)

I also recommend resetting the file modes on:

• .gitignore
• ChangeLog
• DESCRIPTION
• MD5
• NAMESPACE
• R/gamm4.r
• man/gamm4.Rd

Locally I validated the patched version with:

• python_cholmod=TRUE vs FALSE model-level smoke fits (Gaussian smooth-only, Gaussian + random intercept, Poisson + random intercept): all matched to numerical precision
• randomized direct-path getVb(..., woodbury=FALSE) comparisons against the native-R branch: no mismatches after the direct-path fix above
• a singular / semidefinite phi case: the Woodbury Python branch now falls back cleanly instead of hard-failing during factorization

From my side, once the above is applied, this looks ready to merge.

Applied your suggested changes. Removed "# Only system="A" applies CHOLMOD's internal permutation correctly." Since this was misleading to me. I think the previous approach would have been doable as well with the appropriate permutations applied. But your approach works as well and is better.

[fabian-s]

• please resolve my earlier comment on line gamm4.r l. 587
• since the python_cholmod = TRUE branch needs py_require and that was only added in reticulate 1.41, please add the version requirement to DESCRIPTION

[fabian-s]

@simonnwood this is now ready to merge from my perspective.
please let us know if you require additional changes.