PaperlayerView exact revision on arXiv

Source-linked AI summary

Sparse Separable Factor Analysis in the Complex Domain with an Application to Local Field Potential Data

Ian Hultman, Kirtikanth Kalapatapu, Yassine Filali, Rainbo Hultman, Sanvesh Srivastava

arXiv:2608.21551v1stat.MLcs.LGstat.COstat.ME

TL;DR

Covariance estimation for complex-valued arrays must preserve amplitude-phase information and multiway organization, which existing real-domain or vectorization-based approaches do not directly exploit. The paper develops SSFA, combining separable complex factor modeling, sparse loading matrices, and mode-wise PX-EM with phase-preserving soft-thresholding and scale balancing. Simulations show improved covariance estimation over vectorization-based methods, while mouse LFP analyses compare separability structures and impute recordings missing from electrode misplacement.

  • Problem

    Existing covariance methods either ignore complex-array organization or use real-domain embeddings, limiting direct exploitation of amplitude, phase, and multiway structure.

  • Method

    SSFA models each mode-specific covariance with low-rank Hermitian factors and diagonal residuals, estimates sparse complex loadings by mode-wise PX-EM, and balances separable covariance scales.

  • Results

    Across simulation settings, SSFA outperforms vectorization-based methods, including complex PCA, and is applied to separability comparison and model-based LFP imputation.

  • Takeaways & Limitations

    SSFA provides interpretable sparse covariance estimates that retain complex amplitude-phase information for multiway data and supports frequency-dependent LFP dependence analysis and missing-recording imputation.

  • Takeaways & Limitations

    Repeated mode-wise whitening is the main computational bottleneck, and the model assumes proper complex Gaussian data in its foundational formulation.

Abstract

from arXiv · show

Complex-valued arrays arise in signal processing, where scientific interpretation depends on retaining amplitude and phase information. Existing covariance estimation methods either ignore the multiway organization of such data or rely on real-domain embeddings that do not directly exploit their complex structure. We develop sparse separable factor analysis (SSFA), a latent factor model for complex-valued arrays with a separable covariance structure across modes. Each mode-specific covariance matrix is modeled through a low-rank Hermitian factor structure and a diagonal residual covariance matrix. To obtain interpretable estimates, we impose elementwise lasso penalties on the complex loading matrices and estimate the SSFA parameters using a mode-wise parameter-expanded expectation-maximization procedure. The resulting loading updates admit closed-form complex soft-thresholding solutions, which shrink the modulus of each loading while preserving its phase. A separate balancing step resolves the scale nonidentifiability of the separable covariance structure. Simulation studies show that SSFA improves covariance estimation relative to vectorization-based methods, including complex principal component analysis. We apply SSFA to local field potential recordings from mice, where we compare separability structures induced by different groupings of brain region, frequency, and time and perform model-based imputation of recordings missing because of electrode misplacement.

1 Introduction

SSFA addresses covariance estimation for complex-valued arrays by preserving amplitude-phase information and multiway organization while combining separable factor structure with sparse complex loadings. It introduces a lasso-regularized PX-EM framework and applies it to covariance estimation and missing-recording imputation in LFP data.

  • Research gap: Existing SFA models are real-valued with dense loadings, while sparse factor models and complex factor methods primarily target real-valued or vector-valued data.Vectorization also disregards multiway organization and can produce unstable, difficult-to-interpret estimates when dimensions are large relative to sample size.
  • Motivation: Complex-valued LFP and EEG representations retain amplitude and phase information that real-domain embeddings do not directly exploit.Sparse mode-specific loadings can identify contributing brain regions and frequency bands while producing parsimonious low-rank covariance estimates.
  • Method: SSFA models complex-valued arrays with separable covariance across modes and sparse complex loading matrices using low-rank Hermitian factors plus diagonal residual covariances.A balancing condition resolves multiplicative scale nonidentifiability across mode-specific covariance matrices.
  • Estimation: Mode-wise PX-EM whitens other modes, reduces each update to complex matrix factor analysis, and uses closed-form complex soft-thresholding for the penalized loading update.The procedure uses a common regularization grid across modes and a plug-in residual covariance update.
  • Applications: SSFA supports model-based imputation of recordings missing from electrode misplacement and comparison of separability structures involving brain region, frequency, and time.The brain-region–frequency by time matricization yields the most interpretable dependence structure in the application.
  • Contributions: Simulation studies and LFP analyses demonstrate SSFA's practical utility, including covariance-structure comparison and imputation of missing recordings.Across simulation settings, SSFA outperforms vectorization-based methods, including complex PCA.

