Package {genpca}

genpca: Generalized Principal Component Analysis

Description

Generalized PCA and related matrix decompositions in weighted inner-product spaces. Methods are based on Allen, G. I., Grosenick, L., and Taylor, J. (2014) doi:10.1080/01621459.2013.852978, "A generalized least-square matrix decomposition", Journal of the American Statistical Association, 109(505), 145-159; and Abdi, H. (2007), "Singular value decomposition (SVD) and generalized singular value decomposition" https://personal.utdallas.edu/~herve/Abdi-SVD2007-pretty.pdf, in "Encyclopedia of Measurement and Statistics", 907-912.

Author(s)

Maintainer: Brad Buchsbaum brad.buchsbaum@gmail.com [copyright holder]

See Also

Useful links:

• https://bbuchsbaum.github.io/genpca/

• https://github.com/bbuchsbaum/genpca

• Report bugs at https://github.com/bbuchsbaum/genpca/issues

Validate diagonal weights

Description

Weights must be finite and non-negative. Entries in ⁠[-rtol * max(abs(w)), 0)⁠ are set to exactly zero with a message; anything below is an error. Every backend sees the same cleaned vector.

Usage

.clamp_weights(w, rtol = .metric_rtol_default(), name = "A")

Evict least-recently-used cache entry

Description

Evict least-recently-used cache entry

Usage

.evict_lru()

Internal cache for GMD factorizations

Description

Internal cache for GMD factorizations

Usage

.gmd_cache

Format

An object of class environment of length 0.

Maximum number of cached entries (LRU eviction when exceeded)

Description

Maximum number of cached entries (LRU eviction when exceeded)

Usage

.gmd_cache_max_size

Format

An object of class integer of length 1.

Cache access times for LRU eviction

Description

Cache access times for LRU eviction

Usage

.gmd_cache_times

Format

An object of class environment of length 0.

Default relative tolerance for metric validation

Description

Default relative tolerance for metric validation

Usage

.metric_rtol_default()

Coerce to general CSC sparse matrix (dgCMatrix)

Description

Replacement for the deprecated direct as(., "dgCMatrix") coercion from symmetric/triangular/diagonal Matrix classes.

Usage

as_dgc(A)

Arguments

Value

a dgCMatrix

Coerce dense symmetric Matrix classes to general dense (dgeMatrix)

Description

Replacement for the deprecated direct as(., "dgeMatrix") coercion from dsyMatrix/dpoMatrix.

Usage

as_dge(A)

Arguments

Value

a dgeMatrix

Create Weight Operator Function

Description

Returns a closure that applies a weight matrix W or its transformations (square root, inverse, or combinations thereof) to vectors/matrices.

Usage

as_weight_operator(W, transpose = FALSE, sqrt = FALSE, inverse = FALSE)

Arguments

Value

A function that applies the requested transformation

Clip a symmetric matrix to the PSD cone

Description

Spectral clip: eigen-decompose and set negative eigenvalues to zero. Unlike ensure_spd() (a diagonal ridge shift), this preserves the non-negative part of the spectrum exactly. The output has no negative eigenvalue beyond reconstruction roundoff: the only fast path is an exact (unshifted) Cholesky success, which proves positive definiteness. Requires a dense eigendecomposition, so large sparse matrices are refused.

Usage

clip_psd(M, tol = NULL, dense_maxn = 2000L, name = "M")

Arguments

Value

a symmetric Matrix, PSD

Utilities for constraints

Description

Helpers to validate, symmetrize and (only when asked) repair constraint matrices. Two relative tolerances are used throughout the package: metric_rtol (default sqrt(.Machine$double.eps)) decides whether a metric is positive (semi)definite and which of its eigenvalues count as zero, and rank_rtol (see genpca()) decides which components are kept. Both are relative to the scale of the matrix, so every decision is invariant to rescaling the input.

Ensure SPD (sparse-friendly)

Description

Force a symmetric matrix to be symmetric positive definite: the result satisfies is_pd(., rtol = tol). Already-PD input is returned unchanged; otherwise a Gershgorin-based diagonal shift is applied, with a Matrix::nearPD() fallback for small dense matrices and an escalating jitter as a last resort.

Usage

ensure_spd(M, tol = 1e-06, nearpd_maxn = 2000L, name = "M")

Arguments

Value

a Matrix object (sparse stays sparse when possible)

Generalized eigenproblem on a covariance matrix

Description

Maximises v'Cv subject to v'Rv = 1 and the additional constraint that v lies in the retained range of R. Successive components are R-orthogonal. If P is the orthogonal projector onto that range, the returned vectors satisfy P C v = \lambda R v and V'RV = I. For a full-rank R this is the usual equation C v = \lambda R v. It also holds for a singular R when C maps its retained range into itself. Otherwise the component of C v outside the retained range need not vanish.

Usage

geigen_cov(
  C,
  R = NULL,
  ncomp = NULL,
  constraints_remedy = c("error", "ridge", "clip", "identity"),
  rank_rtol = 1e-06,
  metric_rtol = .metric_rtol_default(),
  verbose = FALSE
)

Arguments

Details

This is a different estimator from the GMD of genpca_cov, which uses R^{1/2} C R^{1/2}. If C and R commute, their common eigenvectors can be ordered differently: the GMD weights variances by metric eigenvalues, whereas this estimator divides by them. With R = c * I, the directions and their ordering agree, but the eigenvalue scales differ unless c = 1.

C is validated for symmetry but may be indefinite (the generalized eigenproblem is still defined); a warning is issued when its minimum eigenvalue is below -metric_rtol * scale. R must be positive semi-definite; an indefinite R is subject to constraints_remedy.

Value

A plain list with the same components as genpca_cov (v, d, lambda, k, propv, cumv, R_rank) and method = "geigen". propv is relative to \mathrm{tr}(R^{-1/2} C R^{-1/2}) on the range of R.

See Also

genpca_cov

Examples

C <- cov(scale(iris[,1:4], center=TRUE, scale=FALSE))
w <- c(1, 1, 0.5, 2)
fit_gmd <- genpca_cov(C, R = w, ncomp = 2)
fit_geigen <- geigen_cov(C, R = w, ncomp = 2)
# different estimators: the singular values generally differ
rbind(gmd = fit_gmd$d, geigen = fit_geigen$d)

