  Covariance estimation is essential yet underdeveloped for analyzing
multivariate functional data. We propose a fast covariance estimation method
for multivariate sparse functional data using bivariate penalized splines. The
tensor-product B-spline formulation of the proposed method enables a simple
spectral decomposition of the associated covariance operator and explicit
expressions of the resulting eigenfunctions as linear combinations of B-spline
bases, thereby dramatically facilitating subsequent principal component
analysis. We derive a fast algorithm for selecting the smoothing parameters in
covariance smoothing using leave-one-subject-out cross-validation. The method
is evaluated with extensive numerical studies and applied to an Alzheimer's
disease study with multiple longitudinal outcomes.
