Covariance estimation for derivatives of functional data using an additive penalty in P-splines

1. Problem Statement

Functional Data Analysis (FDA) relies heavily on Functional Principal Component Analysis (FPCA) to reduce dimensionality and capture dominant modes of variation. While standard FPCA analyzes the functions themselves, many applications (e.g., biomechanics, neuroimaging) require analyzing the derivatives of these functions (e.g., velocity, acceleration) to understand underlying dynamics.

The core challenges addressed in this paper are:

• Derivative Estimation Accuracy: Estimating derivatives from noisy functional data is difficult. Standard smoothing methods often lead to "undersmoothing" when calculating derivatives, resulting in noisy and inaccurate estimates.
• Covariance Derivative Estimation: Derivative-based FPCA requires estimating the covariance of the derivatives (or the derivative of the covariance function). Existing methods (e.g., FPCAder in the fdapace package) rely on local linear smoothers, which can be computationally inefficient and less accurate than modern spline-based approaches.
• Multivariate Extension: There is a lack of efficient methods for Multivariate Functional Principal Component Analysis (MFPCA) specifically tailored to derivatives (DMFPCA) to handle multiple correlated functional trajectories simultaneously.

2. Methodology

The authors propose a unified framework combining P-splines, an additive penalty, and the Fast Covariance Estimation (FACE) algorithm.

A. Covariance Estimation with Additive Penalty

The method extends the FACE algorithm (Xiao et al., 2016) by incorporating an additive penalty into the P-spline smoother.

• Smoothing Matrix: Instead of a single penalty term, the smoother matrix $S$ is constructed using two penalty matrices ($P_1, P_2$) with smoothing parameters $\lambda_1$ and $\lambda_2$.
  $$S = B(B^TB + \lambda_1 P_1 + \lambda_2 P_2)^{-1}B^T$$
• Additive Penalty Logic: By defining a total smoothing level $\lambda_+ = \lambda_1 + \lambda_2$ and a relative weight $w$, the method allows for stronger smoothing in specific dimensions. This is crucial for derivative estimation, as it prevents undersmoothing and ensures the resulting covariance surface is sufficiently smooth.
• Computational Efficiency: The authors utilize an eigen-decomposition of the penalty structure to rewrite the smoother in a "sandwich form" ($S = A_s \Sigma_s A_s^T$). This allows the covariance estimator $\hat{K}$ to be computed efficiently as $\hat{K} = B \Theta B^T$, where $\Theta$ is a small $c \times c$ matrix, avoiding large matrix inversions.

B. Derivative of Covariance

Using the properties of B-splines, the $d$-th order derivative of the covariance estimator is derived analytically.

• If $\hat{K}(s,t) = B^T(s)\Theta B(t)$, then the derivative is:
  $$\frac{\partial^d}{\partial s^d}\frac{\partial^d}{\partial t^d}\hat{K}(s,t) = B_{der}^T(s) \Theta B_{der}(t)$$
• Here, $B_{der}$ is a transformed basis matrix derived from the original B-spline basis and a difference operator. This preserves the "sandwich form," ensuring the derivative covariance estimator remains symmetric and positive semi-definite.

C. Derivative-Based FPCA (DFPCA)

• Eigenfunctions: The eigenfunctions of the derivative covariance are estimated by solving a generalized eigenvalue problem involving $\Theta$ and the integral of the derivative basis functions ($G$). This avoids the need for direct eigendecomposition of the high-dimensional derivative covariance surface.
• Scores (BLUP): The authors propose estimating derivative-based scores using Mixed Model Equations (MME) rather than the traditional Best Linear Unbiased Prediction (BLUP) which requires inverting large matrices. This significantly reduces computational complexity.

D. Derivative Multivariate FPCA (DMFPCA)

The framework is extended to multivariate data (multiple correlated functional trajectories).

• It constructs a score matrix from the univariate derivative scores of each component.
• It performs eigendecomposition on the covariance of these scores to derive Multivariate Derivative Principal Components (DMFPCs) and DMFPC-scores.
• This allows for dimension reduction across multiple variables (e.g., joint angles of knee, hip, and ankle) while capturing their joint dynamic variations.

3. Key Contributions

1. Novel Algorithm (DFPCA): Introduction of the FACE algorithm with an additive penalty for estimating the covariance of derivatives. This improves accuracy over local linear smoothers by enforcing stronger smoothing suitable for derivatives.
2. Computational Efficiency: The method retains the computational advantages of the FACE algorithm (sandwich form) while handling derivatives, making it scalable for large datasets.
3. Score Estimation via MME: A new approach to estimating derivative-based scores using Mixed Model Equations, avoiding expensive matrix inversions.
4. DMFPCA Framework: The first proposed method for Derivative Multivariate Functional Principal Component Analysis, enabling the analysis of joint dynamics in multivariate settings.
5. Theoretical Guarantees: The method ensures the estimated derivative covariance is symmetric and positive semi-definite, a critical property for valid FPCA.

4. Results

Simulation Study

The authors compared their method (DFPCA) against the existing FPCAder() function (based on local linear smoothers) under three scenarios: dense/noise-free, dense/noisy, and sparse/noisy.

• Performance Metrics: Relative Error (RE) for eigenvalues, Integrated Squared Error (ISE) for eigenfunctions, and Mean Squared Error (MSE) for scores.
• Findings:
  • DFPCA consistently outperformed FPCAder across all settings.
  • In the noisy dense case, DFPCA achieved a Mean Relative Integrated Squared Error (RMISE) of 3% compared to 5% for FPCAder.
  • In the sparse noisy case, DFPCA achieved 12% RMISE vs 14% for FPCAder.
  • The method was robust to the choice of penalty orders $(l_1, l_2)$, with $(2,3)$ showing slightly better performance.

Real-World Application: Human Movement Data

The method was applied to a dataset of 120 participants performing three locomotion tasks: level walking, incline walking, and stair ascent.

• Data: Joint angles (ankle, hip, knee) measured via wearable sensors.
• Process: DMFPCA was used to estimate joint angular velocities (derivatives) and extract scores.
• Clustering: K-means clustering on the first two DMFPC-scores successfully distinguished the three locomotion tasks with 96% accuracy (95% CI: 0.91–0.99).
  • Level walking was characterized by earlier directional changes and smaller velocity ranges.
  • Stair walking showed larger velocity magnitudes.
• Significance: The DMFPC-scores provided a compact summary of joint dynamics that was highly effective for classification, outperforming methods that might rely solely on raw angle data.

5. Significance and Conclusion

This paper addresses a critical gap in functional data analysis: the accurate and efficient estimation of derivative-based dynamics.

• Methodological Impact: By integrating additive penalties into the FACE algorithm, the authors provide a superior alternative to local linear smoothers for derivative estimation, particularly in noisy environments.
• Practical Impact: The DMFPCA framework enables researchers to analyze complex, multivariate dynamic systems (like human gait) by focusing on the rate of change rather than just the state. This has direct applications in biomechanics, sports science, and clinical diagnostics (e.g., early detection of movement abnormalities or fall risk).
• Limitations: The authors note that the method is primarily designed for densely observed data (though they adapted PACE for sparse data in simulations) and that boundary effects in P-splines may still slightly reduce accuracy at the edges of the domain.

In summary, the proposed method offers a computationally efficient, statistically robust, and highly accurate framework for analyzing the dynamic behavior of functional data through derivative-based principal component analysis.