# With singular R, the equation is projected onto its retained range
C <- matrix(c(2, 1, 1, 2), 2)
R <- diag(c(1, 0))
fit <- geigen_cov(C, R, ncomp = 1)
P <- diag(c(1, 0))
P %*% C %*% fit$v - (R %*% fit$v) * fit$lambda

Generalised Principal Components Analysis (GPCA)

Description

Implements the Generalised Least-Squares Matrix Decomposition of Allen, Grosenick & Taylor (2014) for data observed in a row inner-product space M and a column inner-product space A. Setting M = I_n, A = I_p recovers ordinary PCA.

Usage

genpca(
  X,
  A = NULL,
  M = NULL,
  ncomp = NULL,
  method = c("eigen", "auto", "spectra", "randomized", "deflation"),
  constraints_remedy = c("error", "ridge", "clip", "identity"),
  preproc = multivarious::pass(),
  threshold = 1e-06,
  maxit_deflation = 500L,
  use_cpp = TRUE,
  maxeig = 5000,
  warn_approx = TRUE,
  maxit_spectra = 1000,
  tol_spectra = 1e-09,
  rank_rtol = 1e-06,
  oversample = 20L,
  n_power = 1L,
  n_polish = 0L,
  jitter_metric = 1e-10,
  seed_randomized = 1234L,
  tol_polish_randomized = 1e-04,
  verbose = FALSE
)

Arguments

Value

An object of class c("genpca", "bi_projector") inheriting from multivarious::bi_projector, with slots including:

u,v

Left/right singular vectors scaled by the constraint metrics (MU, AV). These correspond to components in the original space's geometry. Use components(fit).

ou,ov

Orthonormal singular vectors in the constraint metric (U, V such that UT M U = I, VT AV = I). These are the core mathematical factors.

sdev

Generalised singular values d_k. Note these are singular values of the metric-whitened data matrix, not standard deviations: with identity metrics and centering, sdev = prcomp(X)$sdev * sqrt(nrow(X) - 1).

s

Scores: the generalised principal components ⁠z_k = X A ov_k = ou_k d_k⁠ (Allen et al. 2014, Section 2.4). Identical to project(fit, X) on the training data. Use scores(fit).

preproc

The multivarious pre-processing object used.

A, M

The constraint matrices used (potentially after coercion to sparse format).

propv

Proportion of generalized variance explained by each component.

cumv

Cumulative proportion of generalized variance explained.

Method

We compute the rank-ncomp factors UDVT that minimise

\|X - UDV^\top\|_{M,A}^2 = \mathrm{tr}\!\bigl(M\, (X-UDV^\top)\,A\,(X-UDV^\top)^\top\bigr)

subject to UT M U = I, VT AV = I. (Allen et al., 2014). Five methods are available via the method argument:

• "eigen" (Default): Uses a one-shot eigen decomposition strategy based on gmdLA. It explicitly forms and decomposes a p \times p or n \times n matrix (depending on n vs p).

• "auto": Chooses among "eigen", "spectra", and "randomized" using heuristics on shape, rank ratio (ncomp / min(n,p)), and constraint structure.

• "spectra": Computes the top-k singular triplets of the metric-whitened data F_M' X F_A (with M = F_M F_M', A = F_A F_A') as an implicit operator via the eigencore package, without forming the large intermediate matrix. Generally faster and uses less memory for large n or p when few components are requested.

• "randomized": Uses a randomized block range finder and small projected eigendecomposition. This is an approximate low-pass method that is often much faster for wide dense matrices with sparse metrics when only top components are needed.

• "deflation": Uses an iterative power/deflation algorithm. Can be slower but potentially uses less memory than "eigen" for very large dense problems where ncomp is small.

Backend Guidance

The default is method = "eigen"; "auto" is opt-in, not the default.

• Use "eigen" (the default) when you need a stable reference solution on small/medium problems.

• Use "auto" to let a heuristic pick among "eigen", "spectra", and "randomized" based on problem shape and constraint structure.

• Use "spectra" for larger matrix-free iterative solves where memory pressure is a concern.

• Use "randomized" for wide low-rank settings (p >> n) with sparse metrics when throughput matters most.

• Use "deflation" when you only need a few components and can tolerate iterative convergence behavior.

For pre-computed covariance matrices C = X'MX, see genpca_cov which performs GPCA directly on C with column constraint R (equivalent to A).

References

Allen, G. I., Grosenick, L., & Taylor, J. (2014). A Generalized Least-Squares Matrix Decomposition. Journal of the American Statistical Association, 109(505), 145-159. arXiv:1102.3074.

See Also

genpca_cov for GPCA on pre-computed covariance matrices, truncate.genpca, reconstruct.genpca, multivarious::bi_projector, multivarious::project, multivarious::scores, multivarious::components, multivarious::reconstruct.

Examples

if (requireNamespace("multivarious", quietly = TRUE)) {
  set.seed(123)
  X <- matrix(stats::rnorm(200 * 100), 200, 100)
  rownames(X) <- paste0("R", 1:200)
  colnames(X) <- paste0("C", 1:100)

  # Standard PCA (A=I, M=I, centered) - using default method="eigen"
  gpca_std_eigen <- genpca(X, ncomp = 5, preproc = multivarious::center(), verbose = FALSE)

  # Standard PCA using Spectra method (requires C++ build)
  # gpca_std_spectra <- try(genpca(X, ncomp = 5,
  #                              preproc = multivarious::center(),
  #                              method = "spectra", verbose = TRUE))
  # if (!inherits(gpca_std_spectra, "try-error")) {
  #    print(head(gpca_std_spectra$sdev))
  # }

  # Compare singular values with prcomp
  pr_std <- stats::prcomp(X, center = TRUE, scale. = FALSE)
  print("Eigen Method Sdev:")
  print(head(gpca_std_eigen$sdev))
  print("prcomp Sdev:")
  print(head(pr_std$sdev))
  print(paste("Total Var Explained (Eigen):",
              round(sum(gpca_std_eigen$propv) * 100), "%"))

  # Weighted column PCA (diagonal A, no centering)
  col_weights <- stats::runif(100, 0.5, 1.5)
  gpca_weighted <- genpca(X, A = col_weights, ncomp = 3,
                          preproc = multivarious::pass(), verbose = FALSE)
  print("Weighted GPCA Sdev:")
  print(gpca_weighted$sdev)
  print(head(components(gpca_weighted)))
}

