Discriminant Analysis for Multivariate Time Series

Classifying a multivariate time series (MTS)—such as multi-channel EEG recordings, seismic sensor arrays, or multi-asset financial trajectories—requires distinguishing between classes based on both cross-variable correlations and temporal dynamics. Standard Linear Discriminant Analysis (LDA) and Quadratic Discriminant Analysis (QDA) are designed for static vector observations x∈ℝ^p. When applied directly to an MTS matrix X∈ℝ^p×T (with p variables observed over T time steps), vectorizing the data into a single pT-dimensional vector destroys the chronological ordering and triggers the curse of dimensionality whenever pT exceeds the sample size N.

Extending discriminant analysis to multivariate time series relies on four primary paradigms: Matrix-Variate Discriminant Analysis, Frequency-Domain (Spectral) Discriminant Analysis, Time-Domain Autoregressive (VAR) Discriminant Analysis, and Functional Linear Discriminant Analysis.

---

1. Why Vectorized LDA Fails on Multivariate Time Series

Suppose we observe N multivariate time series X₁, …, X_N, each belonging to one of G classes π₁, …, π_G. Each observation X_i is a p×T matrix:

X_i=[x_i,1, x_i,2, …, x_i,T]∈ℝ^p×T

If you flatten X_i into vec(X_i)∈ℝ^pT, classical LDA must estimate a pooled covariance matrix Σ of size pT×pT, containing ½pT(pT+1) free parameters:

• Singularity (pT≫N): Even for a modest 10-channel sensor recorded over 200 time steps (pT=2000), Σ requires estimating over 2 million parameters. Unless N>2000, the sample covariance matrix is singular and non-invertible.
• Loss of Kronecker Structure: Vectorization treats the covariance between channel j at time t and channel k at time s as an unstructured parameter, ignoring that temporal autocorrelation and spatial cross-channel correlation can be factored separately.
• Phase Sensitivity: Slight time shifts or warping across realizations misalign the vectorized features, degrading linear separability.

---

2. Core Paradigms of MTS Discriminant Analysis

Method	Underlying Assumption	Covariance / Parameter Complexity	Best Use Case
Bilinear / 2D-LDA (Kronecker)	Separable spatial (p×p) and temporal (T×T) covariance	O(p2+T2) instead of O(p2T2)	Short, phase-aligned multi-sensor trials (ERP/EEG, gesture recognition)
Spectral (Whittle) Discriminant Analysis	Stationary Gaussian processes with distinct spectral density matrices	p×p Hermitian matrix per frequency ωk	Long stationary signals (seismology, acoustics, rotating machinery)
Parametric VAR Discriminant Analysis	Processes generated by class-specific Vector Autoregressive models	dp2+½p(p+1) per class (for order d)	Econometrics, physiological feedback systems
Functional LDA (fLDA)	Smooth continuous trajectories projected onto basis functions	O(p2K2) where K≪T is the number of basis functions	Irregularly sampled or smooth longitudinal trajectories

Paradigm A: Matrix-Variate and Bilinear Discriminant Analysis (2D-LDA)

Instead of flattening X∈ℝ^p×T, Matrix-Variate Discriminant Analysis models X using a Matrix Normal Distribution ℳ풩_p,T(M_g, Σ_S, Σ_T), where the overall pT×pT covariance matrix decomposes via a Kronecker product:

Cov(vec(X))=Σ_T⊗Σ_S

Here, Σ_S∈ℝ^p×p captures the spatial (cross-channel) covariance, and Σ_T∈ℝ^T×T captures the temporal covariance (often constrained to an AR(1) Toeplitz structure Σ_T,(t,s)=ρ^∣t−s∣).

The Matrix-Variate Mahalanobis Distance

Under equal Kronecker covariances across classes (Σ_S, Σ_T), the linear discriminant score for assigning X to class g with mean matrix M_g∈ℝ^p×T and prior p_g is:

δ_g(X)=−½tr(Σ_S⁻¹(X−M_g)Σ_T⁻¹(X−M_g)^⊤)+log p_g

where tr(⋅) denotes the matrix trace.

Two-Directional LDA (2D-LDA)

When dimension reduction is the primary goal, 2D-LDA seeks two projection matrices—a spatial projection L∈ℝ^p×d₁ and a temporal projection R∈ℝ^T×d₂—that map X to a compact core matrix Y=L^⊤XR∈ℝ^d₁×d₂ by maximizing the generalized Fisher criterion across between-class (S_B) and within-class (S_W) scatter matrices iteratively via alternating least squares.

