mnpca_mrl: Matrix-Normal PCA via Maximum Regularized Likelihood

In genpca: Generalized Principal Component Analysis

Matrix-Normal PCA via Maximum Regularized Likelihood

Description

Fits a rank-ncomp matrix factorization under matrix-normal noise with sparse row/column precision matrices:

Y = XW^\top + E,\quad E \sim \mathcal{MN}(0,\Omega,\Sigma)

using alternating block updates:

1. Alternating least squares updates for X, W with fixed precisions (Theta_row = Omega^{-1}, Theta_col = Sigma^{-1}).

2. Graphical-lasso style precision updates for Theta_row and Theta_col using ADMM, warm starts, and optional block screening.

The precision updates are solved on a residual-free scatter surrogate (rather than the exact residual E = Y - XW^T), so the two block updates do not jointly minimize a single shared objective at each step; as a result the reported objective_path is a convergence diagnostic, not a monotone objective trace (see Details and the objective_path return value).

Usage

mnpca_mrl(
  Y,
  ncomp = min(dim(Y)),
  lambda_row = 0.05,
  lambda_col = 0.05,
  max_outer = 25,
  max_inner = 5,
  tol = 1e-04,
  eps_ridge = 1e-08,
  jitter = 1e-06,
  center = TRUE,
  update_precisions = TRUE,
  warm_start = TRUE,
  gl_maxit = 200,
  gl_tol = 1e-04,
  gl_rho = 1,
  penalize_diagonal = FALSE,
  block_screen = TRUE,
  scale_fix = c("trace", "none"),
  sparsify_threshold = 1e-08,
  as_sparse_precision = TRUE,
  verbose = FALSE
)

Arguments

Y	Numeric matrix (n x p).
ncomp	Target rank (r).
lambda_row	L1 penalty for row precision (Theta_row).
lambda_col	L1 penalty for column precision (Theta_col).
max_outer	Maximum number of outer BCD iterations.
max_inner	Maximum ALS steps per outer iteration.
tol	Relative tolerance used for ALS and objective convergence checks.
eps_ridge	Ridge added to small r x r normal-equation systems.
jitter	Small diagonal jitter used in covariance/precision updates.
center	Logical; center columns of Y before fitting.
update_precisions	Logical; if FALSE, keeps identity precisions and runs weighted ALS only.
warm_start	Logical; warm start precision updates from previous iterate.
gl_maxit	Maximum ADMM iterations per graphical-lasso subproblem.
gl_tol	ADMM convergence tolerance for graphical-lasso subproblems.
gl_rho	ADMM augmented Lagrangian parameter.
penalize_diagonal	Logical; whether to penalize diagonal precision entries in L1 term. Default FALSE.
block_screen	Logical; use thresholded connected components to solve precision subproblems blockwise.
scale_fix	One of "trace" or "none". "trace" normalizes each precision to mean diagonal 1 after updates.
sparsify_threshold	Off-diagonal magnitude threshold used to set tiny precision entries to zero after each graphical-lasso solve.
as_sparse_precision	Logical; store precisions as sparse matrices when many entries are zero.
verbose	Logical; print iteration diagnostics.

Details

The covariance updates use low-rank correction identities and avoid explicit construction of E = Y - XW^T.

This function returns a plain S3 list of class "mnpca_mrl", not a multivarious bi_projector/cross_projector; the scores()/components()/reconstruct() generics used elsewhere in this package do not apply here. Use the returned $X, $W, and $fitted fields directly.

This implementation follows the maximum regularized likelihood (MRL) formulation of MN-PCA, combining low-rank factor updates with sparse precision estimation in row and column spaces, via alternating surrogate updates rather than joint minimization of a single objective at every step.

The main optimization target is:

\frac12\mathrm{tr}\left(\Theta_c (Y-XW^\top)^\top \Theta_r (Y-XW^\top)\right) -\frac{p}{2}\log|\Theta_r|-\frac{n}{2}\log|\Theta_c| + n\lambda_r\|\Theta_r\|_1 + p\lambda_c\|\Theta_c\|_1

with \Theta_r \succ 0, \Theta_c \succ 0. When update_precisions = FALSE, the method reduces to weighted low-rank approximation with fixed identity precisions.

Because the ALS factor updates and the graphical-lasso precision updates are solved against different surrogates of this target (the precision step uses a residual-free scatter matrix rather than the exact residual), the outer alternation is best understood as alternating surrogate updates with a relative-objective-change convergence heuristic, not classical block coordinate descent on a single monotone objective.

Value

An object of class "mnpca_mrl" with components:

X, W

Estimated low-rank factors (n x r, p x r).

Theta_row, Theta_col

Estimated row/column precision matrices.

fitted

Reconstructed matrix on input scale.

fitted_centered

Reconstructed centered matrix used in optimization.

residual_centered

Centered residual matrix.

objective_path

Objective value evaluated after each outer iteration's block updates, before any trace rescaling of the precisions. Because the factor updates (ALS) and precision updates (graphical lasso on a residual-free scatter surrogate) target different surrogates rather than a single shared objective, this path is a convergence heuristic and is not guaranteed to be monotone.

iterations

Number of outer iterations used.

converged

Logical; TRUE if the relative change in objective_path fell below tol before max_outer was reached. This is a convergence heuristic based on relative objective change, not a guarantee of a local optimum.

center

Column centering vector (or NULL).

call

Matched call.

References

Zhang, C., Gai, K., & Zhang, S. (2024). Matrix normal PCA for interpretable dimension reduction and graphical noise modeling. Pattern Recognition, 154, 110591. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.patcog.2024.110591")}

Friedman, J., Hastie, T., & Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3), 432-441.

Examples

set.seed(123)
n <- 20; p <- 12; r <- 3
Y <- matrix(rnorm(n * p), n, p)
fit <- mnpca_mrl(
  Y,
  ncomp = r,
  lambda_row = 0.08,
  lambda_col = 0.08,
  max_outer = 6,
  max_inner = 6,
  verbose = FALSE
)
dim(fit$X)
dim(fit$W)
length(fit$objective_path)

genpca documentation built on Sept. 17, 2026, 1:09 a.m.