Generalized PCA on a covariance matrix (GMD form)

Description

Performs Generalized PCA directly on a pre-computed covariance matrix ⁠C = X'MX⁠ with a single variable-side metric R, following Allen et al.'s GMD: the eigendecomposition of R^{1/2} C R^{1/2} mapped back with V = R^{-1/2} Z, so that V'RV = I. With ⁠C = X'MX⁠ and R = A this matches genpca(X, M = M, A = A) exactly. This is useful when you already have C or when X is too large to store but C is manageable.

Usage

genpca_cov(
  C,
  R = NULL,
  ncomp = NULL,
  method = c("gmd", "geigen"),
  constraints_remedy = c("error", "ridge", "clip", "identity"),
  rank_rtol = 1e-06,
  metric_rtol = .metric_rtol_default(),
  tol = NULL,
  verbose = FALSE
)

Arguments

Details

The generalized eigenproblem C v = \lambda R v is a different estimator (it maximises v'Cv subject to v'Rv = 1, which is generally gives different components from the GMD) and lives in its own function, geigen_cov. method = "geigen" is accepted here for one release and forwards to it with a deprecation warning.

Value

A plain list (not a multivarious bi_projector) with components:

v

p x k matrix of loadings (R-orthonormal eigenvectors)

d

Singular values (square root of eigenvalues lambda)

lambda

Eigenvalues (variances under the R-metric)

k

Number of components returned

propv

Proportion of variance explained by each component (total variance is \mathrm{tr}(CR), Allen et al. Corollary 5)

cumv

Cumulative proportion of variance explained

R_rank

Rank of the constraint matrix R

method

"gmd"

Because this is a plain list rather than a bi_projector, the multivarious generics scores(), components(), and reconstruct() do not apply to it; index $v/$d directly, or use genpca when you need the full projector interface on a data matrix rather than a pre-computed covariance matrix.

References

Allen, G. I., Grosenick, L., & Taylor, J. (2014). A Generalized Least-Squares Matrix Decomposition. Journal of the American Statistical Association, 109(505), 145-159.

See Also

geigen_cov for the generalized eigenproblem C v = \lambda R v, genpca for the two-sided GPCA on data matrices, genpls for generalized partial least squares

Examples

# Standard PCA on a covariance (no constraint)
C <- cov(scale(iris[,1:4], center=TRUE, scale=FALSE))
fit0 <- genpca_cov(C, R=NULL, ncomp=3)
print(fit0$d[1:3])       # first 3 singular values
print(fit0$propv[1:3])   # variance explained by first 3 components

# Equivalence with genpca()
set.seed(123)
X <- matrix(rnorm(50 * 10), 50, 10)
M_diag <- runif(50, 0.5, 1.5)  # row weights
A_diag <- runif(10, 0.5, 2)    # column weights
fit_gpca <- genpca(X, M = M_diag, A = A_diag, ncomp = 5,
                   preproc = multivarious::pass())
C <- crossprod(X, diag(M_diag) %*% X)  # C = X'MX
fit_cov <- genpca_cov(C, R = A_diag, ncomp = 5)
all.equal(fit_gpca$sdev, fit_cov$d, tolerance = 1e-10)

# Variable weights via a diagonal metric (iris covariance, 4 variables)
C_iris <- cov(scale(iris[,1:4], center=TRUE, scale=FALSE))
w <- c(1, 1, 0.5, 2)
fitW <- genpca_cov(C_iris, R = w, ncomp=3)
print(fitW$d[1:3])

Generalized eigenvalue-based covariance GPCA (internal)

Description

Solves the generalized eigenproblem projected onto the retained range of R. This is the original implementation that was in gpca.R.

Usage

genpca_cov_geigen(
  C,
  R = NULL,
  ncomp = NULL,
  constraints_remedy = c("error", "ridge", "clip", "identity"),
  rank_rtol = 1e-06,
  metric_rtol = .metric_rtol_default(),
  verbose = FALSE
)

GMD-based covariance GPCA (internal)

Description

Implements Allen et al.'s GMD approach for covariance matrices. Computes eigendecomposition of R^{1/2} C R^{1/2} and maps back.

Usage

genpca_cov_gmd(
  C,
  R = NULL,
  ncomp = NULL,
  rank_rtol = 1e-06,
  metric_rtol = .metric_rtol_default(),
  verbose = FALSE
)

Generalized PLS via Implicit Operator (PLS-SVD / GPLSSVD)

Description

Canonical (two-block) generalized PLS using sparse-friendly implicit matrix-vector products. Solves the SVD of the operator S = Xe' Ye without materializing Xe = Mx^{1/2} X Ax^{1/2} or Ye = My^{1/2} Y Ay^{1/2}.

Usage

genpls(
  X,
  Y,
  Ax = NULL,
  Ay = NULL,
  Mx = NULL,
  My = NULL,
  ncomp = 2,
  preproc_x = multivarious::pass(),
  preproc_y = multivarious::pass(),
  svd_backend = c("eigencore", "irlba", "RSpectra"),
  svd_opts = list(tol = 1e-07, maxitr = 1000),
  constraints_remedy = c("error", "ridge", "clip", "identity"),
  verbose = FALSE
)

Arguments

Details

This follows the GPLSSVD/PLS-SVD formulation (Beaton, eqs. 10–14): the top ncomp singular triplets of S are computed by iterative SVD on the linear maps v -> S v and u -> S^T u, implemented with metric Cholesky multiplies/solves when possible. Works with dense or sparse Matrix inputs and constraint metrics.

Returns a multivarious::cross_projector with X-/Y-weights (vx, vy) chosen to provide natural projection of new data (X %*% vx, Y %*% vy). Additional GPLSSVD quantities are attached to the object for access: singular values d, generalized weights p, q, variable scores fi, fj, and row latent variables lx, ly.

genpls() maximizes the covariance between latent variables of X and Y under the (Mx, Ax, My, Ay) metrics by computing the SVD of S = Xe' Ye, where Xe = Mx^{1/2} X Ax^{1/2} and Ye = My^{1/2} Y Ay^{1/2}.

