A geometric approach to maximum likelihood estimation of the functional principal components from sparse longitudinal data

In this paper, we consider the problem of estimating the eigenvalues and eigenfunctions of the covariance kernel (i.e., the functional principal components) from sparse and irregularly observed longitudinal data. We approach this problem through a maximum likelihood method assuming that the covariance kernel is smooth and finite dimensional. We exploit the smoothness of the eigenfunctions to reduce dimensionality by restricting them to a lower dimensional space of smooth functions. The estimation scheme is developed based on a Newton-Raphson procedure using the fact that the basis coefficients representing the eigenfunctions lie on a Stiefel manifold. We also address the selection of the right number of basis functions, as well as that of the dimension of the covariance kernel by a second order approximation to the leave-one-curve-out cross-validation score that is computationally very efficient. The effectiveness of our procedure is demonstrated by simulation studies and an application to a CD4 counts data set. In the simulation studies, our method performs well on both estimation and model selection. It also outperforms two existing approaches: one based on a local polynomial smoothing of the empirical covariances, and another using an EM algorithm.

💡 Research Summary

This paper tackles the challenging problem of estimating functional principal components (FPCs) when longitudinal data are observed at a few irregular time points per subject—a setting common in many biomedical and social science studies. The authors assume that the underlying covariance kernel of the stochastic process is smooth and of finite rank r, and that its eigenfunctions can be well approximated by a linear combination of a pre‑specified set of smooth basis functions (e.g., B‑splines). Writing the eigenfunctions as ψ(t)=Φ(t)B, where Φ(t) is an M‑dimensional basis vector and B is an M × r coefficient matrix, they impose the orthonormality constraint BᵀB=I_r. This constraint places B on the Stiefel manifold St(r,M), a Riemannian manifold of matrices with orthonormal columns.

Under the standard Gaussian assumptions for the latent scores ξ_iν and measurement errors ε_ij, the log‑likelihood of the observed data can be expressed explicitly in terms of B, the diagonal matrix of eigenvalues Λ=diag(λ₁,…,λ_r), and the error variance σ². The authors derive the Riemannian gradient and Hessian of the log‑likelihood on the Stiefel manifold, enabling a Newton‑Raphson optimization scheme that respects the manifold geometry. Each iteration consists of (i) a Newton step for B followed by a re‑traction onto the manifold (using the QR‑based retraction of Edelman, Arias, and Smith, 1998), and (ii) closed‑form updates for Λ and σ². This approach guarantees that the estimated covariance operator remains positive semi‑definite and that the eigenfunctions stay orthonormal throughout the optimization.

Choosing the rank r and the number of basis functions M is critical. Traditional leave‑one‑curve‑out cross‑validation (LOCO‑CV) is computationally prohibitive because it requires refitting the model n times. The authors propose a second‑order Taylor approximation of the LOCO‑CV score, exploiting the fact that the gradient of the log‑likelihood vanishes at the optimum. The resulting approximate CV criterion can be evaluated in O(n r²) time, making model selection feasible even for moderately large data sets.

Simulation studies cover a range of sparsity levels (2–5 observations per curve) and signal‑to‑noise ratios. The proposed geometric maximum‑likelihood estimator (MLE) consistently outperforms two benchmark methods: (1) an EM algorithm for FPCA (James, Hastie, and Sugar, 2000) and (2) a local polynomial smoothing of empirical covariances (Yao, Müller, and Wang, 2005). The EM method does not enforce orthonormality during iteration, requiring a post‑hoc eigen‑decomposition and often converging slowly; it also fails to satisfy the zero‑gradient condition needed for the CV approximation. The local polynomial approach can produce non‑positive‑definite covariance estimates and even negative error‑variance estimates. In contrast, the geometric MLE yields positive‑definite covariance operators, strictly positive σ², and more accurate eigenvalue/eigenfunction estimates, as measured by mean‑squared error and L² distances.

The methodology is applied to a real HIV data set consisting of CD4+ T‑cell counts measured irregularly over time for 200 patients, with an average of 4.3 measurements per patient. The approximate CV selects r = 2 and M = 7. The first two estimated eigenfunctions capture a clear upward trend and a later decline, respectively, providing clinically interpretable patterns of immune recovery and deterioration. The estimated covariance surface is smooth and respects the required positive semi‑definiteness, unlike the competing methods that either oversmooth or produce artifacts.

Beyond FPCA, the authors discuss how the same geometric framework can be extended to other matrix‑valued parameters subject to orthonormality constraints, such as factor loadings in multivariate mixed‑effects models or rotation matrices in shape analysis. They argue that exploiting the intrinsic Riemannian geometry of the parameter space yields estimators with better statistical efficiency and computational stability.

In summary, the paper makes three major contributions: (1) a reduced‑rank representation of the covariance kernel that leverages smoothness of eigenfunctions; (2) a Newton‑Raphson algorithm on the Stiefel manifold that enforces orthonormality and positive‑definiteness throughout estimation; and (3) a fast, second‑order approximation to LOCO‑CV for simultaneous selection of the number of eigenfunctions and basis dimension. The combination of these ideas provides a powerful, theoretically grounded, and practically efficient solution for functional principal component analysis in the sparsely observed longitudinal data regime.