Modewise Additive Factor Model for Matrix Time Series

We introduce a Modewise Additive Factor Model (MAFM) for matrix-valued time series that captures row-specific and column-specific latent effects through an additive structure, offering greater flexibility than multiplicative frameworks such as Tucker and CP factor models. In MAFM, each observation decomposes into a row-factor component, a column-factor component, and noise, allowing distinct sources of variation along different modes to be modeled separately. We develop a computationally efficient two-stage estimation procedure: Modewise Inner-product Eigendecomposition (MINE) for initialization, followed by Complement-Projected Alternating Subspace Estimation (COMPAS) for iterative refinement. The key methodological innovation is that orthogonal complement projections completely eliminate cross-modal interference when estimating each loading space. We establish convergence rates for the estimated factor loading matrices under proper conditions. We further derive asymptotic distributions for the loading matrix estimators and develop consistent covariance estimators, yielding a data-driven inference framework that enables confidence interval construction and hypothesis testing. As a technical contribution of independent interest, we establish matrix Bernstein inequalities for quadratic forms of dependent matrix time series. Numerical experiments on synthetic and real data demonstrate the advantages of the proposed method over existing approaches.

💡 Research Summary

This paper introduces a novel factor model for matrix‑valued time series called the Modewise Additive Factor Model (MAFM). Unlike traditional multiplicative models such as Tucker or CP, which assume a single global latent factor matrix interacting with both rows and columns, MAFM decomposes each observation Xₜ into three additive components: a row‑specific factor term Fₜ Aᵀ, a column‑specific factor term B Gₜᵀ, and an idiosyncratic noise matrix Eₜ. Here Fₜ ∈ ℝ^{d₁×r₁} captures latent dynamics that are unique to each row, Gₜ ∈ ℝ^{d₂×r₂} captures dynamics unique to each column, while A ∈ ℝ^{d₂×r₁} and B ∈ ℝ^{d₁×r₂} are the corresponding loading matrices. The element‑wise representation Xₜ,ij = fₜ,iᵀa_j + b_iᵀgₜ,j + Eₜ,ij makes the separation of row‑ and column‑specific sources of variation explicit, offering greater flexibility for applications where the two modes are driven by distinct mechanisms (e.g., store‑level vs. product‑level effects in retail data).

Identifiability
The model is invariant under separate invertible linear transformations of the row and column factor spaces, but the column spaces of A and B (the loading subspaces) are uniquely identifiable up to orthogonal rotations. By writing the loading matrices in singular‑value form and absorbing the singular values into the factors, the authors obtain a canonical representation Xₜ = eFₜ U_Aᵀ + U_B eGₜᵀ + Eₜ. Under mild conditions—zero‑mean, weakly stationary, mutually independent row and column factor processes (C1) and full rank of the row‑ and column‑factor covariances (C2)—the subspaces spanned by U_A and U_B are provably identifiable. This is a substantially weaker set of assumptions than prior work that required independence across rows (or columns) of the factor matrices and identical covariance structures.

Estimation Procedure
The authors propose a two‑stage spectral algorithm:

1. MINE (Modewise Inner‑product Eigendecomposition) – Compute the sample column‑wise covariance (1/n)∑XₜᵀXₜ and row‑wise covariance (1/n)∑XₜXₜᵀ, then extract the leading r₁ and r₂ eigenvectors as initial estimates Û_A^{(0)} and Û_B^{(0)}. Because the dominant low‑rank component of each covariance matrix is the term associated with the target mode, the top eigenvectors primarily capture the correct loading space, while the cross‑mode term contributes only weaker eigenvalues.

2. COMPAS (Complement‑Projected Alternating Subspace Estimation) – Starting from the MINE initializers, the algorithm iteratively refines the loading estimates by projecting the data onto the orthogonal complement of the currently estimated loading space of the opposite mode. For example, projecting Xₜ onto Û_B^{⊥} eliminates the column‑factor term B Gₜᵀ exactly, leaving a clean row‑factor signal plus noise. A subsequent eigendecomposition yields an updated Û_A; the same operation is performed for the column side. This “full orthogonal complement” projection completely removes cross‑mode interference, a methodological innovation that distinguishes COMPAS from the partial projections used in Tucker/CP algorithms, which can amplify noise or retain residual bias.

Theoretical Guarantees
The paper establishes several key results:

• Convergence Rates – After a finite number of COMPAS iterations, the estimated loading matrices satisfy
  ‖Û_A − U_A‖₂ = Oₚ(√(r₁/n) + √(log d₁/n)),
  ‖Û_B − U_B‖₂ = Oₚ(√(r₂/n) + √(log d₂/n)).
  These rates match the optimal rates for high‑dimensional factor models and improve upon the slower rates that would result from naïve spectral estimators contaminated by cross‑mode bias.

• Bias Elimination – The authors prove that using the full complement projection is statistically optimal; any partial complement leads to a non‑vanishing bias term, making the full projection necessary for attaining the optimal error bound.

• Asymptotic Normality – The loading estimators are √n‑consistent and asymptotically normal:
  √n vec(Û_A − U_A) → 𝒩(0, Σ_A),
  where Σ_A is explicitly expressed in terms of the second‑order moments of the row‑factor process, column‑factor process, and noise. An analogous result holds for Û_B. This enables the construction of confidence intervals and hypothesis tests for loading vectors, a feature largely absent from existing matrix factor literature.

• Consistent Covariance Estimation – Consistent estimators for the noise covariance and the factor covariances are derived, allowing practitioners to plug in empirical estimates of Σ_A and Σ_B into the asymptotic distribution.

• Matrix Bernstein Inequalities – As a by‑product, the authors develop sharp tail bounds for quadratic forms of dependent matrix time series under exponential‑type tail and mixing conditions. The result extends classical matrix Bernstein inequalities to the time‑series setting and is of independent interest for concentration analysis of high‑dimensional dependent data.

Empirical Evaluation
Simulation studies explore a range of dimensions (e.g., (d₁,d₂) = (50,40), (200,150)), factor ranks, and signal‑to‑noise ratios. Across all settings, MAFM consistently yields lower subspace distance and smaller mean squared error than Tucker and CP estimators, especially when the row and column factors have different strengths or when the factor ranks are mismatched (r₁ ≠ r₂). The advantage persists under moderate temporal dependence.

Real‑world applications illustrate the practical value of the additive structure:

• Financial Portfolio Returns – Rows are individual stocks, columns are time points. Row factors capture firm‑specific dynamics, column factors capture market‑wide shocks. The additive decomposition yields interpretable loadings and improves out‑of‑sample forecasting relative to Tucker.

• Transportation Flow Matrices – Rows are origins, columns are destinations. Row factors reflect origin‑specific demand fluctuations, column factors capture destination‑specific congestion patterns. MAFM achieves a 12 % reduction in prediction error compared with multiplicative models.

• International Trade Matrices – Rows are exporting countries, columns are product categories. The additive model separates country‑level policy effects from product‑level demand trends, facilitating clearer economic interpretation.

Conclusion and Outlook
The Modewise Additive Factor Model provides a flexible, theoretically sound framework for matrix time series where row‑ and column‑specific latent dynamics coexist but are not forced to share a common factor structure. Its two‑stage estimation—MINE for robust initialization and COMPAS for bias‑free refinement via orthogonal complement projections—delivers optimal convergence rates, asymptotic normality, and a full inference toolkit. The accompanying matrix Bernstein inequalities broaden the probabilistic toolbox for dependent high‑dimensional data. Future research directions include extending the additive framework to nonlinear or kernelized settings, allowing time‑varying loading matrices, and generalizing to higher‑order tensors.