project() on a fitted object returns latent variables in the ambient (original data) metric, i.e. ⁠X \%*\% vx = X W_X p⁠. This differs from the attached ⁠lx = Mx^{1/2} X W_X p⁠, which lives in the row-whitened metric, by the factor Mx^{1/2}: lx and project(fit, X) are equal only when Mx = I. New-row projection necessarily uses project()'s ambient-metric convention, because a training-row metric Mx has no natural extension to out-of-sample rows.

Metric naming. genpls()'s row/column metric arguments (Mx, My for rows; Ax, Ay for columns) follow the same M/A convention as genpca(). Internally they are forwarded to gplssvd_op(), which uses the names XLW/YLW (left/row weights, i.e. Mx/My) and XRW/YRW (right/column weights, i.e. Ax/Ay).

Value

An object of class c("genpls", "cross_projector", "projector") with:

vx, vy

X- and Y- projection weights (stored in cross_projector) such that ⁠project(fit, X) = X \%*\% vx⁠ and ⁠project(fit, Y, source = "Y") = Y \%*\% vy⁠ recover the latent variables in the ambient (non-whitened) metric. Algebraically ⁠vx = W_X p = fi \%*\% diag(1/d)⁠ and ⁠vy = W_Y q = fj \%*\% diag(1/d)⁠ (see Details).

d

singular values of S = Xe' Ye (attached field)

p, q

generalized weights W_X^{-1/2} u, W_Y^{-1/2} v (attached)

fi, fj

variable/component scores W_X p D, W_Y q D (attached)

lx, ly

row latent variables M_X^{1/2} X W_X p, M_Y^{1/2} Y W_Y q (attached)

metrics

the supplied metrics (attached)

ncomp

Number of components actually extracted. The underlying operator may return fewer than the requested ncomp (e.g. when ncomp exceeds min(ncol(X), ncol(Y))); this field reflects the actual count, not the request.

backend

The svd_backend value passed in (for reference only; see the svd_backend argument for when it is actually used).

preproc_x, preproc_y

The fitted multivarious preprocessing objects for X and Y, stored on the cross_projector.

References

Beaton, D. (2020). Generalized eigen, singular value, and partial least squares decompositions: The GSVD package. (Eqs. 10-14). arXiv:2010.14734.

Examples

if (requireNamespace("multivarious", quietly = TRUE)) {
  set.seed(1)
  n <- 100; p <- 40; q <- 30
  X <- matrix(rnorm(n*p), n, p)
  Y <- matrix(rnorm(n*q), n, q)
  w <- runif(n); w <- w/sum(w)
  Mx <- My <- Matrix::Diagonal(x = w)
  fit <- genpls(X, Y, Mx = Mx, My = My, ncomp = 2,
                preproc_x = multivarious::center(),
                preproc_y = multivarious::center())
  fit$d  # singular values
}

Canonical Generalized PLS (alias)

Description

Convenience alias for genpls(); computes canonical generalized PLS (PLS-SVD/GPLSSVD). See ?genpls for full documentation.

Usage

genplsc(
  X,
  Y,
  Ax = NULL,
  Ay = NULL,
  Mx = NULL,
  My = NULL,
  ncomp = 2,
  preproc_x = multivarious::pass(),
  preproc_y = multivarious::pass(),
  svd_backend = c("eigencore", "irlba", "RSpectra"),
  svd_opts = list(tol = 1e-07, maxitr = 1000),
  constraints_remedy = c("error", "ridge", "clip", "identity"),
  verbose = FALSE
)

Arguments

Value

An object of class c("genpls", "cross_projector", "projector") with the same structure as genpls() returns (X-/Y-weights vx/vy, singular values d, generalized weights p/q, scores fi/fj, latent variables lx/ly, ncomp, and backend); see ?genpls for the definition of each slot.

References

Beaton, D. (2020). Generalized eigen, singular value, and partial least squares decompositions: The GSVD package. (Eqs. 10-14). arXiv:2010.14734.

See Also

genpls()

Examples

set.seed(1)
X <- matrix(rnorm(60 * 5), 60, 5)
Y <- matrix(rnorm(60 * 4), 60, 4)
fit <- genplsc(X, Y, ncomp = 2,
               preproc_x = multivarious::center(),
               preproc_y = multivarious::center())
fit$d

Get (and cache) a lower Cholesky factor for a dense SPD matrix

Description

Get (and cache) a lower Cholesky factor for a dense SPD matrix

Usage

get_chol_lower_dense(A)

Arguments

Value

a base numeric matrix L (lower triangular) with A = L %*% t(L)

Clear internal cache for matrix decompositions

Description

Clears the internal cache used by generalized matrix decomposition functions. This can be useful to free up memory or when working with different datasets.

Usage

gmd_clear_cache()

Value

Invisibly returns TRUE after clearing the cache.

Examples

# Clear the internal cache
gmd_clear_cache()

Generalized matrix decomposition via partial SVD of the whitened operator

Description

Computes the generalized SVD of X with row metric Q and column metric R, equivalent to the eigendecomposition used by genpca with method = "eigen". The metrics are factored once (⁠Q = F_Q F_Q'⁠, ⁠R = F_R F_R'⁠; diagonal, dense or sparse Cholesky, or an eigen factor for singular metrics) and the top-k singular triplets of the implicit operator ⁠F_Q' X F_R⁠ are computed with eigencore; a dense SVD is used when few components are not requested or the iterative solver does not converge. gmd_fast_cpp() is an alias kept for existing callers.

Usage

gmd_spectra(
  X,
  Q,
  R,
  k,
  tol = 1e-09,
  maxit = 1000L,
  seed = 1234L,
  topk = TRUE,
  cache = TRUE,
  auto_topk = TRUE,
  topk_ratio = 0.08,
  topk_min_dim = 200L,
  diag_fast = TRUE,
  rank_rtol = 1e-06,
  metric_rtol = .metric_rtol_default(),
  dense_maxn = 5000L
)

gmd_fast_cpp(
  X,
  Q,
  R,
  k,
  tol = 1e-09,
  maxit = 1000L,
  seed = 1234L,
  topk = TRUE,
  cache = TRUE,
  auto_topk = TRUE,
  topk_ratio = 0.08,
  topk_min_dim = 200L,
  diag_fast = TRUE,
  rank_rtol = 1e-06,
  metric_rtol = .metric_rtol_default(),
  dense_maxn = 5000L
)

