Fast covariance estimation for high-dimensional functional data

We propose two fast covariance smoothing methods and associated software that scale up linearly with the number of observations per function. Most available methods and software cannot smooth covariance matrices of dimension \(J>500\); a recently introduced sandwich smoother is an exception but is not adapted to smooth covariance matrices of large dimensions, such as \(J= 10{,}000\). We introduce two new methods that circumvent those problems: (1) a fast implementation of the sandwich smoother for covariance smoothing; and (2) a two-step procedure that first obtains the singular value decomposition of the data matrix and then smoothes the eigenvectors. These new approaches are at least an order of magnitude faster in high dimensions and drastically reduce computer memory requirements. The new approaches provide instantaneous (a few seconds) smoothing for matrices of dimension \(J=10{,}000\) and very fast (\(<\)10 min) smoothing for \(J=100{,}000\). R functions, simulations, and data analysis provide ready to use, reproducible, and scalable tools for practical data analysis of noisy high-dimensional functional data.

This work was supported by Grant Number R01EB012547 from the National Institute of Biomedical Imaging And Bioengineering and Grant Number R01NS060910 from the National Institute of Neurological Disorders and Stroke. This work represents the opinions of the researchers and not necessarily that of the granting organizations.

Authors and Affiliations

1. Department of Biostatistics, Johns Hopkins University, Baltimore, MD, USA

   Luo Xiao, Vadim Zipunnikov & Ciprian Crainiceanu

2. Department of Statistical Science and School of Operations Research and Information Engineering, Cornell University, Ithaca, NY, USA

   David Ruppert

Corresponding author

Correspondence to Luo Xiao.

Appendix 1: Proofs

Proof of Proposition 1

The design matrix \(\mathbf{B}\) is of full rank (Xiao et al. 2012). Hence \(\mathbf{B}^T\mathbf{B}\) is invertible and \(\mathbf{A}_S\) is of rank \(c\). \({\varvec{\Sigma }}_S\) is a diagonal matrix with all elements greater than 0 and \({\widetilde{\mathbf{Y}}}\) is of rank at most \(\min (c,I)\). Hence \({\widetilde{\mathbf{K}}}= \mathbf{A}_S\left( I^{-1}{\varvec{\Sigma }}_S{\widetilde{\mathbf{Y}}}{\widetilde{\mathbf{Y}}}^T{\varvec{\Sigma }}_S\right) \mathbf{A}_S^T\) has a rank at most \(\min (c,I)\) and the proposition follows. \(\square \)

Proof of Proposition 2

First of all, \(\text{ tr }(\mathbf {S}) = \text{ tr }({\varvec{\Sigma }}_S)\) which is easy to calculate. We now compute \(\sum _{i=1}^I \Vert \mathbf{Y}_i - \mathbf{S}\mathbf{Y}_i\Vert ^2\). Because \( \Vert \mathbf{Y}_i - \mathbf{S}\mathbf{Y}_i\Vert ^2 = \mathbf{Y}_i^T(\mathbf{S}-\mathbf{I}_J)^2\mathbf{Y}_i = \text{ tr }\{(\mathbf{S}-\mathbf{I}_J)^2\mathbf{Y}_i\mathbf{Y}_i^T\}\),

It can be shown that \(\mathbf{S}^2 = \mathbf{A}_S {\varvec{\Sigma }}_S^2\mathbf{A}_S^T\). Hence \(\text{ tr }(\mathbf{S}^2\mathbf{Y}\mathbf{Y}^T) = \text{ tr }(\mathbf{Y}^T\mathbf{S}^2\mathbf{Y}) = \text{ tr }({\widetilde{\mathbf{Y}}}^T{\varvec{\Sigma }}_S^2{\widetilde{\mathbf{Y}}})=\text{ tr }({\varvec{\Sigma }}_S^2{\widetilde{\mathbf{Y}}}{\widetilde{\mathbf{Y}}}^T)\). Similarly, we derive \(\text{ tr }(\mathbf{S}\mathbf{Y}\mathbf{Y}^T) = \text{ tr }({\varvec{\Sigma }}_S{\widetilde{\mathbf{Y}}}{\widetilde{\mathbf{Y}}}^T)\). We have \(\text{ tr }(\mathbf{Y}\mathbf{Y}^T) = \Vert \mathbf{Y}\Vert _F^2\). It follows that

\(\square \)

Proposition 3