2 Motivating application

The motivating application uses Fourier-transformed mouse LFP recordings to estimate frequency- and time-dependent dependence among brain regions while retaining complex coefficients. SSFA produces coherence and phase-offset estimates, compares separability structures, and supports imputation of recordings lost to electrode misplacement.

  • Data and goals: Fourier-transformed LFP recordings form complex arrays indexed by brain region, frequency, and time window.The scientific goal is to characterize frequency-specific dependence between brain regions and its evolution over time.
  • Data representation: Modeling Fourier coefficients directly preserves their complex representation and avoids separately constructing amplitude and dependence summaries.For 13 regions and 100 frequencies, this uses 1300 complex coefficients rather than 9100 real-valued power and coherence features per time window.
  • Dependence measures: Estimated covariance matrices yield coherence and phase-offset matrices, with coherence measuring association strength and phase offsets describing relative phase.Phase offsets are less informative when the corresponding coherence is small and should be interpreted jointly with coherence.
  • Cross-frequency dependence: Figure 2 compares within-frequency coherence at 1 and 2 Hz with their canonical cross-coherence, which is substantially weaker.This pattern is consistent with the expected near-orthogonality of distinct Fourier frequencies.
  • Imputation: SSFA supports imputation of missing brain-region recordings under a missing-completely-at-random assumption, avoiding information loss from discarding affected regions or animals.The missingness arises when electrodes are misplaced or record from unintended locations.

3 Methodology for complex-valued vector data

The vector-data foundation is a proper complex factor model whose covariance combines a low-rank loading component with diagonal residual variation. Sparse PX-EM estimation adds complex lasso regularization while preserving phase information in loading updates.

  • Factor model: The complex factor model represents observations as x = Λz + e, with low-dimensional latent factors and diagonal residual errors.Here, x is observed, z is latent, e is residual, and k ≪ p.
  • Complex Gaussian model: Proper complex Gaussian modeling uses Hermitian covariance and zero pseudo-covariance, implying circular symmetry under global phase rotations.This restriction is connected to unitary nonidentifiability in factor representations.
  • Marginal covariance: Marginalizing the latent factors gives Σ = ΛΛ⋆ + Ψ, combining a low-rank factor covariance with diagonal residual covariance.The matrix-variate formulation stacks independent observations into X, Z, and E.
  • PX-EM: PX-EM introduces a working latent covariance, updates working parameters, and maps them back to the original loading parameterization without changing the observed covariance.The parameter-reduction map uses the Cholesky factor of the working latent covariance.
  • Sparse extension: Sparse estimation imposes an elementwise complex lasso penalty on the reduced loading matrix while retaining the PX-EM E-step and working-parameter updates.The penalty targets loading moduli rather than separately penalizing real and imaginary components.
  • Loading update: Complex soft-thresholding shrinks the modulus of each nonzero loading while preserving its phase, and all loading entries can be updated simultaneously.The residual covariance is then updated using the sparse loading estimate, preserving its sparsity structure.
  • Residual update: The plug-in residual covariance update uses unexplained empirical variance and produced more stable sparse solutions and rank selection in numerical experiments.The resulting procedure is called the sparse PX-EM procedure.

4 Methodology for complex-valued array data