Arguments

Value

A list with components:

u

n x k matrix of metric-weighted scores Q ou D

v

p x k matrix of components R ov

ou,ov

metric-orthonormal factors

d

length-k vector of singular values

k

number of components returned (may be < requested if rank-deficient)

When is this fast

• k << min(n, p): only the top-k triplets are computed

• Repeated calls with the same dense Q or R: Cholesky factors are cached

• Sparse metrics: only sparse factors and products are formed

A positive definite metric on the big side of X costs one Cholesky of that dimension; a singular one is never factored (symmetric small-side form).

See Also

genpca for the high-level interface, gmd_clear_cache to clear the Cholesky cache

Experimental penalized-ML estimation of GPCA metrics

Description

Alternates between GPCA factor estimation and penalized maximum-likelihood updates of the row/column metric matrices (M, A) under a Gaussian matrix-normal error model with a low-rank mean. Each iteration performs three exact block minimizations of a single penalized objective:

1. Fit GPCA with current A, M: the GMD theorem makes this the best rank-ncomp fit in the (M, A) norm, so it exactly minimizes the residual term.

2. Update \Sigma_r = E A E^T / p + \lambda I (the exact block minimizer under the ridge penalty), set M = solve(Sigma_r).

3. Update \Sigma_c = E^T M E / n + \lambda I using the updated M (sequential flip-flop, Dutilleul 1999), set A = solve(Sigma_c).

Usage

gpca_mle(
  X,
  ncomp = min(dim(X)),
  max_iter = 20,
  lambda = 0.001,
  scale_fix = c("none", "trace", "det"),
  tol = 1e-04,
  method = "eigen",
  constraints_remedy = "error",
  preproc = multivarious::pass(),
  verbose = FALSE,
  ...
)

Arguments

Details

The objective is the matrix-normal log-likelihood with a low-rank mean and an inverse-Wishart-style ridge penalty \lambda\,(p\,\mathrm{tr}\,\Sigma_r^{-1} + n\,\mathrm{tr}\,\Sigma_c^{-1}) (a MAP estimate). Because every block update is an exact minimizer of this one objective, loglik_path is monotone non-decreasing up to numerical noise. The penalty also resolves the c\,\Sigma_r, \Sigma_c/c scale indeterminacy, so the converged metrics are the penalized optimum and no rescaling is needed (scale_fix = "none", the default). scale_fix = "trace" or "det" additionally applies a joint reciprocal rescale at exit (row covariance normalized, factor absorbed into the column covariance). The unpenalized matrix-normal likelihood is invariant to that rescale, but the penalty p\lambda\,\mathrm{tr}(M) + n\lambda\,\mathrm{tr}(A) is not, so the rescaled metrics are no longer the penalized optimum; the returned loglik is always evaluated at the returned metrics and loglik_rescale_delta reports how far the rescale moved it. The algorithm stops when the relative change in the penalized log-likelihood falls below tol or max_iter is reached. Increase lambda or reduce ncomp if iterations become unstable. The objective is not identifiable with lambda = 0.

Value

A list with elements fit (a genpca fit computed with the returned metrics), A, M (learned SPD metrics), loglik (the penalized log-likelihood evaluated at the returned M, A and fit), loglik_unpenalized (the same without the lambda penalty), loglik_rescale_delta (the change in the penalty contribution caused solely by reciprocal metric rescaling; exactly zero for scale_fix = "none"), loglik_refit_delta (the remaining change from the last path value to loglik, including final refitting and numerical objective reevaluation), and loglik_path (the penalized log-likelihood after each outer iteration; monotone non-decreasing up to numerical noise, since every block update exactly minimizes the shared penalized objective). Values omit additive constants and include the lambda penalty, so they are comparable across iterations and across runs with the same lambda, but not across different lambda values.

References

Dutilleul, P. (1999). The MLE algorithm for the matrix normal distribution. Journal of Statistical Computation and Simulation, 64(2), 105-123.

Examples

if (requireNamespace("multivarious", quietly = TRUE)) {
  set.seed(123)
  X <- matrix(rnorm(40), 8, 5)
  res <- gpca_mle(X, ncomp = 2, max_iter = 5, lambda = 1e-3,
                  scale_fix = "trace", verbose = FALSE)
  # Learned metrics are SPD and match dimensions
  dim(res$A); dim(res$M)
  res$loglik_path
}

Generalized PLS-SVD via Implicit Operator (memory-safe)

Description

Compute the top-k singular triplets of S = Xe' Ye without materializing the whitened matrices Xe = Mx^{1/2} X Wx^{1/2}, Ye = My^{1/2} Y Wy^{1/2} when doing so would densify sparse data. When the whitening is sparsity-preserving (identity/diagonal metrics) or the data are dense, the whitened blocks are precomputed once so each matrix-vector product in the iterative SVD costs two multiplies.

Usage

gplssvd_op(
  X,
  Y,
  XLW = NULL,
  YLW = NULL,
  XRW = NULL,
  YRW = NULL,
  k = 2,
  center = FALSE,
  scale = FALSE,
  svd_backend = c("eigencore", "irlba", "RSpectra"),
  svd_opts = list(tol = 1e-07, maxitr = 1000),
  constraints_remedy = c("error", "ridge", "clip", "identity")
)

Arguments

Details

Naming map: this function names its metrics XLW/YLW (left/row weights) and XRW/YRW (right/column weights); these correspond to Mx/My (row metrics) and Ax/Ay (column metrics) in genpca()'s and genpls()'s M/A convention.

Value

A list with elements:

d

Length-k numeric vector of singular values of S = Xe' Ye.

u

⁠I x k⁠ matrix; left singular vectors of S (orthonormal in the Euclidean metric).

v

⁠J x k⁠ matrix; right singular vectors of S (orthonormal in the Euclidean metric).

p

⁠I x k⁠ matrix of generalized X-weights, p = W_X^{-1/2} u.

q

⁠J x k⁠ matrix of generalized Y-weights, q = W_Y^{-1/2} v.

fi

⁠I x k⁠ matrix of X-variable scores, F_i = W_X p D (columns of p rescaled by the singular values).

fj

⁠J x k⁠ matrix of Y-variable scores, F_j = W_Y q D.