Paradigm B: Frequency-Domain (Spectral) Discriminant Analysis

When time series are long and stationary, class differences rarely lie in the time-aligned mean M_g (which is often zero for all classes), but rather in the autocovariance structure Γ_g(h)=Cov_g(x_t, x_t−h). Pioneered by Shumway and Kakizawa, Spectral Discriminant Analysis transforms the MTS into the frequency domain via the Discrete Fourier Transform (DFT):

J(ω_k)=(1/√T)∑_t=1^Tx_te^−i2πω_kt,   ω_k=k/T

For a stationary multivariate process, the DFT vectors J(ω_k)∈ℂ^p at distinct Fourier frequencies ω_k are asymptotically independent complex normal random variables with zero mean and covariance equal to the p×p spectral density matrix f_g(ω_k):

f_g(ω)=∑_h=−∞^∞Γ_g(h)e^−i2πωh

Whittle Likelihood & Kullback-Leibler Discrimination

Because frequencies decouple asymptotically, the high-dimensional pT×pT time-domain covariance problem reduces to T/2 independent p×p Hermitian matrix operations. Using the Whittle log-likelihood or the Chernoff / Kullback-Leibler spectral divergence, a new time series with periodogram I(ω_k)=J(ω_k)J^*(ω_k) is assigned to the class g that minimizes the spectral discriminant criterion:

D_KL(I‖f_g)=(1/M)∑_k=1^M[tr(I(ω_k)f_g⁻¹(ω_k))+log det f_g(ω_k)]−log p_g

This is mathematically equivalent to Quadratic Discriminant Analysis (QDA) operating on the frequency-domain covariance matrices.

Paradigm C: Time-Domain Parametric (VAR) Discriminant Analysis

When the data-generating mechanism is well-modeled by a Vector Autoregressive process of order d (VAR(d)), each class g is parameterized by its transition matrices {Φ_g,1, …, Φ_g,d} and white-noise innovation covariance Σ_ε,g:

x_t=μ_g+∑_ℓ=1^dΦ_g,ℓx_t−ℓ+ε_g,t,   ε_g,t∼풩(0, Σ_ε,g)

To classify an unlabeled multivariate series X, you pass it through each class’s estimated VAR(d) filter to compute the one-step-ahead prediction residuals ε̂_g,t and assign X to the class maximizing the conditional Gaussian log-likelihood:

ℒ_g(X)=−((T−d)/2)log det(Σ̂_ε,g)−½∑_t=d+1^Tε̂_g,t^⊤Σ̂_ε,g⁻¹ε̂_g,t+log p_g

Paradigm D: Functional Linear Discriminant Analysis (fLDA)

When multivariate trajectories are smooth and potentially sampled at irregular time points t_i,1, …, t_i,Ti, Functional Data Analysis (FDA) treats each sample as a p-dimensional vector of continuous functions 풳_i(t)=[X_i,1(t), …, X_i,p(t)]^⊤ on t∈[0, 풯].

1. Basis Expansion: Project each channel onto a truncated basis of K functions (such as B-splines, wavelets, or Fourier bases) φ(t)=[φ₁(t), …, φ_K(t)]^⊤, yielding a coefficient matrix C_i∈ℝ^p×K.
2. Multivariate Functional Principal Component Analysis (MFPCA): Extract the top M eigenfunctions ψ_m(t) of the cross-covariance operator to regularize the infinite-dimensional space.
3. Regularized Fisher Projection: Perform penalized LDA directly on the basis coefficients or MFPCA scores using a roughness penalty λ∫[D²β(t)]²dt to guarantee smooth discriminant weight functions β(t).

---

3. Selecting the Right Approach in Practice

1. Do classes differ in their mean wave shape or in their oscillation/correlation structure?
   • If classes differ by their time-locked average trajectory (e.g., an evoked brain response after a stimulus), use Matrix-Variate LDA / 2D-LDA or Functional LDA.
   • If all classes have a zero mean and differ by frequency rhythms, volatility, or cross-channel lead-lag relationships (e.g., seizure vs. baseline EEG, fault detection in turbines), mean-based LDA will fail completely; use Spectral (QDA) Discriminant Analysis or VAR Discriminant Analysis.
2. Regularization is Mandatory: Even after Kronecker or spectral reduction, when p is large, shrink the spatial or spectral covariance matrices using Ledoit-Wolf shrinkage (Σ̂_shrunk=(1−γ)Σ̂+γτI_p) or ℓ₁-penalized sparse discriminant analysis to prevent overfitting on small sample sizes.