Matrix Variate Bilinear MFA

• Matrix Variate Bilinear MFA is a latent factor model for matrix data that preserves two-dimensional structure using separate row and column loading matrices.
• It employs a two-step spectral decomposition to estimate low-rank factors efficiently, ensuring robust dimension reduction and prediction.
• Extensions include regression models like LaGMaR, mixture models for clustering, and robust methods to handle skewed or heavy-tailed data.

Matrix variate bilinear MFA refers to a class of latent factor models for matrix-valued data in which the signal is generated through a bilinear low-rank structure, typically for dimension reduction, feature extraction, prediction, or model-based clustering. Unlike vectorization approaches that discard the inherent two-dimensional geometry of the data, these models preserve the matrix structure through a decomposition employing separate loading matrices for rows and columns. This paradigm underpins a range of methodologies—from classical unsupervised matrix-variate bilinear factor analyzers to recent regression-based and robust inference extensions.

1. Core Bilinear Matrix Factor Analysis Model

The foundational model observes $n$ independent random matrices $X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$, and expresses each as a bilinear generative process: $X_i = R\,Z_i\,C^{\top} + E_i,$ where $R\in\mathbb{R}^{p_1\times k_1}$ and $C\in\mathbb{R}^{p_2\times k_2}$ are low-rank row and column loading matrices ($k_1\ll p_1$, $k_2\ll p_2$), $Z_i\in\mathbb{R}^{k_1\times k_2}$ is the latent factor matrix for observation $i$, and $E_i$ represents idiosyncratic noise. The typical distributional assumption is matrix-variate normality,

$X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$0

but extensions to skewed and heavy-tailed settings utilize mixtures or matrix-variate $X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$1 distributions (Gallaugher et al., 2017, Gallaugher et al., 2018, Ma et al., 2024).

Consistency of the bilinear form is fundamental: under high-dimensional asymptotics and “pervasive” factors, estimators of the factor scores $X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$2 and loading spaces $X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$3 are consistent up to orthogonal transformation (Zhang et al., 2022).

2. Principal Component–Based Estimation

The link between bilinear MFA and high-dimensional principal component analysis is operationalized by leveraging the following two-step spectral decomposition:

• Row Direction: Form the empirical “column-wise covariance”

$X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$4

and obtain its $X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$5 leading eigenvectors as $X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$6.

• Column Direction: Similarly, use

$X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$7

to estimate $X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$8.

Each latent $X_1,\ldots,X_n \in \mathbb{R}^{p_1\times p_2}$9 is then reconstructed via

$X_i = R\,Z_i\,C^{\top} + E_i,$0

Factor numbers $X_i = R\,Z_i\,C^{\top} + E_i,$1 are identified by a ratio-of-eigenvalues criterion, reducing the need for extensive parameter tuning seen in penalized vector regression (Zhang et al., 2022). This preserves the two-dimensional geometry and enables a tissue-specific factorization in imaging or spatiotemporal data applications.

3. Extensions: Regression, Mixtures, and Robustness

3.1. Generalized Regression and LaGMaR

The latent generalized matrix regression (LaGMaR) model incorporates matrix-variate bilinear MFA scores $X_i = R\,Z_i\,C^{\top} + E_i,$2 as predictors of an outcome variable $X_i = R\,Z_i\,C^{\top} + E_i,$3, potentially exponential family-distributed, via the generalized linear model: $X_i = R\,Z_i\,C^{\top} + E_i,$4 Here, $X_i = R\,Z_i\,C^{\top} + E_i,$5, $X_i = R\,Z_i\,C^{\top} + E_i,$6 denotes Frobenius inner product, and $X_i = R\,Z_i\,C^{\top} + E_i,$7 are additional covariates. The approach achieves dimension reduction without large-scale penalization, and Kullback–Leibler–consistent prediction—despite $X_i = R\,Z_i\,C^{\top} + E_i,$8 and $X_i = R\,Z_i\,C^{\top} + E_i,$9 being identified only up to rotation (Zhang et al., 2022).

3.2. Mixture Models and Cluster Analysis

Mixture models extend bilinear MFA to clustering and classification. Each component $R\in\mathbb{R}^{p_1\times k_1}$0 has its own mean $R\in\mathbb{R}^{p_1\times k_1}$1, loadings $R\in\mathbb{R}^{p_1\times k_1}$2, and uniqueness parameters; the marginal density is

$R\in\mathbb{R}^{p_1\times k_1}$3

with identifiability up to rotations of the factor spaces (Gallaugher et al., 2017, Gallaugher et al., 2019). Parsimonious mixtures impose various constraints on loadings and uniquenesses, leading to an 8×8 grid of models (Gallaugher et al., 2019). Skewed and heavy-tailed mixtures are realized through variance-mean mixture (e.g., matrix-variate skew-$R\in\mathbb{R}^{p_1\times k_1}$4 or GH components) (Gallaugher et al., 2018).

3.3. Robust Bilinear Factor Analysis