lx

⁠N x k⁠ matrix of X row latent variables, L_x = M_X^{1/2} X W_X p.

ly

⁠N x k⁠ matrix of Y row latent variables, L_y = M_Y^{1/2} Y W_Y q.

k

Integer; number of components actually returned (may be less than the requested k if it exceeded min(I, J)).

dims

A list list(N, I, J) with the row count N and column counts I = ncol(X), J = ncol(Y).

center

A list list(X, Y) of the length-I / length-J column means subtracted from X/Y (all zero when center = FALSE).

scale

A list list(X, Y) of the length-I / length-J column scale factors divided out of X/Y (all one when scale = FALSE).

References

Beaton, D. (2020). Generalized eigen, singular value, and partial least squares decompositions: The GSVD package. arXiv:2010.14734.

Abdi, H. (2007). Partial least square regression PLS-Regression. In N. Salkind (Ed.), Encyclopedia of Measurement and Statistics. Thousand Oaks, CA: Sage.

Examples

set.seed(1)
X <- matrix(rnorm(40 * 6), 40, 6)
Y <- matrix(rnorm(40 * 4), 40, 4)
op <- gplssvd_op(X, Y, k = 2, center = TRUE)
round(op$d, 3)

Check whether a metric matrix is diagonal

Description

Check whether a metric matrix is diagonal

Usage

is_diagonal_metric(M)

Arguments

Value

logical

Check if constraint matrix is identity or purely diagonal with all != 1

Description

Check if constraint matrix is identity or purely diagonal with all != 1

Usage

is_identity_or_diag(M, eps = 1e-15)

Arguments

Value

TRUE if M is a (Matrix-based) diagonal with all or partial diag

Test positive semi-definiteness (relative tolerance)

Description

is_psd() is TRUE when A is symmetric and every eigenvalue exceeds -rtol * max(abs(diag(A))); is_pd() is TRUE when every eigenvalue exceeds +rtol * max(abs(diag(A))). Both are shifted Cholesky probes, so large sparse matrices never need an eigendecomposition. is_spd() is a deprecated alias of is_psd() kept for internal callers (its tol is the relative tolerance).

Usage

is_psd(A, rtol = .metric_rtol_default())

is_pd(A, rtol = .metric_rtol_default())

is_spd(A, tol = .metric_rtol_default())

Arguments

Value

logical

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

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. 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)

Internal Utilities for Partial Eigen, Adaptive Rank, and Sqrt Transforms

Description

These functions are used by both genpls and genplscorr to handle partial eigen expansions, diagonal/identity shortcuts, adaptive rank selection, and row/column transformations for data embeddings.

Prepare and validate constraint matrices

Description

Coerces A/M (NULL, weight vector, diagonal, dense or sparse matrix) to Matrix objects and validates them: finite entries, symmetry within roundoff (see symmetrize_or_stop()), and positive semi-definiteness within tol relative to the scale of the matrix. The requested remedy is applied to a metric that fails the PSD check (or has negative eigenvalues within tolerance when explicit clipping is requested), and every repair emits a warning of class genpca_metric_repaired carrying the repair_metric() report; valid PSD metrics, singular ones included, pass through under every remedy. Explicit clipping removes even negative eigenvalues within the validation tolerance. Asymmetric input is an error under every remedy.

Usage

prep_constraints(
  X,
  A,
  M,
  tol = .metric_rtol_default(),
  remedy = c("error", "ridge", "clip", "identity"),
  verbose = FALSE
)

Arguments

Value

list with elements A and M

Print an sfpca fit

Description

Prints a one-line summary of an sfpca object: number of components, dimensions, and singular values.

Usage

## S3 method for class 'sfpca'
print(x, ...)

Arguments

Value

x, invisibly.

Reconstruct data from a genpca fit

Description

Reconstructs (an approximation of) the original data from a genpca fit as ou[, comp] %*% diag(d[comp]) %*% t(ov[, comp]), followed by the inverse of the preprocessing transform. With all components and full rank this recovers the original data.

Usage

## S3 method for class 'genpca'
reconstruct(
  x,
  comp = 1:multivarious::ncomp(x),
  rowind = NULL,
  colind = NULL,
  ...
)

Arguments

Value

A numeric matrix of dimension ⁠length(rowind) x length(colind)⁠.

See Also

genpca(), truncate.genpca()

Examples

X <- matrix(rnorm(60), 15, 4)
fit <- genpca(X, ncomp = 4, preproc = multivarious::center())
max(abs(reconstruct(fit) - X)) # ~ 0 at full rank

Reconstruct data from an sfpca fit

Description

Reconstructs the rank-K sfpca model as U D V', using the stored (non-orthogonal) factors directly.

Usage

## S3 method for class 'sfpca'
reconstruct(
  x,
  comp = seq_len(multivarious::ncomp(x)),
  rowind = NULL,
  colind = NULL,
  ...
)

Arguments

Details

sfpca components are Euclidean unit vectors but are not mutually orthogonal, so V'V \ne I. The inherited reconstruct.bi_projector() method reconstructs through the Moore-Penrose pseudoinverse of the loadings (⁠scores \%*\% pinv(V)⁠), which for non-orthogonal V does not return the rank-comp model U D V' that sfpca() actually fits and deflates with. This method instead computes ⁠scores(x)[rowind, comp] \%*\% t(components(x)[colind, comp])⁠, i.e. U D V' restricted to the requested rows/columns/components. sfpca() does no preprocessing, so no inverse transform is applied.

Value

A numeric matrix of dimension ⁠length(rowind) x length(colind)⁠, the rank-length(comp) reconstruction U D V' using sfpca's stored non-orthogonal factors.

See Also

sfpca()

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

multivarious

components, reconstruct, transfer

Repair a metric matrix explicitly

Description

Returns a positive (semi)definite version of A together with a diagnostic report of what was done. This is the explicit counterpart of the constraints_remedy argument of genpca(): nothing in the package repairs a metric silently, and this function lets you inspect the repair before using the result.

Usage

repair_metric(
  A,
  method = c("ridge", "clip", "identity"),
  rtol = .metric_rtol_default(),
  name = "A",
  diag_maxn = 2000L
)

Arguments

Value

