Project Features onto a Spectral Basis

Projects a feature matrix onto a given spectral basis. Handles potential transposition of the feature matrix and uses the appropriate orthogonal projection formula based on whether the basis is orthonormal. Sparse inputs are processed with `Matrix` methods so that large matrices need not be fully densified. The returned matrix is always dense.

Usage

project_features_to_spectral_space(
  feature_matrix,
  U_basis,
  tol_orthonormal = sqrt(.Machine$double.eps),
  assume_orthonormal = FALSE
)

Arguments

feature_matrix

A numeric matrix representing features. It can be in `V_p x C` format (parcels x features/conditions) or `C x V_p` format. The function will attempt to orient it correctly based on `U_basis`. If `feature_matrix` is a sparse `Matrix`, operations use sparse methods when possible (coercion occurs only for small intermediate matrices). Very large sparse inputs (`prod(dim(.)) > 1e7`) trigger an error to avoid accidental densification.

U_basis

A numeric matrix representing the spectral basis, typically with dimensions `V_p x k_dims_basis` (parcels x basis dimensions). Sparse `U_basis` inputs are handled with `Matrix` methods. Extremely large sparse bases also trigger an error to prevent densification.

tol_orthonormal

A numeric tolerance to check for orthonormality of `U_basis`. `max(abs(crossprod(U_basis) - I)) / k_dims_basis` is compared against this tolerance. Default is `sqrt(.Machine$double.eps)`.

assume_orthonormal

Logical. If `TRUE`, `U_basis` is assumed to be orthonormal, and the check (including rank check and `UtU` calculation) is skipped, directly using the faster projection `crossprod(U_basis, features)`. Default is `FALSE`.

Value

A numeric (dense) matrix containing the projected features. If `feature_matrix` was oriented to `V_p x C` (parcels x features/conditions), the output will be `k_dims_basis x C` (basis dimensions x features/conditions). If `feature_matrix` was `C x V_p` and was transposed for projection, the output will be `C x k_dims_basis` to maintain the original feature orientation as rows.

Examples

V_p <- 10
k_dims <- 3
C <- 5
U_basis_ortho <- svd(matrix(rnorm(V_p * k_dims), V_p, k_dims))$u
U_basis_non_ortho <- matrix(rnorm(V_p * k_dims), V_p, k_dims)

# Feature matrix: V_p x C (parcels x conditions)
features1 <- matrix(rnorm(V_p * C), V_p, C)
proj1_ortho <- project_features_to_spectral_space(features1, U_basis_ortho)
print(dim(proj1_ortho))
#> [1] 3 5
proj1_non_ortho <- project_features_to_spectral_space(features1, U_basis_non_ortho)
print(dim(proj1_non_ortho))
#> [1] 3 5

# Sparse input example
features_sp <- Matrix::rsparsematrix(V_p, C, 0.2)
proj_sp <- project_features_to_spectral_space(features_sp, U_basis_ortho)
#> Error in project_features_to_spectral_space(features_sp, U_basis_ortho): feature_matrix must be a numeric matrix or Matrix object.
print(class(proj_sp))
#> Error: object 'proj_sp' not found

# Ambiguous square matrix error
V_p_sq <- 7
U_basis_sq <- matrix(rnorm(V_p_sq*k_dims), V_p_sq, k_dims)
features_sq_ambiguous <- matrix(rnorm(V_p_sq*V_p_sq), V_p_sq, V_p_sq)
try(project_features_to_spectral_space(features_sq_ambiguous, U_basis_sq))
#> Error in project_features_to_spectral_space(features_sq_ambiguous, U_basis_sq) : 
#>   feature_matrix dimensions (7 x 7) are ambiguous as both match V_p_basis (7). Please ensure feature_matrix is oriented as V_p x C (parcels x features/conditions). If unsure, or if it is C x V_p where C == V_p, please transpose it first before calling.

# Rank deficient U_basis warning
U_rank_def <- matrix(rnorm(V_p * k_dims), V_p, k_dims)
if (k_dims > 1) U_rank_def[, k_dims] <- U_rank_def[, 1] # Make last col same as first
try(project_features_to_spectral_space(features1, U_rank_def))
#> Warning: U_basis appears rank deficient (rank 2 / 3 dimensions). Projection results may be unreliable or reflect a lower-dimensional space.
#> Error in value[[3L]](cond) : 
#>   Error in qr.solve, likely UtU is singular or ill-conditioned (U_basis may be rank-deficient, or k_dims_basis > V_p_basis and U_basis not full rank). Original error: singular matrix 'a' in solve