The computation time of FACE is \(O(IJc +Jc^2+c^3 + ck_0)\), where \(k_0\) is the number of iterations needed for selecting the smoothing parameter (see Sect. 3.2), and the total required computer memory is \(O(JI+I^2+Jc+c^2+k_0)\) memory units.

Proof of Proposition 3

We need to compute or store the following quantities: \(\mathbf{X}\), \(\mathbf{B}\), \(\mathbf{B}^T\mathbf{B}\), \((\mathbf{B}^T\mathbf{B})^{-1/2}, \mathbf{P}, (\mathbf{B}^T\mathbf{B})^{-1/2}\mathbf{P}(\mathbf{B}^T\mathbf{B})^{-1/2}, \mathbf{A}_S, {\widetilde{\mathbf{Y}}}, \mathbf{A}, \mathbf{U}\), and \(\mathbf{A}_S \mathbf{A}\). For the computational complexity, \(\mathbf{B}^T\mathbf{B}\), \(\mathbf{A}_S = \mathbf{B}(\mathbf{B}^T\mathbf{B})^{-1/2}\mathbf{U}\), and \(\mathbf{A}_S\mathbf{A}\) require \(O(Jc^2)\) computations; \((\mathbf{B}^T\mathbf{B})^{-1/2}\), \(\mathbf{P}\), \((\mathbf{B}^T\mathbf{B})^{-1/2}\mathbf{P}(\mathbf{B}^T\mathbf{B})^{-1/2}\), \(\mathbf{A}\), and \(\mathbf{U}\) require \(O(c^3)\) computations; \({\widetilde{\mathbf{Y}}}= \mathbf{A}_S^T\mathbf{Y}\) requires \(O(JIc)\) computations. So in total, \(O(JIc+Jc^2+c^3)\) computations are required. For the memory burden, the loading of \(\mathbf{Y}\) requires \(O(JI)\) memory units, computer of \(\mathbf{B}\) and \(\mathbf{A}_S\mathbf{A}\) requires \(O(Jc)\) memory units, and other objects require \(O(c^2)\) memory units. \(\square \)

Proof of Theorem 1 We have \({\widehat{{\varvec{\xi }}}}_i \!=\! J^{-1/2}(\mathbf{A}_S{\widehat{\mathbf{A}}}_N)^T\mathbf{Y}_i\! =\! J^{-1/2}{\widehat{\mathbf{A}}}_N^T(\mathbf{A}_S^T\mathbf{Y}_i)=J^{-1/2}{\widehat{\mathbf{A}}}_N^T{\widetilde{\mathbf{Y}}}_i\).\(\square \)

Proof of Theorem 2 Let \({\widetilde{\mathbf{A}}}_N\) denote the first \(N\) columns of \(\mathbf{A}_S\mathbf{A}\), then \({\widetilde{\mathbf{A}}}_N = \mathbf{A}_S{\widehat{\mathbf{A}}}\). The estimated BLUPs for \({\varvec{\xi }}_i\) (Ruppert et al. 2003) is

The inverse matrix in the above equality can be replaced by the following [Seber (2007), page 309, equality b(i)],

It follows that

\(\square \)

Appendix 2: Empirical covariance operators for \(K_X\) and \(K_U\)

Let \(I\) denote the number of pairs of cases and controls. For simplicity, we assume estimates of \(\mu _A(t)\) and \(\mu _C(t)\) have been subtracted from \(Y_{iA}\) and \(Y_{iC}\), respectively. Let \(\mathbf{Y}_{iA} = (Y_{iA}(t_1),\ldots , Y_{iA}(t_T))^T\) and \(\mathbf{Y}_{iC} = (Y_{iC}(t_1),\ldots , Y_{iC}(t_J))^T\). By Zipunnikov et al. (2011), we have estimates of the covariance operators,

and

Let \(\mathbf{Y}_A = [\mathbf{Y}_{1A},\ldots , \mathbf{Y}_{nA}]\), \(\mathbf{Y}_C = [\mathbf{Y}_{1C}, \ldots , \mathbf{Y}_{nC}]\) and \(\mathbf{Y}= [\mathbf{Y}_A,\mathbf{Y}_C]\). Then \(\mathbf{Y}\) is of dimension \(J\times 2I\). It can be shown that \({\widehat{\mathbf{K}}}_X = \mathbf{Y}\mathbf{H}_X\mathbf{Y}^T\) and \({\widehat{\mathbf{K}}}_U = \mathbf{Y}\mathbf{H}_U\mathbf{Y}^T\), where

Reprints and permissions