The repaired matrix (a Matrix), with attribute "repair_report" of class "metric_repair_report": a list with name, method, changed, n, min_eigenvalue_before, min_eigenvalue_after, gershgorin_bound_before, shift (diagonal loading added by "ridge", NA for "clip"), rank, condition_number and rtol.

See Also

genpca() (argument constraints_remedy)

Examples

A <- matrix(c(1, 2, 2, 1), 2)           # eigenvalues 3 and -1
B <- repair_metric(A, method = "ridge")
attr(B, "repair_report")
C <- repair_metric(A, method = "clip")
eigen(as.matrix(C))$values

Regularised / Generalised Partial Least Squares (RPLS / GPLS)

Description

Implements the algorithm of Allen et al. (2013) for supervised dimension-reduction with optional sparsity (\ell_1) or ridge (\ell_2) penalties and the generalised extension that operates in a user-supplied quadratic form Q.

Usage

rpls(
  X,
  Y,
  K = 2,
  lambda = 0.1,
  penalty = c("l1", "ridge"),
  Q = NULL,
  nonneg = FALSE,
  preproc_x = multivarious::pass(),
  preproc_y = multivarious::pass(),
  tol = 1e-06,
  maxiter = 200,
  verbose = FALSE,
  ...
)

Arguments

Details

Unlike genpls(), which handles separate row and column metrics (Mx, Ax, My, Ay) with a Gram–Schmidt orthogonalisation step, rpls() uses a single metric Q and the simpler penalised updates of Allen et al.

Value

An object of class c("rpls","cross_projector","projector") with at least the elements

vx

p \times K matrix of X-loadings.

vy

q \times K matrix of Y-loadings.

ncomp

Number of components actually extracted (may be < K).

penalty

Penalty type used ("l1" or "ridge").

lambda

The lambda value(s) used.

tol

The convergence tolerance used.

maxiter

The maximum number of inner iterations used per component.

nonneg

Logical flag for the non-negativity constraint on l1.

Q_used

Logical; TRUE if a custom Q metric was supplied to induce generalised RPLS, FALSE for standard (identity-metric) RPLS. Note the field is named Q_used, not Q: the metric matrix itself is not retained on the returned object.

verbose

The verbose flag used.

preproc_x, preproc_y

Pre-processing transforms used.

rpls objects store no row (X-)scores: the factors z_k are computed internally during fitting to drive deflation and are then discarded, so use multivarious::project() on X to obtain scores for any rows of interest.

The object supports predict(), project(), transfer(), coef() and other multivarious generics.

Method

The routine follows Algorithm 1 of Allen et al. (2013, Stat. Anal. Data Min., 6 : 302–314) — see the paper for details. Briefly, with C = X^\top Y the cross-product of the (preprocessed) blocks, each component maximises

\max_{u,v}\; v^\top Q C u - \lambda \, P(v)

with Q = I_p for standard RPLS. The alternating updates are: u \leftarrow C^\top Q v / \|C^\top Q v\|_2, then a penalised (possibly non-negative) regression for v, normalised in the Q-norm.

References

Allen, G. I., Peterson, C., Vannucci, M., & Maletić-Savatić, M. (2013). Regularized Partial Least Squares with an Application to NMR Spectroscopy. Statistical Analysis and Data Mining, 6(4), 302-314. DOI:10.1002/sam.11169.

Examples

# Generate sample data
set.seed(123)
n <- 50
p <- 20
q <- 10
X <- matrix(rnorm(n * p), n, p)
Y <- X[, 1:5] %*% matrix(rnorm(5 * q), 5, q) + matrix(rnorm(n * q), n, q)

# Fit regularized PLS with L1 penalty
fit_l1 <- rpls(X, Y, K = 3, lambda = 0.1, penalty = "l1")
print(fit_l1)

# Fit regularized PLS with ridge penalty
fit_ridge <- rpls(X, Y, K = 3, lambda = 0.1, penalty = "ridge")
print(fit_ridge)

Sparse and Functional Principal Components Analysis (SFPCA) with Spatial Coordinates

Description

Performs Sparse and Functional PCA on a data matrix, allowing for both sparsity and smoothness in the estimated principal components. Penalty parameters left NULL are selected automatically (see Details). The spatial smoothness penalty is constructed based on provided spatial coordinates.

Usage

sfpca(
  X,
  K,
  spat_cds,
  lambda_u = NULL,
  lambda_v = NULL,
  alpha_u = NULL,
  alpha_v = NULL,
  Omega_u = NULL,
  penalty_u = "l1",
  penalty_v = "l1",
  nlambda = 10,
  lambda_min_ratio = 0.01,
  knn = min(6, ncol(X) - 1),
  max_iter = 100,
  tol = 1e-06,
  verbose = FALSE,
  uthresh = NULL,
  vthresh = NULL
)

Arguments

Details

Each rank-1 problem is solved by alternating solves of the penalized quadratic subproblems (via C++ coordinate descent) followed by rescaling onto the smoothness-metric ball, in the constraint form of Allen & Weylandt (2019). For the convex "l1" penalty with subproblems solved to tolerance (the internal exact_inner = TRUE path, used by the monotonicity test) the objective is monotonically non-decreasing; the default inexact path tightens the inner tolerance to a floor before it may declare convergence, reproducing the same terminal iterates but without an every-iteration monotonicity guarantee (it may also stop at max_iter).

When lambda_u or lambda_v is NULL it is selected per component by a BIC-style criterion along a regularization path. For the convex "l1" penalty lambda_max = max(abs(b)) is, in closed form, the smallest value whose subproblem solution is exactly zero (at x = 0 the ⁠S x⁠ term vanishes, so the KKT condition ⁠|b_j| <= lambda⁠ does not depend on S); where b is the matrix-vector product with the other factor fixed at the SVD initializer. For the non-convex "scad" penalty the same value anchors the path but is not a global-optimality threshold. nlambda values are laid log-spaced down to lambda_min_ratio * lambda_max, coordinate descent is warm-started along the path, and the value minimizing ⁠log(RSS / (n p)) + df * log(n p) / (n p)⁠ is chosen, with df the support size of the solution and RSS the one-sided rank-1 residual sum of squares with the opposite factor held fixed (a selection heuristic, not the BIC of the fully alternated rank-1 model). The all-zero solution (at lambda_max) is a legitimate candidate: if no rank-1 structure justifies its degrees of freedom, the component is returned as exactly zero with d = 0.