The array-valued extension applies complex separable covariance modeling to LFP arrays and estimates each mode-specific covariance through low-rank Hermitian factors plus diagonal residuals. Mode-wise whitening enables PX-EM updates, while balancing resolves scale ambiguity and supports common regularization across modes.

  • Application model: Each mouse recording is a complex r × f × t array indexed by brain region, frequency, and time window.The target is dependence across regions and frequencies and its evolution over time.
  • Separable representations: Different matricizations and separability assumptions provide distinct scientific interpretations of dependence across brain region, frequency, and time.The rf × t representation is more parsimonious and interpretable when Fourier coefficients are approximately uncorrelated across time windows.
  • Array-normal model: The complex array normal model specifies a Kronecker covariance structure across modes, with each mode-specific covariance modeled as Σj = ΛjΛj⋆ + Ψj.Λj is a low-rank complex loading matrix and Ψj is diagonal.
  • Mode-wise estimation: For a fixed mode, whitening all other modes reduces the array problem to a matrix-variate factor-analysis problem.This mode-wise strategy avoids direct estimation of the full vectorized covariance and avoids introducing multiple latent components simultaneously.
  • Sparse PX-EM: SSFA applies an elementwise lasso penalty to each mode-specific loading matrix and updates mode parameters sequentially through PX-EM.The cycle repeats across modes until convergence, with E-step, penalized loading, and plug-in residual covariance updates.
  • Loading sparsity: Complex soft-thresholding preserves the phase of each nonzero mode-specific loading while shrinking its modulus.The update is applied entrywise to the reduced loading estimate.
  • Identifiability and tuning: The Kronecker covariance is identifiable only up to multiplicative rescaling, so SSFA equalizes the smallest diagonal residual entries across rescaled modes.The resulting unique balancing vector also places mode-specific parameters on a common scale for a shared penalty parameter.

5 Simulated data analysis

The simulations evaluate SSFA across vector-, matrix-, and array-valued complex observations, comparing covariance and loading-space recovery with vectorization-based competitors. SSFA’s advantage is strongest for higher-order arrays, while its mode-specific estimates and effective-rank recovery improve with sample size and array structure.

  • 5.1 Simulation setup: Six settings vary observation order and dimension, covering vector-, matrix-, and three-dimensional array-valued complex data.Each setting specifies mode-specific complex loading matrices and diagonal residual covariance matrices.
  • 5.1 Simulation setup: The robustness study generates data with the target Hermitian covariance but generally nonzero pseudo-covariance, departing from SSFA’s circular complex normal assumption.Latent real-valued standard normal variables are transformed with mode-specific complex Cholesky factors.
  • 5.2 Competing methods: SSFA is compared with empirical covariance, complex PCA, and real-domain sparse PCA after vectorization, using EBIC-based tuning and repeated simulation summaries.Sample sizes range from n ∈ {5, 10, 20, 40}, with n = 100 additionally used for vector-valued settings; each combination is replicated 100 times.
  • 5.4 Simulation results: SSFA and competitors have similar full covariance errors for vector-valued data, but SSFA has much smaller errΣ for matrix- and three-dimensional array-valued observations.The errors generally decrease as n increases, consistent with benefits from exploiting separability and mode-specific factor structure.
  • 5.4 Simulation results: Complex PCA and SSFA have comparable loading-space accuracy in vector-valued simulations, whereas sparse PCA has substantially larger err2Λ1 errors.The result indicates that real-domain sparsity does not recover the complex loading structure as accurately as direct complex-domain estimation.
  • 5.4 Simulation results: Mode-specific covariance and loading errors generally decrease with n, while effective-rank estimates improve substantially from vector-valued to higher-order array observations.For matrix- and three-dimensional arrays, estimated ranks are close to the true ranks, especially as dimension and sample size increase.

6 Real-world data analyses

The real-data analyses use complex factor models to impute missing brain-region recordings and compare SSFA covariance structures across brain region, frequency, and time. The results show that separability choices affect estimated dependence, while vectorized SSFA retains complex-valued and frequency-localized structure relative to PCA competitors.

  • Data and imputation: Nine of 15 mice had recordings missing from one or more brain regions because of electrode misplacement.
  • Data and imputation: The imputation evaluation masked one region or one pair of regions in six mice with complete recordings and predicted the omitted recordings from the remaining regions.
  • SSFA covariance structures: The three SSFA formulations differed in whether brain region and frequency were combined, separated, or paired with independent time-window observations.
  • Data and imputation: Larger-rank unpenalized vectorized complex factor models outperformed mean and empirical-covariance imputation, while k = 250 offered only marginal improvement over k = 150 at substantially higher computational cost.The subsequent analyses therefore used k = 150 for imputation.
  • SSFA covariance structures: The fully separable model imposes one brain-region coherence and phase-offset structure across frequencies and time windows, whereas X mat allows brain-region dependence to vary with frequency.
  • SSFA covariance structures: The estimated time-mode loadings indicated weaker dependence during phase 1 and stronger dependence during phase 2 without using phase labels during fitting.
  • Coherence and phase estimates: As frequency increased, coherence generally decreased and phase relationships among coherent region pairs became more localized.
  • Method comparisons: Vectorized SSFA retained stronger coherence and phase structure at 20 Hz than complex PCA and sparse PCA, while recovering similar lower-frequency patterns to complex PCA.