Recent work embeds the bilinear structure in the matrix-variate $R\in\mathbb{R}^{p_1\times k_1}$5 distribution: $R\in\mathbb{R}^{p_1\times k_1}$6 enabling robust inference even under contamination or heavy tails. The $R\in\mathbb{R}^{p_1\times k_1}$7BFA model attains a breakdown point substantially higher than vectorized $R\in\mathbb{R}^{p_1\times k_1}$8FA, as the joint decomposition of row/column covariances reduces the effective dimension determining robustness (Ma et al., 2024).

4. Algorithmic Implementation

Matrix variate bilinear MFA estimation routinely consists of the following pipelines:

• Spectral Estimation: Two leading eigendecompositions for row/column covariance matrices. This step is tuning-free aside from selection of $R\in\mathbb{R}^{p_1\times k_1}$9.
• EM or AECM: For mixture or robust models, alternating expectation-conditional maximization (AECM) cycles treat different latent missing data (e.g., factors, cluster allocations, scale variables) in blockwise manner.
• Closed-form Updates: Many parameter updates admit explicit formulas, with the majority of computational cost in linear algebra operations (e.g., eigenanalysis, matrix multiplications).

A summary of computational steps for LaGMaR (Zhang et al., 2022):

Step	Operation	Complexity
Compute $C\in\mathbb{R}^{p_2\times k_2}$0, $C\in\mathbb{R}^{p_2\times k_2}$1	Empirical covariances	$C\in\mathbb{R}^{p_2\times k_2}$2
Eigen-decompose	Row/column spectral decompositions	$C\in\mathbb{R}^{p_2\times k_2}$3
Extract $C\in\mathbb{R}^{p_2\times k_2}$4	Matrix multiplications	$C\in\mathbb{R}^{p_2\times k_2}$5
GLM fit on latent scores	Standard GLM methods	depends on $C\in\mathbb{R}^{p_2\times k_2}$6

Mixures and robust models proceed similarly, with additional EM/AECM steps dictated by the structure of latent variables (e.g., scale factors $C\in\mathbb{R}^{p_2\times k_2}$7 in $C\in\mathbb{R}^{p_2\times k_2}$8BFA) (Ma et al., 2024, Gallaugher et al., 2019).

5. Theoretical Properties and Consistency

The central theoretical guarantees for matrix variate bilinear MFA include bilinear-form consistency, coefficient (regression) consistency up to rotation, and prediction consistency:

• Bilinear-form consistency: For some orthogonal $C\in\mathbb{R}^{p_2\times k_2}$9,

$k_1\ll p_1$0

as $k_1\ll p_1$1 if factors are “pervasive” (Zhang et al., 2022).

• Coefficient consistency: The regression coefficient $k_1\ll p_1$2 (or its mixture analog) is consistently estimated up to the same rotations.
• Prediction consistency: Predicted outcomes from the fitted model converge in probability to the true conditional mean $k_1\ll p_1$3 even though $k_1\ll p_1$4 are only identified up to orthogonal transforms.
• Robustness: In $k_1\ll p_1$5BFA, the breakdown point is at least $k_1\ll p_1$6 or $k_1\ll p_1$7, which dominates the $k_1\ll p_1$8 bound of vectorized $k_1\ll p_1$9FA, implying substantial robustness gains for matrix-valued data (Ma et al., 2024).

Rates of convergence for estimated factors follow those of vector factor models, typically $k_2\ll p_2$0 where $k_2\ll p_2$1 (Zhang et al., 2022).

6. Applications, Strengths, and Limitations

Matrix variate bilinear MFA and its extensions are widely adopted for:

• Imaging and medical diagnosis: LaGMaR was motivated by 2D CT image biomarkers for COVID-19 status prediction, offering dimension reduction that preserves spatial structure without costly penalization (Zhang et al., 2022).
• High-dimensional clustering/classification: PMMVBFA and MMVBFA deliver accurate clustering and semi-supervised classification in scenarios such as MNIST and face recognition (Gallaugher et al., 2017, Gallaugher et al., 2019).
• Robust inference: $k_2\ll p_2$2BFA attains higher resilience to outliers and heavy tails in financial and biomedical contexts, where classical Gaussian-based models break down (Ma et al., 2024).

Key strengths:

• Structural respect for 2D matrix geometry, avoiding flattening-induced information loss.
• Tuning-free or minimal-tuning estimation in leading PCA-based approaches.
• Closed-form and computationally efficient implementation (especially for unsupervised and regression variants).

Limitations:

• Requires strong low-rank separability in signal; if bilinear factor structure is violated, estimation may fail.
• Weakly correlated noise and pervasive factor assumptions are necessary for consistency.
• Ratio-of-eigenvalues factor selection may not reliably distinguish weak factors.

A plausible implication is that future work on matrix-variate bilinear MFA will focus on relaxing separation/model assumptions, advancing robustifications, and scaling high-throughput algorithms for very large $k_2\ll p_2$3 encountered in contemporary imaging and genomics.