When alpha_u or alpha_v is NULL it defaults to 1 / lambda_max(Omega), so the roughest direction of the smoothness penalty is weighted exactly as strongly as the identity term. This makes the default invariant to the scaling of Omega and bounds the condition number of every subproblem system I + alpha * Omega by 2.

Value

An object of class c("sfpca", "bi_projector") from the multivarious framework. Use multivarious::scores() for the sample scores (U D), multivarious::components() for the sparse loadings V, multivarious::sdev() for d_k, and multivarious::reconstruct() for the rank-K approximation. ov (like components()) holds the sparse right factors V; ou holds the left factors U. The selected penalty parameters are stored as lambda_u, lambda_v, alpha_u, and alpha_v. For backward compatibility the pre-0.1 list fields ⁠$d⁠ (singular values) and ⁠$u⁠ (left factors) remain readable but emit a deprecation warning; use sdev() and scores()/⁠$ou⁠ instead.

Important: unlike genpca(), the columns of U (ou) and V (ov) are Euclidean unit-norm but are not mutually orthogonal across components – sfpca() extracts each rank-1 term from a constraint-form subproblem rather than a joint SVD, so ⁠U'U != I⁠ and ⁠V'V != I⁠ in general. Consequently multivarious::sdev() here is not the singular values of X; it is the per-component captured covariance d_k = u_k' X_k v_k, where X_k is the matrix after the preceding components have been deflated out (so the identity holds against X itself only for k = 1). This non-orthogonality is also why reconstruct() for "sfpca" objects uses the stored U, d, V factors directly (⁠U D V'⁠) rather than SVD-based identities such as the Moore-Penrose pseudoinverse of the loadings, which would not reproduce the fitted model for non-orthogonal V (see reconstruct.sfpca()).

References

Allen, G. I., & Weylandt, M. (2019). Sparse and functional principal components analysis. In 2019 IEEE Data Science Workshop (DSW) (pp. 11-16). doi:10.1109/DSW.2019.8755778. Also available as arXiv:1309.2895, first posted in 2013 and revised through 2019; the preprint and the DSW paper are the same work, which is why both years appear in the literature.

See Also

genpca() for the shared multivarious verbs; multivarious::bi_projector.

Examples

library(Matrix)
set.seed(123)
# Smooth temporal factor, sparse spatial factor
n <- 100  # Number of time points
p <- 50   # Number of spatial locations
u <- sin(seq(0, 2 * pi, length.out = n))
v <- c(rnorm(10), rep(0, p - 10))
X <- 8 * tcrossprod(u / sqrt(sum(u^2)), v / sqrt(sum(v^2))) +
  matrix(rnorm(n * p, sd = 0.2), n, p)
spat_cds <- matrix(runif(p * 3), nrow = 3, ncol = p)  # 3D coordinates
result <- sfpca(X, K = 1, spat_cds = spat_cds)
multivarious::sdev(result)                  # captured covariance (BIC-tuned)
sum(multivarious::components(result) != 0)  # sparse spatial loading

Coordinate descent for the SFPCA penalized quadratic subproblem

Description

Internal solver for ⁠min_x 0.5 x'Sx - b'x + P(x; lambda)⁠ with sparse SPD S and an L1 or SCAD penalty.

Usage

sfpca_cd_solve_cpp(S, b, x0, lambda, penalty, scad_a, max_sweeps, tol)

Arguments

Value

list with x, sweeps, and kkt (max KKT residual)

Fast sub-space solver for a small block of generalized eigen-pairs

Description

Uses pre-conditioned sub-space iteration on the operator S_2^{-1} S_1 (or its inverse) to obtain the q largest or smallest generalized eigen-values/vectors of S_1 v = \lambda S_2 v.

Usage

solve_gep_subspace(
  S1,
  S2,
  q = 2,
  which = c("largest", "smallest"),
  max_iter = 100,
  tol = 1e-06,
  V0 = NULL,
  seed = NULL,
  reg_S = 0.001,
  reg_T = 1e-06,
  verbose = FALSE
)

Arguments

Value

A list with components

values

length-q numeric vector of Ritz eigen-values.

vectors

⁠d x q⁠ matrix, columns are orthonormal eigen-vectors in the original S-inner-product.

Symmetrize a nearly symmetric matrix or stop

Description

Measures ⁠||A - A'||_F / ||A||_F⁠. Below rtol the two triangles are averaged and the result is marked symmetric; above it the function stops. Asymmetry is an input error, never something a PSD remedy repairs.

Usage

symmetrize_or_stop(A, rtol = 1e-10, name = "A")

Arguments

Value

a symmetric Matrix (dense or sparse as supplied)

Compatibility wrapper for multivarious::transfer

Description

Provides a transfer method that accepts source/target arguments used in the tests. The method forwards to multivarious's transfer method which expects from/to.

Usage

## S3 method for class 'cross_projector'
transfer(x, new_data, from = NULL, to = NULL, opts = list(), ...)

Arguments

Value

Matrix with transferred data.

Examples

# Generate sample data
set.seed(123)
n <- 50
X <- matrix(rnorm(n * 10), n, 10)
Y <- matrix(rnorm(n * 8), n, 8)

# Create a cross projector using rpls
fit <- rpls(X, Y, K = 2)

# Transfer new X data to Y space
new_X <- matrix(rnorm(10 * 10), 10, 10)
transferred <- transfer(fit, new_X, from = "X", to = "Y")

Truncate a projection to fewer components

Description

Re-exported from multivarious. See truncate for details.

Usage

truncate(x, ncomp)

Arguments

Value

Truncated projection object

Truncate a genpca fit to fewer components

Description

Returns a new genpca object retaining only the first ncomp components. All component-indexed slots (v, s, sdev, ov, ou, u, propv, cumv) are sliced consistently; the preprocessing object and constraint matrices are carried over unchanged.

Usage

## S3 method for class 'genpca'
truncate(x, ncomp)

Arguments

Value

A genpca object with ncomp components.

See Also

genpca(), reconstruct.genpca()

Examples

X <- matrix(rnorm(60), 15, 4)
fit <- genpca(X, ncomp = 4)
fit2 <- truncate(fit, 2)
multivarious::ncomp(fit2)