7 Discussion

The discussion identifies scalability and flexibility extensions for SSFA and contrasts it with Tucker-structured factor models. Repeated mode-wise whitening is the main computational bottleneck, while Tucker models target dimension reduction rather than mode-specific covariance estimation.

  • Extensions and limitations: Repeated mode-wise whitening is identified as SSFA’s main computational bottleneck because updated covariance estimates require rewhitening using the remaining modes.
  • Extensions and limitations: Prewhitening each mode once with diagonal covariance estimates could avoid repeated whitening, enable parallel mode-specific updates, and reduce computational cost.
  • Extensions and limitations: Alternating parameter updates with conditional imputation could incorporate partially observed arrays directly into SSFA model fitting.
  • Alternative factor structures: A Tucker-structured factor model uses one low-dimensional latent core tensor and mode-specific loading matrices, making it better suited to dimension reduction than mode-specific covariance estimation.

Supplementary Material: Sparse Separable Factor Analysis in the Complex Domain with an Application to Local Field Potential Data

The supplementary material defines the proper complex array normal distribution through vectorization and a Kronecker-product covariance. It also gives the mode-wise matricization and density representation used for array-valued modeling.

  • Array normal distribution: A proper complex array normal variable has a vectorized complex normal distribution with mean vec(M) and covariance V_d ⊗ · · · ⊗ V_1.
  • Array normal distribution: The distribution is denoted by Y ∼C N c1×···×cd (M, V_1, . . . , V_d), with Hermitian positive-definite mode-specific covariance matrices.
  • Mode-wise representation: For mode j, the mode-j matricization uses the product dimension c_−j and the Kronecker covariance V_−j formed from all modes other than j.
  • Mode-wise representation: The array density can be written using the complex matrix-normal density applied to any mode-j matricization.

B Factor analysis for complex-valued vectors

The supplementary factor-analysis derivation begins with an expanded complex working model and then returns to a reduced parameterization. This enables sparsity to be imposed directly on complex loading matrices through complex soft-thresholding.

  • PX-EM parameterization: The PX-EM derivation first estimates the expanded parameters ˜θ = (˜Λ, Γ, Ψ) under the working model.
  • Reduced parameterization: The reduced parameterization sets Λ = ˜ΛR with R = chol(Γ), allowing sparsity to be imposed directly on Λ through a complex soft-thresholding operator.
  • PX-EM parameterization: The working model represents X as ˜Λ˜Z + E, with a complex array-normal latent factor ˜Z and residual E.
  • Complete-data likelihood: The complete-data likelihood factors into f(X | ˜Z; ˜θ)f(˜Z; ˜θ).
  • Complete-data likelihood: Posterior expectations of complete-data sufficient statistics are used in the expanded-parameter likelihood updates.

Proposition B.1 (Posterior expectations of the complete-data sufficient statistics)

The proposition derives posterior expectations of the complete-data sufficient statistics under the expanded complex factor model, then uses them to formulate the PX-EM updates. The reduced loading objective has an unpenalized minimizer given by A(t), while sparsity is imposed through complex soft-thresholding.

  • Posterior expectations: The conditional distribution of the latent factor matrix provides the posterior means and covariances needed for the complete-data sufficient statistics.These statistics enter the expected complete-data objective used in the PX-EM procedure.
  • Posterior expectations: The vectorized working model expresses observations as complex factor scores transformed by the loading matrix plus independent complex-normal residuals.The factor scores and residuals have covariances In ⊗ Γ and In ⊗ Ψ, respectively.
  • Reduced PX-EM objective: At each iteration, fixing Ψ(t) reduces the loading update to minimizing the loading-dependent expected complete-data objective.The positive definiteness of Ψ(t) yields the unique unpenalized minimizer Λ(t+1) = A(t).
  • Sparse update: Elementwise modulus penalties separate across loading entries and produce a unique complex soft-thresholding update.The operator preserves the phase of nonzero entries while reducing their modulus by the threshold.
  • Sparse update: The sparse PX-EM procedure updates the residual diagonal covariance after replacing the unpenalized loading estimate with the sparse estimate.The residual update depends on the sparse loading through Λ(t+1)Λ(t+1)⋆ and retains the resulting sparsity structure.

C Factor analysis for complex-valued arrays

SSFA extends complex factor analysis to array-valued observations by assigning each mode a low-rank loading structure with diagonal residual covariance. A mode-wise PX-EM algorithm whitens the other modes, matricizes the array, and applies sparse factor updates successively.

  • Model: The separable complex array model assigns each mode covariance Σj = ΛjΛj⋆ + Ψj, with complex loadings and positive diagonal residual covariance.The mode-specific parameters are θj = (Λj, Ψj).
  • Model: Stacking observations along a new mode gives an array-valued complex normal model with observation-mode covariance In.This preserves separability across the original modes and the observation mode.
  • Mode-wise algorithm: Each SSFA PX-EM iteration performs d mode-specific updates, holding the remaining modes fixed during each update.The parameter values immediately before and after mode j define the sequential within-iteration updates.
  • Mode-wise algorithm: Before updating mode j, the algorithm whitens every mode except mode j and the observation mode, then matricizes the resulting array.The conditional likelihood becomes equivalent to a vector-valued complex factor model.
  • Mode-wise algorithm: The mode-specific loading update uses complex soft-thresholding, followed by a diagonal residual covariance update based on unexplained covariance.The remaining modes’ parameters stay fixed while the mode-j parameters are updated.

D Identifiability of the mode-specific covariance matrices

The separable covariance factors are identifiable only up to positive multiplicative rescaling, so SSFA imposes a balancing constraint on residual covariances. Equalizing the minimum diagonal entries selects a unique representative without changing the overall Kronecker covariance.

  • Scale nonidentifiability: The covariance of vec(X) is Σd ⊗ · · · ⊗ Σ1, but the individual mode-specific factors are identifiable only up to multiplicative rescaling.The Kronecker product itself remains identifiable under admissible positive rescalings.
  • Balancing constraint: The balancing constraint equates the minimum diagonal entries of the rescaled residual covariance matrices across modes.This is imposed among admissible positive rescalings.
  • Balancing constraint: The balanced tuple of rescaled residual covariance matrices is unique.The unique constants are determined from the inverse minimum diagonal entries, up to the common product constraint.
  • Identifiability result: The corresponding rescaled mode covariances preserve the full Kronecker covariance and form the unique admissible balanced tuple.The uniqueness extends from diagonal residual matrices to the mode-specific covariance matrices.

F.1 Null-space projection error and principal angles

The null-space projection error measures discrepancy between true and estimated loading spaces using rank loss and principal-angle disagreement. Its derivation relates the projector-based error to the geometry of the two column spaces.

  • Error measure: The null-space projection error is proposed as an accuracy measure for estimated loading spaces.The definition compares the column spaces of the true and estimated loading matrices.
  • Error measure: The measure accounts for both the rank of the true and estimated loading matrices and their shared dimensions.It uses rj, brj, and rmin,j = min(rj, brj).
  • Interpretation: The remaining sum measures angular disagreement between the true and estimated spaces over dimensions shared by both.The angles are the principal angles between the two column spaces.
  • Interpretation: When rmin,j = 0, the principal-angle sum is interpreted as zero.This convention handles cases with no shared dimensions.
  • Interpretation: The term rj − rmin,j represents loss from underestimating the true loading-space dimension.This component is separate from angular disagreement over shared dimensions.

F.2 Weighted null-space projection error

The proposition characterizes weighted null-space projection error through angles between the true loading directions and the estimated loading space. The weighting is determined by the squared singular values of the true loading matrix.

  • The weighted null-space projection error equals a weighted average of squared sines measuring angles between true loading vectors and the estimated loading space.The weights are proportional to the squared singular values of the true loading matrix.
  • For each mode j, the true loading matrix Λ_j is compared with its estimate bΛ_j through the estimated column-space projector.The projector is defined using the Moore–Penrose pseudoinverse.
  • A compact singular value decomposition of Λ_j supplies the left singular vectors and positive singular values used in the error decomposition.Only the r_j nonzero singular directions are included.
  • The angle ω_i^(j) is measured between each true left singular vector u_i^(j) and an orthonormal basis for the estimated loading space.The angles lie between 0 and π/2.
  • The derivation uses the orthonormality of the right singular vectors V_j and the unit norm of each left singular vector.
Loading 2608.21551v1…