# Functional Principal Subspace Sampling for Large Scale Functional Data
Analysis

Shiyuan<EMAIL_ADDRESS>[ Xiaomeng
<EMAIL_ADDRESS>[ Center for Applied Statistics,
Institute of Statistics and Big Data, Renmin University of China, Bejing,
100872, China.  Department of Statistics, Texas A&M University, College
Station, TX, 77840, USA.

(2020)

###### Abstract

Functional data analysis (FDA) methods have computational and theoretical
appeals for some high dimensional data, but lack the scalability to modern
large sample datasets. To tackle the challenge, we develop randomized
algorithms for two important FDA methods: functional principal component
analysis (FPCA) and functional linear regression (FLR) with scalar response.
The two methods are connected as they both rely on the accurate estimation of
functional principal subspace. The proposed algorithms draw subsamples from
the large dataset at hand and apply FPCA or FLR over the subsamples to reduce
the computational cost. To effectively preserve subspace information in the
subsamples, we propose a functional principal subspace sampling probability,
which removes the eigenvalue scale effect inside the functional principal
subspace and properly weights the residual. Based on the operator perturbation
analysis, we show the proposed probability has precise control over the first
order error of the subspace projection operator and can be interpreted as an
importance sampling for functional subspace estimation. Moreover,
concentration bounds for the proposed algorithms are established to reflect
the low intrinsic dimension nature of functional data in an infinite
dimensional space. The effectiveness of the proposed algorithms is
demonstrated upon synthetic and real datasets.

62R10,

62H25,

62D99,

functional principal component analysis,

functional linear regression,

randomized algorithm,

operator perturbation theory,

###### keywords:

[class=MSC]

###### keywords:

††volume: 0††issue: 0

2010.00000

and

###### Contents

  1. 1 Introduction
    1. 1.1 Review of the Importance Sampling
    2. 1.2 Limitation of the Importance Sampling
    3. 1.3 Contributions of This Work
  2. 2 Functional Principal Subspace Sampling for FPCA
    1. 2.1 Randomized FPCA
    2. 2.2 Functional Principal Subspace Sampling Probability
    3. 2.3 Implementation of FunPrinSS
  3. 3 Theoretical Analysis of Randomized FPCA with FunPrinSS
    1. 3.1 Operator Perturbation Theory
    2. 3.2 Main Results
    3. 3.3 Theoretical Justification of Algorithm 3
    4. 3.4 Choice of $R$
  4. 4 Extenstion: Randomized Functional Linear Regression
    1. 4.1 Randomized FLR
    2. 4.2 Analysis of Randomized FLR with FunPrinSS
  5. 5 Simulation Study
    1. 5.1 Data Generation
    2. 5.2 Randomized FPCA
    3. 5.3 Randomized FLR
  6. 6 Real Data Analysis
    1. 6.1 Results for Randomized FPCA
    2. 6.2 Results for Randomized FLR
  7. A Comparison with sampling probabilities in the multivariate context
  8. B Concentration Inequality for Compact Operators
  9. C Proof for Section 3
    1. C.1 Proof of Lemma 6
    2. C.2 Proof of Theorem 8
    3. C.3 Proof of Lemma 17
    4. C.4 Proof of Lemma 18
  10. D Proof for Section 4
    1. D.1 Proof of Lemma 13
    2. D.2 Proof of Theorem 14
    3. D.3 Proof of Lemma 22
    4. D.4 Proof of Lemma 20
    5. D.5 Proof of Lemma 21
  11. E Proof of Proposition 9
  12. F Additional Results for Simulation Study

##  1 Introduction

The modern era brings us the opportunities and challenges of big data, where
both sample size and dimension increase out of the capacity of classical
statistical methods. Examples of the large scale datasets include: millions of
spectra collected by astronomical surveys (e.g. LAMOST of Zhao et al. (2012)
and SDSS of Eisenstein et al. (2011)); petabyte amounts of hyperspectral
images recorded by airborne or satellite remote sensing (Liu et al., 2018);
brain medical images through the functional magnetic resonance imaging (fMRI
Tian, 2010; Turkbrowne, 2013). A common characteristic of these datasets is
that their observations can be viewed as functions over some continuous domain
(e.g. temporal or spatial domain), and functional data analysis (FDA, Kokoszka
and Reimherr, 2017) can serve as a designated toolbox. FDA treats functional
observations as realizations of some stochastic process with low intrinsic
dimension in an infinite dimensional Hilbert space. This means observations
can be effectively approximated by splines, wavelets or linear combinations of
a few functional principal components. In contrast to classical multivariate
methods, FDA takes these features into consideration to achieve optimal
convergence rates. Despite their demonstrated value in theory and application,
the current FDA methods lack the scalability for modern large scale datasets
(Hadjipantelis and Müller, 2018).

Among the current approaches to making statistical methods scalable to large
datasets, randomized algorithms have recently gained popularity. The
randomized algorithms approach the problem by constructing a substantially
smaller sketch of the large scale data at hand. Then, existing methods are
applied to the sketch to reduce computation cost. When the sketch keeps the
most relevant information, the result computed from the sketch should remain
close to the result from the original dataset. There are a few ways to
construct a sketch of a large data set. One approach is to draw subsamples
with respect to some carefully designed probability, which will select
informative samples with larger probability. These sampling probabilities
include the importance sampling in matrix multiplication (Drineas, Kannan and
Mahoney, 2006a), the leverage sampling for least squares regression (Drineas
et al., 2012; Ma, Mahoney and Yu, 2015), the subspace sampling for low rank
matrix construction (Drineas, Mahoney and Muthukrishnan, 2006), etc. Effective
sampling probabilities have also been proposed for logistic regression (Wang,
Zhu and Ma, 2018; Wang, 2019) and generalized linear models (Zhang, Ning and
Ruppert, 2021) to minimize the asymptotic variance of the estimator. An
alternative way towards sketching is to mix the original data with a random
projection and draw samples from the projected data (Drineas et al., 2011;
Wang, Gittens and Mahoney, 2017).

Since the seminal works of Drineas, Kannan and Mahoney (2006a, b, c) for
matrix multiplication and approximation, the idea of randomized algorithm has
been successfully applied to optimization (Pilanci and Wainwright, 2015,
2017), low rank matrix estimation (Halko, Martinsson and Tropp, 2011; Drineas,
Mahoney and Muthukrishnan, 2006), least squares regression (Drineas et al.,
2012), nonparametric kernel regression (Yang et al., 2017), etc. These
algorithms are able to yield comparatively accurate results at reduced
computational and storage costs. See the references Drineas and Mahoney
(2018); Woodruff (2014) for an overview. Most theoretical analysis of these
randomized algorithms is conducted from the algorithmic perspective, where the
analysis is carried out conditionally on an arbitrarily fixed dataset. Some
recent works (Ma, Mahoney and Yu, 2015; Raskutti and Mahoney, 2016; Wang,
Gittens and Mahoney, 2017) also draw analysis from the statistical
perspective, where statistical properties such as bias and average prediction
error are considered.

The current literature of randomized algorithms focuses on multivariate
statistical methods which are not directly applicable to the functional data
setting. Their theoretical results also does not naturally extend to
functional data in an infinite dimensional space. He and Yan (2020) recently
studied randomized estimation of functional covariance operator. Still, most
functional data methods remain to be explored via the randomized algorithm
approach for scalability. Among these methods, functional principal component
analysis (FPCA) (e.g. James, Hastie and Sugar, 2000; Yao, Müller and Wang,
2005; Peng and Paul, 2009) is one of the most important. Like the classical
principal component analysis, FPCA identifies a few orthonormal functional
principal components (FPC) explaining most of the variability of a dataset.
The computation and theoretical analysis of functional data in an infinite-
dimensional space are usually facilitated by FPCA. Indeed, a large amount of
FDA methods have been proposed based on FPC. For example, when fitting
functional linear regression, (Hall et al., 2007; Kato, 2012) first identify a
few FPCs of the predictor, and then the model is fitted via the leading FPCs.
For functional two-sample test, the test statistic can be constructed based on
FPCs (Horváth, Kokoszka and Reeder, 2013). The work of Delaigle et al. (2010)
relies on FPCs to define probability density for random functions. More
examples can be found in Horváth and Kokoszka (2012).

Due to the fundamental role of FPCA, we study a randomized FPCA algorithm and
extend it to functional linear regression (FLR) with scalar response in this
work. Our work focuses on the fully observed functional data case, where the
function value is known at any point of its domain. In practical applications,
for example, these could be functions recorded with a fixed high frequency
over a time interval or recorded over a dense spatial grid. This type of data
has been widely collected by astronomical spectral surveys, remote sensing,
etc, as discussed at the beginning of this section. Specifically, we will view
functional observations as indivisible elements in a Hilbert space, and the
sketched data is constructed based on subsampling. Our basic algorithm is
motivated by the work of He and Yan (2020) for which we now provide a review.

###  1.1 Review of the Importance Sampling

Let $\mathscr{H}_{X}$ be a Hilbert space equipped with an inner product
$\langle\cdot,\cdot\rangle$ and an induced norm $\|\cdot\|$. For example, when
the space $\mathscr{H}_{X}=L_{2}(T)$ is the set of all square integrable
functions over a compact interval $T$, the inner product is $\langle
x,x^{\prime}\rangle=\int_{T}x(t)x^{\prime}(t)\,\mathrm{d}t$ and the norm is
$\|x\|=\big{(}\int_{T}x^{2}(t)\,\mathrm{d}t\big{)}^{1/2}$ for
$x,x^{\prime}\in\mathscr{H}_{X}$. Given a random element $x\in\mathscr{H}_{X}$
with zero mean and finite second moment $\mathbb{E}\|x\|^{2}$, its covariance
operator $\mathscr{C}_{XX}$: $\mathscr{H}_{X}\mapsto\mathscr{H}_{X}$ is a
mapping such that $u\in\mathscr{H}_{X}$ is mapped to $\mathbb{E}\langle
u,x\rangle x$. Equivalently, we can write
$\mathscr{C}_{XX}=\mathbb{E}\left(x\otimes x\right)$, where $\otimes$ is
tensor product such that $\left(x\otimes x\right)u=\langle x,u\rangle\times x$
for any $u\in\mathscr{H}_{X}$. Covariance operator generalizes the concept of
covariance matrix in multivariate statistics. For any $u,v\in\mathscr{H}_{X}$,
it holds that $\langle\mathscr{C}_{XX}u,v\rangle=\mathbb{E}(\langle
x,u\rangle\langle x,v\rangle)$. We can see that
$\langle\mathscr{C}_{XX}u,v\rangle$ quantifies the correlation between the
random scalars $\langle x,u\rangle$ and $\langle x,v\rangle$. A comprehensive
introduction to these FDA concepts can be found in Hsing and Eubank (2015).

Input: Dataset $\\{x_{n}\\}_{n=1}^{N}$; sampling probability
$\\{p_{n}\\}_{n=1}^{N}$; subsample size $C$.

Output: $\widetilde{\mathscr{C}}_{XX}$.

1: for $c=1,\cdots,C$ do

2: Sample $\tilde{x}_{c}$ from $\\{x_{n}\\}_{n=1}^{N}$ and get $\tilde{p}_{c}$
with probability
$\mathbb{P}\big{(}(\tilde{x}_{c},\tilde{p}_{c})=(x_{n},p_{n})\big{)}=p_{n}.$

3: end for

4: Compute
$\widetilde{\mathscr{C}}_{XX}=\frac{1}{C}\sum_{c=1}^{C}\frac{1}{N\tilde{p}_{c}}\tilde{x}_{c}\otimes\tilde{x}_{c}$.

Algorithm 1 Randomized Covariance Operator Estimation

Suppose we have $N$ (not necessarily independent) realizations
$\\{x_{n}\\}_{n=1}^{N}$ of the random element $x$. The covariance operator
$\mathscr{C}_{XX}$ can be estimated by the empirical covariance operator
$\widehat{\mathscr{C}}_{XX}=(1/N)\sum_{n=1}^{N}x_{n}\otimes x_{n}$. The
empirical covariance operator involves the calculation of $N$ tensor products
and their summation. This computational cost is overwhelming for large sample
size. For example, when each $x_{n}\in L_{2}(T)$ is digitally stored as a high
dimensional vector of length $L$ over a dense grid of $T$, the cost of
computing $\widehat{\mathscr{C}}_{XX}$ scales as $O(NL^{2})$. In a typical
astronomical survey, $L$ is in the order of thousands, and $N$ has the
magnitude of millions. Under this scenario, computing a covariance operator
using all available data for scientific research is unrealistic.

Based on the randomized matrix multiplication of Drineas, Kannan and Mahoney
(2006a), the algorithm of He and Yan (2020) estimates the covariance operator
from a sketch of much smaller size. The procedure is summarized in Algorithm
1. Given a full dataset $\\{x_{n}\\}_{n=1}^{N}$ and a sampling probability
$\\{p_{n}\\}_{n=1}^{N}$, it draws with replacement a subsample
$\\{\tilde{x}_{c}\\}_{c=1}^{C}$ of size $C$ from $\\{x_{n}\\}_{n=1}^{N}$.
Then, the subsampled empirical covariance operator is computed as

$\widetilde{\mathscr{C}}_{XX}=\frac{1}{C}\sum_{c=1}^{C}\frac{1}{N\tilde{p}_{c}}\tilde{x}_{c}\otimes\tilde{x}_{c}.$
(1.1)

Suppose $c$-th subsample is indexed by $i_{c}$ in the original dataset, we
have set $\tilde{x}_{c}=x_{i_{c}}$ and $\tilde{p}_{c}=p_{i_{c}}$
correspondingly in the above expression. It is obvious that, conditional on
the full dataset $\\{x_{n}\\}_{n=1}^{N}$, the subsampled
$\widetilde{\mathscr{C}}_{XX}$ is an unbiased estimator of the full sample
covariance operator $\widehat{\mathscr{C}}_{XX}$ with any strictly positive
probability $\\{p_{n}\\}_{n=1}^{N}$. It is shown in He and Yan (2020) that the
optimal sampling probability, which minimizes the expected squared Hilbert-
Schmidt norm of the subsampling error
$\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}$, is of the form
$p_{n}\propto\|x_{n}\|^{2}$. This sampling probability will be referred as the
importance sampling (IMPO) throughout this work.

###  1.2 Limitation of the Importance Sampling

Recall the purpose of this work is to develop randomized functional principal
component analysis (FPCA) and its extension to functional linear regression
(FLR). At first glance, the importance sampling (IMPO) of He and Yan (2020) is
already a plausible sampling strategy. By direct application of perturbation
theory, the eigenfunctions (functional principal components) of
$\widetilde{\mathscr{C}}_{XX}$ and $\widehat{\mathscr{C}}_{XX}$ should not be
far away from each other, provided the related eigengap is positive and the
difference $\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}$ is small.
However, for FPCA and FLR, this importance sampling probability is far from
being optimal.

To see the reasons, consider a toy principal component analysis example in
$\mathbb{R}^{4}$. We generate $N=1000$ observations by
$\mathbf{x}_{n}=\sigma_{1}\xi_{n1}\boldsymbol{\theta}_{1}+\sigma_{2}\xi_{n2}\boldsymbol{\theta}_{2}+\sigma_{3}\xi_{n3}\boldsymbol{\theta}_{3}$,
where the principal component directions
$\boldsymbol{\theta}_{1}=(1,1,0,0)^{T}/\sqrt{2}$,
$\boldsymbol{\theta}_{2}=(1,-1,0,0)^{T}/\sqrt{2}$ and
$\boldsymbol{\theta}_{3}=(0,0,1,1)^{T}/\sqrt{2}$ are scaled by
$\sigma_{1}=10$, $\sigma_{2}=2$, $\sigma_{3}=0.1$, respectively. The scores
$\xi_{n1},\xi_{n2},\xi_{n3}\in\mathbb{R}$ are drawn from the standard normal
distribution. Figure 1 plots the generated points by their first two
coordinates. Suppose we want to estimate the first two principal components
$\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}$ from subsamples, and the
sampling is taken with respect to the IMPO sampling probability (He and Yan,
2020). To illustrate the effect of IMPO, we color points with larger
$p_{n}\propto\|x_{n}\|^{2}=\sigma_{1}^{2}\xi_{n1}^{2}+\sigma_{2}^{2}\xi_{n2}^{2}+\sigma_{3}^{2}\xi_{n3}^{2}$
by darker blue in the left penal of Figure 1. It is evident this probability
over-emphasizes the observations along the first direction
$\boldsymbol{\theta}_{1}$ while neglects the observations along
$\boldsymbol{\theta}_{2}$. When $C$ subsamples are drawn according to IMPO,
the resulting subsamples will likely contain abundant information for
estimating $\boldsymbol{\theta}_{1}$, but fail to determine
$\boldsymbol{\theta}_{2}$ accurately.

In the above example, the first two scaling factors only have a moderately
large ratio $\sigma_{1}^{2}/\sigma_{2}^{2}=25$, but their effect in
determining $p_{n}$ is not negligible. To resolve the issue, we may consider
to remove the scaling factor to get a new sampling probability
$p_{n}\propto\xi_{n1}^{2}+\xi_{n2}^{2}$. The points in the right penal of
Figure 1 are colored accordingly. Notice this probability puts equal sampling
weights for observations along $\boldsymbol{\theta}_{1}$ and
$\boldsymbol{\theta}_{2}$. The subsamples, which are drawn by this
probability, are more likely to contain balanced information for estimating
the first two principal components. Moreover, this idea can be generalized to
infinite dimensional space to estimate the leading $R(>0)$ functional
principal components.

Figure 1: An toy example for principal component analysis in $\mathbb{R}^{4}$.
The $N=1000$ points
$\mathbf{x}_{n}=(x_{n1},x_{n2},x_{n3},x_{n4})^{T}\in\mathbb{R}^{2}$ are
plotted by their first two coordinates $(x_{n1},x_{n2})$.. The red lines mark
the first and second principal components $\boldsymbol{\theta}_{1}$ and
$\boldsymbol{\theta}_{2}$. The left panel colors the data points by
$p_{n}\propto\|x_{n}\|^{2}$, while the right panel by
$p_{n}\propto(\xi_{n1}^{2}+\xi_{n2}^{2})$. Points with larger sampling
probability has darker blue color.

###  1.3 Contributions of This Work

This work develops randomized algorithms for functional principal component
analysis (FPCA) and functional linear regression (FLR) with scalar response.
Our methodological and theoretical contributions are summarized in the
following paragraphs.

In Section 2, a randomized algorithm is developed for FPCA. Motivated by the
example from Section 1.2, we propose the functional principal subspace
sampling probability (FunPrinSS) to effectively preserve subspace information
in subsamples. The FunPinSS probability removes the eigenvalue scale effect
inside the functional principal subspace; and at the same time, it properly
weights the complementary subspace residual. The latter part is necessary
because we need to control the interaction between the leading and the
remaining eigenfunctions in the infinite dimensional space. As the exact value
of the FunPrinSS probability is unknown unless we execute full sample
computation, a fast and theoretically justifiable algorithm is proposed for
practical implementation.

For the theoretical analysis of the randomized FPCA algorithm, the task of
precisely estimating the leading $R$ eigenfunctions is translated into
estimating the projection operator of their spanned subspace. The analysis is
conducted from an algorithmic perspective, which is conditional on the full
dataset and imposes minimal assumption (on eigengap). With operator
perturbation theory, we are able to decompose the subsample estimation error
into a first order term and a higher order term. The proposed FunPrinSS
probability is then justified from two aspects:

  1. (i)

From a heuristic aspect, we show the proposed probability enjoys the property
of precisely controlling the first order error. In fact, the FunPrinSS
probability can be interpreted as an importance sampling for the first order
error. This means FunPrinSS for FPCA plays a role similar to that of the IMPO
sampling (He and Yan, 2020) in minimizing the difference
$\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}$ for covariance
operator.

  2. (ii)

From a rigorous aspect, we develop concentration bounds for the proposed
algorithm. The constants of the derived bounds reflect the low intrinsic
dimension nature of functional data. That is, the subsample uncertainty is
characterized by the dimension of the functional principal subspace and the
intrinsic dimension of the residual covariance operator. These dimension
parameters are usually small for functional data, and as a result, our theory
indicates the subsample size need not be very large to control the subsample
error within a desired accuracy level.

The randomized FPCA algorithm and its analysis naturally extend to the
functional linear regression setting in Section 4. For FLR as an inverse
problem in the Hilbert space, the involved operator inverse is not well
defined due to its infinitely many eigenvalues decreasing to zero. We study an
estimator whose operator inverse is truncated after its leading $R$ functional
principal components. It follows that the accuracy of the FLR estimator relies
on the precision of the subsampled estimates of functional principal
components. Based on the FunPrinSS probability, a randomized algorithm for FLR
is proposed.

For the theoretical justification the randomized FLR algorithm, we conduct a
careful analysis of the subsample prediction error and decompose it into three
components: the functional principal subspace estimation error, the truncated
inverse operator estimation error, and an orthogonal residual error. The error
decomposition reveals distinct sources of uncertainty during subsampling. Like
FPCA, our FLR theoretical bound is further established via the operator
perturbation expansion of the related error components.

The proposed randomized algorithms are closely related to the survey based
functional data analysis methods (see e.g. Cardot and Josserand, 2011; Cardot,
Degras and Josserand, 2013; Degras, 2014; Lardin-Puech, Cardot and Goga,
2014). These works study the estimation of the mean (or total) function,
covariance function, or eigenfunctions for a population (with possibly unknown
size). Classical survey sampling strategies (e.g. simple sampling without
replacement and stratified sampling) are employed to draw functional samples.
Most of the existing works focus on estimating the mean or total function.
When auxiliary variables are available, Cardot, Goga and Lardin (2013)
construct a model-assisted estimator of the mean function. For eigenfunction
estimation, Cardot et al. (2010) establish the asymptotic design unbiasedness
and asymptotic variance of the estimator for general survey sampling
probabilities. However, no specific strategy is proposed to reduce the
asymptotic variance of the eigenfunction estimator. To our best knowledge,
functional linear regression has not been considered and analyzed yet in the
survey based functional data analysis framework.

We conclude this section by introducing notations used throughout this work.
Let $\mathscr{J}:\mathscr{H}_{X}\mapsto\mathscr{H}_{X}$ be a bounded linear
operator defined on $\mathscr{H}_{X}$. The operator norm is defined by
$\|\mathscr{J}\|:=\sup_{u\in\mathscr{H}_{X}:\,\|u\|\leq 1}\|\mathscr{J}u\|$.
Assume that $\\{e_{j}\\}_{j=1}^{\infty}$ is a complete orthonormal system for
$\mathscr{H}_{X}$, the Hilbert–Schmidt norm for $\mathscr{J}$ is defined by
$\|\mathscr{J}\|_{HS}^{2}:=\sum_{j=1}^{\infty}\langle\mathscr{J}e_{j},e_{j}\rangle^{2}$.
Besides, its intrinsic dimension is computed via
$\text{intdim}(\mathscr{J}):=\mathrm{tr}(\mathscr{J})/\|\mathscr{J}\|$, where
the trace of $\mathscr{J}$ is defined as
$\mathrm{tr}(\mathscr{J}):=\sum_{j=1}^{\infty}\langle\mathscr{J}e_{j},e_{j}\rangle$.
In this work, the $N$-dimensional Euclidean space $\mathbb{R}^{N}$ is equipped
with a normalized inner product
$\langle\mathbf{a},\mathbf{b}\rangle_{N}=(1/N)\sum_{i=1}^{N}a_{i}b_{i}$ for
any vectors $\mathbf{a},\mathbf{b}\in\mathbb{R}^{N}$ and its induced norm is
defined as $\|\mathbf{a}\|_{N}=\big{(}\sum_{i=1}^{N}a_{i}^{2}/N\big{)}^{1/2}$.

The rest of the article is organized as follows. In Section 2, we develop our
randomized algorithm for FPCA. We propose the functional principal subspace
probability and a fast two-step algorithm to estimate it. The theoretical
analysis of the randomized FPCA algorithms is conducted in Section 3. In
Section 4, the algorithm and theoretical results are extended to functional
linear regression with scalar response. Simulation experiments are conducted
in Section 5 to compare distinct sampling proposals. The algorithm performance
is further examined over an astronomical real dataset in Section 6.

##  2 Functional Principal Subspace Sampling for FPCA

Functional principal component analysis (FPCA, Horváth and Kokoszka, 2012;
Ramsay, 2004) is the cornerstone for many functional data analysis methods. To
reduce the computation cost of FPCA over a large dataset, we develop a
randomized algorithm for FPCA in Section 2.1. To retain the most relevant
principal subspace information after subsampling, we propose the functional
principal subspace sampling (FunPrinSS) probability in Section 2.2. A
practical implementation of our algorithm is discussed in Section 2.3.

###  2.1 Randomized FPCA

FPCA generalizes the classical multivariate principal component analysis to
identify the modes of variations among a dataset. Classical multivariate
principal components can be computed from the eigen-decomposition of a
covariance matrix. Similarly, functional principal components can be estimated
via the eigenvalue-eigenfunction decomposition of the empirical covariance
operator
$\widehat{\mathscr{C}}_{XX}=\sum_{r=1}^{\infty}\hat{\sigma}_{r}^{2}\hat{\theta}_{r}\otimes\hat{\theta}_{r}$,
where $\hat{\sigma}_{1}^{2}\geq\hat{\sigma}_{2}^{2}\geq\cdots$ are non-
negative eigenvalues, and $\hat{\theta}_{r}$’s are the corresponding
eigenfunctions (or called functional principal components). According to the
Karhunen-Loève (KL) expansion, each sample $x_{n}$ can be represented by

$x_{n}=\sum_{r=1}^{\infty}\hat{\sigma}_{r}\hat{\xi}_{nr}\hat{\theta}_{r},$
(2.1)

where $\hat{\xi}_{nr}$ is the $r$-th principal component score of the $n$-th
sample. It is computed by $\hat{\xi}_{nr}=\langle
x_{n},\hat{\theta}_{r}\rangle/\hat{\sigma}_{r}$ when $\hat{\sigma}_{r}^{2}>0$,
and we set $\hat{\xi}_{nr}=1$ when $\hat{\sigma}_{r}^{2}=0$. Let
$\boldsymbol{\xi}_{r}=(\hat{\xi}_{1r},\cdots,\hat{\xi}_{Nr})^{T}$ denote the
vector of the $r$-th score for all samples. It satisfies the normalizing
property $\|\boldsymbol{\xi}_{r}\|_{N}=1$. From the KL expansion and the
orthonormality of the eigenfunctions, we can also find that
$\|x_{n}\|^{2}=\sum_{r=1}^{\infty}\sigma_{r}^{2}\hat{\xi}_{nr}^{2}$.

Input: Dataset $\\{x_{n}\\}_{n=1}^{N}$; sampling probability
$\\{p_{n}\\}_{n=1}^{N}$; subsample size $C$.

Output: $\tilde{\sigma}_{r}$ and $\tilde{\theta}_{r}$ for $r=1,\cdots,R$

1: for $c=1,\cdots,C$ do

2: Sample $\tilde{x}_{c}$ from $\\{x_{n}\\}_{n=1}^{N}$ and get $\tilde{p}_{c}$
with probability
$\mathbb{P}\big{(}(\tilde{x}_{c},\tilde{p}_{c})=(x_{n},p_{n})\big{)}=p_{n}.$

3: end for

4: Compute
$\widetilde{\mathscr{C}}_{XX}=\frac{1}{C}\sum_{c=1}^{C}\frac{1}{N\tilde{p}_{c}}\tilde{x}_{c}\otimes\tilde{x}_{c}$.

5: Compute the eigenvalue-eigenvector decomposition
$\widetilde{\mathscr{C}}_{XX}=\sum_{r=1}^{\infty}\tilde{\sigma}_{r}^{2}\tilde{\theta}_{r}\otimes\tilde{\theta}_{r}$.

Algorithm 2 Randomized FPCA

When the eigenvalues $\hat{\sigma}_{r}^{2}$’s have a fast decay toward zero,
the covariance operator can be well approximated by the leading eigenpairs,
$\widehat{\mathscr{C}}_{XX}\approx\sum_{r=1}^{R}\hat{\sigma}_{r}^{2}\hat{\theta}_{r}\otimes\hat{\theta}_{r}$
for some $R$. Besides, each sample can be approximated by
$x_{n}\approx\sum_{r=1}^{R}\hat{\sigma}_{r}\hat{\xi}_{nr}\hat{\theta}_{r}$. In
practice, obtaining the leading $R$ functional principal components
$\hat{\theta}_{1},\cdots,\hat{\theta}_{R}$ relies on the computation and
decomposition of the empirical $\widehat{\mathscr{C}}_{XX}$. As noted in
Section 1, computing $\widehat{\mathscr{C}}_{XX}$ is an expensive procedure
for a large dataset. We propose to estimate the full sample
$\hat{\theta}_{1},\cdots,\hat{\theta}_{R}$ based on a subsampled dataset. In
particular, we can apply Algorithm 1 to get a subsampled covariance operator
$\widetilde{\mathscr{C}}_{XX}$, and perform eigenvalue-eigenfunction
decomposition of $\widetilde{\mathscr{C}}_{XX}$ instead of
$\widehat{\mathscr{C}}_{XX}$. This is summarized in Algorithm 2 for some
sampling probability $\\{p_{n}\\}_{n=1}^{N}$, where the probability satisfies
$p_{n}\geq 0$ and $\sum_{n=1}^{N}p_{n}=1$. The algorithm will output
subsampled eigenfunction estimates
$\tilde{\theta}_{1},\cdots,\tilde{\theta}_{R}$. The effectiveness of Algorithm
2 relies on suitably choosing the sampling probability, which will be
discussed in details in Section 2.2.

###  2.2 Functional Principal Subspace Sampling Probability

Algorithm 2 leaves the sampling probability $\\{p_{n}\\}_{n=1}^{N}$
unspecified. As noted in Section 1.2, though $p_{n}\propto\|x_{n}\|^{2}$ is
optimal in controlling the error
$\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}$ in terms of Hilbert-
Schmidt norm, it is far from being optimal for the principal subspace
estimation. The problem arises because the eigenvalues
$\hat{\sigma}_{r}^{2}$’s place unequal weights on distinct principal component
directions for
$p_{n}\propto\|x_{n}\|^{2}=\sum_{r=1}^{\infty}\sigma_{r}^{2}\hat{\xi}_{nr}^{2}$.
To remove the influence of eigenvalues in determing the sampling probability,
we propose the functional principal subspace sampling (FunPrinSS) probability,
which is

$p_{n}^{\mathrm{Exact}}=\frac{\sum_{r=1}^{R}\hat{\xi}_{nr}^{2}+\|(I-\widehat{\mathscr{P}}_{R})x_{n}\|^{2}/\hat{\sigma}_{R}^{2}}{\sum_{m=1}^{N}\big{[}\sum_{r=1}^{R}\hat{\xi}_{mr}^{2}+\|(I-\widehat{\mathscr{P}}_{R})x_{m}\|^{2}/\hat{\sigma}_{R}^{2}\big{]}},$
(2.2)

where $I$ is the identity operator (i.e., $Ix=x$ for any
$x\in\mathscr{H}_{X}$), and
$\widehat{\mathscr{P}}_{R}=\sum_{r=1}^{R}\hat{\theta}_{r}\otimes\hat{\theta}_{r}$
is the projection operator of the subspace spanned by
$\hat{\theta}_{1},\cdots,\hat{\theta}_{R}$. For the $n$-th sample, the
numerator of (2.2) includes two parts: (i) the sum of squared scores in the
principal subspace $\sum_{r=1}^{R}\hat{\xi}_{nr}^{2}$ ; and (ii) the squared
norm of $(I-\widehat{\mathscr{P}}_{R})x_{n}$, which is the observation
projected out of the principal subspace. The second part of the numerator is
necessary because, in addition to the first $R$ scores, all the remaining
scores will affect the operator perturbation
$\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}$ (see Lemma 6 of
Section 3.2 and the followed discussion). However, it is not practical to
compute all the remaining scores in the infinite dimensional function space
for a large dataset. For this reason, we adopt the squared norm of the
projected residual
$\|(I-\widehat{\mathscr{P}}_{R})x_{n}\|^{2}/\hat{\sigma}_{R}^{2}$ in (2.2) to
account for the effect of the remaining scores.

The proposed sampling probability (2.2) is denoted as “Exact” because its
value is determined by the full sample eigenfunctions
$\hat{\theta}_{1},\cdots,\hat{\theta}_{R}$ and the related scores
$\hat{\xi}_{nr}$. In practice, its exact value is unknown unless we carry out
full sample FPCA computation. Fortunately, practitioners only need an
approximate value of (2.2) as input for Algorithm 2. A probability
$\\{p_{n}\\}_{n=1}^{N}$ is regarded as an approximation to (2.2) if

$p_{n}\geq p_{n}^{\mathrm{Exact}}/\beta$ (2.3)

for some fixed $\beta\geq 1$. We call a sampling probability
$\\{p_{n}\\}_{n=1}^{N}$ satisfying (2.3) nearly exact FunPrinSS probability.
Section 2.3 will develop a fast algorithm to obtain an estimate
$\\{\hat{p}_{n}\\}_{n=1}^{N}$ fulfilling the lower bound requirement (2.3)
with high probability. Our theoretical analysis of the randomized FPCA will
also be based on any sampling probability satisfying (2.3).

The normalizing constant for the sampling probability (2.2) has an interesting
interpretation. Define
$\mathscr{R}:=\sum_{r=R+1}^{\infty}\hat{\sigma}_{r}^{2}\big{(}\hat{\theta}_{r}\otimes\hat{\theta}_{r}\big{)}$
as the residual operator obtained by removing the first $R$ eigenpairs from
$\widehat{\mathscr{C}}_{XX}$, i.e.
$\mathscr{R}=\widehat{\mathscr{C}}_{XX}-\sum_{r=1}^{R}\hat{\sigma}_{r}^{2}\big{(}\hat{\theta}_{r}\otimes\hat{\theta}_{r}\big{)}$.
Let $\Delta_{R}$ be the intrinsic dimension of $\mathscr{R}$ such that

$\Delta_{R}:=\mathrm{intdim}(\mathscr{R})=\mathrm{tr}(\mathscr{\mathscr{R}})/\|\mathscr{R}\|=\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}/\hat{\sigma}_{R+1}^{2},$
(2.4)

if $\hat{\sigma}_{R+1}^{2}>0$; and $\Delta_{R}=0$ if
$\hat{\sigma}_{R+1}^{2}=0$. Then, for the normalizing constant in the
denominator of (2.2), we have

$\displaystyle\sum_{m=1}^{N}\big{[}\sum_{r=1}^{R}\hat{\xi}_{mr}^{2}+\|(I-\widehat{\mathscr{P}}_{R})$
$\displaystyle
x_{m}\|^{2}/\hat{\sigma}_{R}^{2}\big{]}=\sum_{m=1}^{N}\big{[}\sum_{r=1}^{R}\hat{\xi}_{mr}^{2}+\sum_{s=R+1}^{N}\hat{\sigma}_{s}^{2}\hat{\xi}_{ms}^{2}/\hat{\sigma}_{R}^{2}\big{]}$
$\displaystyle=N\big{(}R+\sum_{s=R+1}^{N}\hat{\sigma}_{s}^{2}/\hat{\sigma}_{R}^{2}\big{)}\leq
N\big{(}R+\Delta_{R}\big{)},$ (2.5)

where the second equality uses the score normalizing property
$\|\boldsymbol{\xi}_{r}\|_{N}=1$ for
$\boldsymbol{\xi}_{r}=(\hat{\xi}_{1r},\cdots,\hat{\xi}_{Nr})^{T}$. The last
inequality uses that $\hat{\sigma}_{R}^{2}\geq\hat{\sigma}_{R+1}^{2}$. The
above states that the normalizing constant in (2.2) is upper bounded by
$N(R+\Delta_{R})$. In particular, the quantity $R+\Delta_{R}$ is a summation
of the principal subspace dimension $R$ and the residual intrinsic dimension
$\Delta_{R}$.

###  2.3 Implementation of FunPrinSS

Input: The dataset $\\{x_{n}\\}_{n=1}^{N}$; the subsample size $C$; the number
$R$.

Output: The estimated $\\{p_{n}\\}_{n=1}^{N}$.

1: Compute
$p_{n}^{\prime}=\alpha\cdot(1/N)+(1-\alpha)\cdot\|x_{n}\|^{2}/\big{(}\sum_{m=1}^{N}\|x_{n}\|^{2}\big{)}$
for each $n$.

2: Apply Algorithm 2 with $\\{p^{\prime}_{n}\\}_{n=1}^{N}$ and $C$ to get the
pilot estimates $(\tilde{\sigma}^{\prime}_{r})^{2}$ and
$\tilde{\theta}^{\prime}_{r}$ for $r\leq R$.

3: for $n=1,\cdots,N$ do

4: Compute
$\hat{\xi}^{\prime}_{nr}=\langle\tilde{\theta}^{\prime}_{r},x_{n}\rangle/\tilde{\sigma}^{\prime}_{r}$
for $r=1,\cdots,R$.

5: Compute
$\hat{p}_{n}=\sum_{r=1}^{R}(\hat{\xi}^{\prime}_{nr})^{2}+\|(I-\widetilde{\mathscr{P}}^{\prime}_{R})x_{n}\|^{2}/(\tilde{\sigma}^{\prime}_{R})^{2}$.

6: end for

7: Normalize $\\{\hat{p}_{n}\\}_{n=1}^{N}$ to become a proper probability.

Algorithm 3 Estimating the Functional Principal Subspace Probability

Recall the exact value of (2.2) is not available unless we obtain
$\hat{\theta}_{1},\cdots,\hat{\theta}_{R}$ with full sample computation. This
would contradict the goal of reducing computational cost via subsampling. In
practice, a two-step procedure can be carried out:

>   1. Step 1.
>
> Obtain an estimate $\\{\hat{p}_{n}\\}_{n=1}^{N}$ of the exact sampling
> probability (2.2) via the randomized procedure in Algorithm 3.
>
>   2. Step 2.
>
> Plug the estimated probability $\\{\hat{p}_{n}\\}_{n=1}^{N}$ into Algorithm
> 2, which will gives the final estimator of the functional principal
> components.
>
>

Specifically, in the first step, Algorithm 3 estimates the sampling
probability (2.2) as follows. At first (Line 1–2), it executes Algorithm 2
with the sampling probability

$p_{n}^{\prime}=\alpha(1/N)+(1-\alpha)\|x_{n}\|^{2}/\Big{(}\sum_{m=1}^{N}\|x_{n}\|^{2}\Big{)},$

for some $\alpha\in[0,1]$. The probability $\\{p_{n}^{\prime}\\}_{n=1}^{N}$ is
a mixture of the uniform sampling probability and the importance sampling
probability (He and Yan, 2020). This provides a pilot estimate of the first
$R$ eigenvalues $(\tilde{\sigma}^{\prime}_{r})^{2}$ and eigenfunctions
$\tilde{\theta}^{\prime}_{r}$ in Line 2 of Algorithm 3. At the same time, a
pilot estimation of the projection operator
$\widehat{\mathscr{P}}^{\prime}_{R}=\sum_{r=1}^{R}\hat{\theta}^{\prime}_{r}\otimes\hat{\theta}^{\prime}_{r}$
is obtained. With these pilot estimates, we compute the score estimates
$\hat{\xi}^{\prime}_{nr}=\langle\tilde{\theta}^{\prime}_{r},x_{n}\rangle/\tilde{\sigma}^{\prime}_{r}$
for each sample in Line 4 of Algorithm 3. Then, the probability estimator
$\hat{p}_{n}\propto\big{(}\sum_{r=1}^{R}(\hat{\xi}_{nr}^{\prime})^{2}+\|(I-\widehat{\mathscr{P}}^{\prime}_{R})x_{n}\|^{2}/(\hat{\sigma}^{\prime}_{R})^{2}\big{)}$
is computed for each sample in Line 5. Algorithm 3 will output the normalized
$\\{\hat{p}_{n}\\}_{n=1}^{N}$. In Section 3.3, it will be shown that the
probability $\\{\hat{p}_{n}\\}_{n=1}^{N}$ computed by Algorithm 3 satisfies
the lower bound requirement (2.3) with high probability.

######  Remark 1.

Two-step procedures can often be found in the literature of randomized
algorithms. In Drineas et al. (2012), the leverage scores are approximately
computed by sketching the design matrix. In the work of Wang, Zhu and Ma
(2018); Zhang, Ning and Ruppert (2021) for generalized linear regression, the
sampling probabilities are also computed from pilot regression estimators.

######  Remark 2.

The computational complexity of the above two-step procedure can be analyzed
when each functional observation is digitally recorded as a high dimensional
vector of length $L$. In this case, the eigenfunction can be obtained by
directly applying SVD to the data matrix. The full sample computational
complexity is $\mathscr{O}(NL\min\\{N,L\\})$. On the other hand, the proposed
Algorithm 2 costs a complexity of order $\mathscr{O}(NRL+CL\min\\{C,L\\})$,
where computing the sampling probabilities via Algorithm 3 costs
$\mathscr{O}(NRL)$.

Compared with the full sample computation, the proposed algorithm has much
smaller complexity when $C\ll N$ and $R\ll L$. Specifically, the relation
$R\ll L$ is a characteristic of functional data, whose intrinsic dimension is
much smaller than the ambient dimension. Recall from Section 1, in practical
applications such as astronomical spectral processing, $L$ is in the order of
thousands, and $N$ has the magnitude of millions. Some work (Connolly et al.,
1994) has found $R=3$ principal components are enough to capture most
variability of the dataset. Theorem 8 suggests the subsample size $C$ need not
to be very large for accurate subsample estimation, due to the small dimension
parameter $R+\Delta_{R}$ for functional data.

##  3 Theoretical Analysis of Randomized FPCA with FunPrinSS

This section develops a theoretical justification of Algorithm 2 when the
employed sampling probability satisfies (2.3). The derivation is based on the
operator perturbation theory reviewed in Section 3.1, and the main results
follow in Section 3.2. After that, the theoretical justification of Algorithm
3 is presented in Section 3.3. We will show its estimated probability
satisfies the lower bound requirement (2.3) with high probability. Finally, we
discuss practical choice of the parameter $R$ in Section 3.4.

Our analysis formulates estimating the leading $R$ eigenfunctions as
estimating the projection operator
$\widehat{\mathscr{P}}_{R}=\sum_{r=1}^{R}\hat{\theta}_{r}\otimes\hat{\theta}_{r}$.
Algorithm 2 delivers a randomized estimate
$\widetilde{\mathscr{P}}_{R}=\sum_{r=1}^{R}\tilde{\theta}_{r}\otimes\tilde{\theta}_{r}$,
which is computed from the sketched data. We will focus on bounding their
difference $\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}$. The
theoretical analysis is carried out from an algorithmic perspective. In other
words, the concentration bounds are established conditional on the full
dataset, and the derivation exploits the independence of the subsampling draws
within our algorithms. The full dataset can be almost arbitrary and the
theoretical bounds provide worst-case guarantee. Neither specific distribution
nor independence assumption is imposed on the full dataset. Only a single
condition on the eigenvalues is required.

######  Condition 3.

The empirical eigenvalues of $\widehat{\mathscr{C}}_{XX}$ satisfy the
following relation:

$\hat{\sigma}_{1}^{2}>\hat{\sigma}_{2}^{2}>\cdots>\hat{\sigma}_{R}^{2}>\hat{\sigma}_{R+1}^{2}\geq\hat{\sigma}_{R+2}^{2}\geq\cdots.$

######  Remark 4.

A strictly positive eigengap
$g_{R}=\hat{\sigma}_{R}^{2}-\hat{\sigma}_{R+1}^{2}>0$ is the minimal
assumption for the identifiability of $\widehat{\mathscr{P}}_{R}$. To simplify
presentation, we further assume the leading $R$ eigenvalues
$\hat{\sigma}_{1}^{2},\cdots,\hat{\sigma}_{R}^{2}$ are distinct with
multiplicity one. If there exists repeated eigenvalues, the same theoretical
results can be obtained using a similar argument but with complications in
notation. Lastly, note we always have $\hat{\sigma}_{r}^{2}=0$ for $r>N$,
because the empirical $\widehat{\mathscr{C}}_{XX}$ is computed from $N$
samples, i.e. $\text{rank}\big{(}\widehat{\mathscr{C}}_{XX}\big{)}\leq N$.

###  3.1 Operator Perturbation Theory

We set up the tools from the operator perturbation theory (Hsing and Eubank,
2015) for further analysis of the randomized FPCA. For the operator
$\widehat{\mathscr{C}}_{XX}$, its resolvent is defined as
$\mathscr{R}_{\widehat{\mathscr{C}}_{XX}}(\eta):=(\widehat{\mathscr{C}}_{XX}-\eta
I)^{-1}$, where $I$ is the identity mapping. Let
$\Gamma_{R}:=\\{\eta\in\mathbb{C}:\
\mathrm{dist}(\eta,\,[\hat{\sigma}_{R}^{2},\hat{\sigma}_{1}^{2}])=g_{R}/2\\}$
represent the boundary of a disk in a complex plane. Every point on
$\Gamma_{R}$ has equal distance
$g_{R}/2=(\hat{\sigma}_{R}-\hat{\sigma}_{R+1})/2$ to the interval
$[\hat{\sigma}_{R}^{2},\hat{\sigma}_{1}^{2}]\subset\mathbb{R}$ on the positive
real axis. Then, the subspace projection operator can be expressed as a
contour integration $\widehat{\mathscr{P}}_{R}=-\frac{1}{2\pi
i}\oint_{\Gamma_{R}}\mathscr{R}_{\widehat{\mathscr{C}}_{XX}}(\eta)\,\mathrm{d}\eta$.

Denote $\mathscr{E}=\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}$
as the covariance operator approximation error. Then, using the result of
Koltchinskii et al. (2016), the difference of the projection operators can be
decomposed as
$\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}=L_{R}(\mathscr{E})+S_{R}(\mathscr{E})$,
where

$\displaystyle L_{R}(\mathscr{E}):$ $\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma_{R}}\mathscr{R}_{\widehat{\mathscr{C}}_{XX}}(\eta)\mathscr{E}\mathscr{R}_{\widehat{\mathscr{C}}_{XX}}(\eta)\mathrm{d}\eta,\,$
(3.1) $\displaystyle S_{R}(\mathscr{E}):$
$\displaystyle=\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}-L_{R}(\mathscr{E}).$
(3.2)

Notice the first term $L_{R}(\mathscr{E})$ is linear in $\mathscr{E}$, while
$S_{R}(\mathscr{E})$ is the residual error. Because
$\mathbb{E}(\mathscr{E})=0$ where the expectation is taken with respect to the
subsampling procedure, the error $L_{R}(\mathscr{E})$ in (3.1) also has zero
expectation $\mathbb{E}L_{R}(\mathscr{E})=0$.

Recall $g_{R}=\hat{\sigma}_{R}^{2}-\hat{\sigma}_{R+1}^{2}(>0)$ is the
eigengap. In the case
$\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|<g_{R}/3$, it
holds that

$\max_{r}|\tilde{\sigma}_{r}^{2}-\hat{\sigma}_{r}^{2}|\leq\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|<g_{R}/3.$

This implies that the first $R$ eigenvalues of $\widetilde{\mathscr{C}}_{XX}$
are strictly in the interior of the disk $\Gamma_{R}$. Therefore, the relation
$\widetilde{\mathscr{P}}_{R}=-\frac{1}{2\pi
i}\oint_{\Gamma_{R}}\mathscr{R}_{\widetilde{\mathscr{C}}_{XX}}(\eta)\,\mathrm{d}\eta$
holds for the subsampled estimator with
$\mathscr{R}_{\widetilde{\mathscr{C}}_{XX}}(\eta)=(\widetilde{\mathscr{C}}_{XX}-\eta
I)^{-1}$. In this case, the residual approximation error $S_{R}(\mathscr{E})$
has the integral expression

$S_{R}(\mathscr{E})=-\frac{1}{2\pi i}\oint_{\Gamma_{R}}\sum_{k\geq
2}[-\mathscr{R}_{\widehat{\mathscr{C}}_{XX}}(\eta)\mathscr{E}]^{k}\mathscr{R}_{\widehat{\mathscr{C}}_{XX}}(\eta)\mathrm{d}\eta.$

From the above, we can see that $L_{R}(\mathscr{E})$ can be regarded as the
first order approximation error for
$\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}$, while
$S_{R}(\mathscr{E})$ as the higher order error. In fact, according to Lemma 2
of Koltchinskii et al. (2016), the two terms $L_{R}(\mathscr{E})$ and
$S_{R}(\mathscr{E})$ can be bounded as in the following lemma.

######  Lemma 5.

For the full sample covariance operator $\widehat{\mathscr{C}}_{XX}$ and its
subsampled counterpart $\widetilde{\mathscr{C}}_{XX}$, denote
$\mathscr{E}=\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}$ as their
difference. Consider the decomposition
$\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}=L_{R}(\mathscr{E})+S_{R}(\mathscr{E})$
with $L_{R}(\mathscr{E})$ and $S_{R}(\mathscr{E})$ defined in (3.1) and (3.2),
respectively. It holds that

$(i)\quad\|L_{R}(\mathscr{E})\|\leq\big{[}1+(\hat{\sigma}_{1}^{2}-\hat{\sigma}_{R}^{2})/(\pi
g_{R})\big{]}\times\|\mathscr{E}\|;\quad(ii)\quad\|S_{R}(\mathscr{E})\|\leq
K_{R}\left(\|\mathscr{E}\|/g_{R}\right)^{2}.$

where $K_{R}:=15[1+2(\hat{\sigma}_{1}^{2}-\hat{\sigma}_{R}^{2})/(\pi g_{R})]$.

###  3.2 Main Results

We are now ready to derive our main theorem for the randomized FPCA in
Algorithm 2 with the nearly exact sampling probability (2.3). Before
presenting the formal result, we provide a heuristic justification of the
proposed functional principal subspace sampling probability (2.2) and (2.3).
The decomposition in (3.1) and (3.2) suggests that, to minimize the loss
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|$, we need to have
precise control over the first order error $L_{R}(\mathscr{E})$; at the same
time, the magnitude of $S_{R}(\mathscr{E})$ should also be contained with high
probability. Our proposed FunPrinSS probability backs up this intuition. In
fact, the probability (2.2) or (2.3) can be interpreted as an importance
sampling probability to control the first order error $L_{R}(\mathscr{E})$.
The argument is based on the following fact.

######  Lemma 6.

The first order error term $L_{R}(\mathscr{E})$ can be expressed as
$L_{R}(\mathscr{E})=(1/C)\sum_{c=1}^{C}\mathscr{Z}_{c}/(N\tilde{p}_{c})$. Each
$\mathscr{Z}_{c}$ is related to one subsampled $\tilde{x}_{c}$ and that

$\mathscr{Z}_{c}=\sum_{r=1}^{R}\sum_{s=R+1}^{\infty}\frac{\hat{\sigma}_{r}\hat{\sigma}_{s}}{\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}}\times(\tilde{\xi}_{cr}\tilde{\xi}_{cs})\times\big{[}\hat{\theta}_{r}\otimes\hat{\theta}_{s}+\hat{\theta}_{s}\otimes\hat{\theta}_{r}\big{]}.$

In the above,
$\tilde{\xi}_{cs}=\langle\tilde{x}_{c},\hat{\theta}_{s}\rangle/\hat{\sigma}_{s}$
when $\hat{\sigma}_{s}>0$ and $\tilde{\xi}_{cs}=1$ when $\hat{\sigma}_{s}=0$.

######  Remark 7.

Note the quantity $\tilde{\xi}_{cs}$ associated with $\tilde{x}_{c}$ is
computed with the full sample eigenvalue $\hat{\sigma}_{s}$ and the full
sample eigenfunction $\hat{\theta}_{s}$. Suppose the $c$-th subsample is the
$i_{c}$-th observation in the original dataset (i.e.,
$\tilde{x}_{c}=x_{i_{c}}$), and the original $x_{i_{c}}$ has full sample
scores $\hat{\xi}_{i_{c},1},\hat{\xi}_{i_{c},2},\cdots$. Recall we have set
$\tilde{p}_{c}=p_{i_{c}}$ as the sampling probability of the original
$i_{c}$-th observation. Similarly, we denote for the $c$-th subsample
$\tilde{\xi}_{cs}=\hat{\xi}_{i_{c},s}$ ($s=1,2,\cdots$) as the corresponding
full sample score.

Lemma 6 reveals that $L_{R}(\mathscr{E})$ can be expressed as a summation of
$C$ terms, and each summand $\mathscr{Z}_{c}$ is related to one subsampled
$\tilde{x}_{c}$. For the importance sampling, the sampling probability is
proportional to the size of each summand. We can show the proposed probability
(2.2) or (2.3) serve this purpose by computing the squared norm of each
summand $\mathscr{Z}_{c}$ in Lemma 6 as

$\displaystyle\|\mathscr{Z}_{c}\|_{HS}^{2}=2\sum_{r=1}^{R}$
$\displaystyle\sum_{s=R+1}^{\infty}\frac{\hat{\sigma}_{r}^{2}\hat{\sigma}_{s}^{2}}{[\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}]^{2}}\cdot\tilde{\xi}_{cr}^{2}\tilde{\xi}_{cs}^{2}$
$\displaystyle=$ $\displaystyle\quad
2\sum_{r=1}^{R}\sum_{s=R+1}^{\infty}\underbrace{\frac{\hat{\sigma}_{r}^{4}}{[\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}]^{2}}}_{f_{rs}^{2}}\cdot\tilde{\xi}_{cr}^{2}(\hat{\sigma}_{s}^{2}\tilde{\xi}_{cs}^{2}/\hat{\sigma}_{r}^{2}).$
(3.3)

Define $f_{rs}:=(1-q_{r,s})^{-1}$ with
$q_{r,s}=\hat{\sigma}_{s}^{2}/\hat{\sigma}_{r}^{2}$ for the above equation.
For $r<s$, it is obvious that $q_{r,s}\leq 1$ and $f_{rs}\geq 1$. Consider the
special case where the leading $R$ principal components dominate the overall
signal, i.e., $\hat{\sigma}_{R}^{2}\gg\hat{\sigma}_{R+1}^{2}$. In this case,
we have the approximation $f_{rs}\approx 1$ for $r(\leq R)$ and $s(>R)$. It
follows approximately that

$\displaystyle\|\mathscr{Z}_{c}\|_{HS}^{2}$ $\displaystyle\approx
2\sum_{r=1}^{R}\sum_{s=R+1}^{\infty}\tilde{\xi}_{cr}^{2}\big{(}\hat{\sigma}_{s}^{2}\tilde{\xi}_{cs}^{2}/\hat{\sigma}_{r}^{2}\big{)}\stackrel{{\scriptstyle(i)}}{{\leq}}2\Big{(}\sum_{r=1}^{R}\tilde{\xi}_{cr}^{2}\Big{)}\Big{(}\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}\tilde{\xi}_{cs}^{2}/\hat{\sigma}_{R}^{2}\Big{)}$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\leq}}\frac{1}{2}\Big{(}\sum_{r=1}^{R}\tilde{\xi}_{cr}^{2}+\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}\tilde{\xi}_{cs}^{2}/\hat{\sigma}_{R}^{2}\Big{)}^{2}$
$\displaystyle=\frac{1}{2}\Big{(}\sum_{r=1}^{R}\tilde{\xi}_{cr}^{2}+\|(I-\widehat{\mathscr{P}}_{R})\tilde{x}_{c}\|^{2}/\hat{\sigma}_{R}^{2}\Big{)}^{2}.$
(3.4)

In the above, (i) uses that $\hat{\sigma}_{r}^{2}\geq\hat{\sigma}_{R}^{2}$ for
$r\leq R$, and (ii) uses the inequality $4ab\leq(a+b)^{2}$ for
$a,b\in\mathbb{R}$. As expected, the last line of (3.4) indicates the
approximate proportional relation
$\tilde{p}_{c}^{\mathrm{Exact}}\propto\|\mathscr{Z}_{c}\|_{HS}$ for the
FunPrinSS probability (2.2). This resembles the importance sampling
$\tilde{p}_{c}\propto\|\tilde{x}_{c}\otimes\tilde{x}_{c}\|_{HS}=\|\tilde{x}_{c}\|^{2}$
for the covariance operator estimation (1.1) in He and Yan (2020). From here,
we find the proposed FunPrinSS probability can be viewed as an importance
sampling probability to control $L_{R}(\mathscr{E})$ when
$\hat{\sigma}_{R}^{2}\gg\hat{\sigma}_{R+1}^{2}$.

The above heuristic justification is developed in the case
$\hat{\sigma}_{R}^{2}\gg\hat{\sigma}_{R+1}^{2}$. Under the general eigenvalue
setting of Condition 3, to control the first order term $L_{R}(\mathscr{E})$,
we define

$G_{R}:=(1-q_{R,R+1})^{-1}=\hat{\sigma}_{R}^{2}/(\hat{\sigma}_{R}^{2}-\hat{\sigma}_{R+1}^{2}),$

and find it as an upper bound for the ratio factor $f_{rs}$ in (3.2). That is,
for $r\leq R<s$, we have

$f_{rs}=(1-q_{r,s})^{-1}\leq(1-q_{R,R+1})^{-1}=G_{R}.$

This relation holds because $f(q)=1/(1-q)$ is an increasing function for
$q\in(0,1)$. Theorem 8 below formalizes this idea and provides theoretical
guarantees for the subsampled FPCA in Algorithm 2 with the nearly exact
sampling probability (2.3). Note all the probability bounds in this work are
derived conditional on the full dataset and with respect to the subsampling.

######  Theorem 8.

Under Condition 3, define
$G_{R}=\hat{\sigma}_{R}^{2}/(\hat{\sigma}_{R}^{2}-\hat{\sigma}_{R+1}^{2})$ and
define

$V=\max\\{G_{R}^{2}Z^{2}/(2\beta),\,Z\\}\quad\text{and}\quad
L=\max\\{G_{R}Z/\sqrt{2},\,Z+1\\},$

with $Z=\beta(R+\Delta_{R})$. Then, for sampling probability
$\\{p_{n}\\}_{n=1}^{N}$ satisfying $p_{n}\geq p_{n}^{\text{Exact}}/\beta$ with
$\beta\geq 1$, and for $\epsilon$ satisfying $\epsilon\cdot
C\geq\sqrt{CV}+L/3$, it holds that

$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|\leq\epsilon+K_{R}\cdot\hat{\sigma}_{1}^{4}\epsilon^{2}/g_{R}^{2},$
(3.5)

with probability at least
$1-12(R+\Delta_{R})\exp\Big{(}-\frac{C\epsilon^{2}/2}{V+L\epsilon/3}\Big{)}.$

In the above theorem, $G_{R}$ can be viewed as a measurement of the eigengap
$\hat{\sigma}_{R}^{2}-\hat{\sigma}_{R+1}^{2}$. The dimension $R+\Delta_{R}$
discussed in (2.5) also appears in multiple places above. On the right hand
side of the error bound (3.5), the first term $\epsilon$ bounds the first
order error $L_{R}(\mathscr{E})$ in (3.1), while the second term
$K_{R}\cdot\hat{\sigma}_{1}^{4}\epsilon^{2}/g_{R}^{2}$ bounds the high order
$S_{R}(\mathscr{E})$ in (3.2).

According to (3.5) of Theorem 8, the subsample error
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|$ can be controlled
within a precision level $\epsilon$ with high probability. For a large enough
$C$ and with probability at least $0.9$, the bound (3.5) holds for the
precision level $\epsilon$ of magnitude order
$\mathscr{O}\big{(}(R+\Delta_{R})\log^{1/2}[120(R+\Delta_{R})]/C^{1/2}\big{)}$.
The precision level $\epsilon$ is inversely proportional to $C^{1/2}$, and is
proportional to the dimension parameter $R+\Delta_{R}$ with an additional
logarithm factor. The nature of functional data usually indicates the
dimension parameter $R+\Delta_{R}$ is small, though the ambient space
dimension is infinite. This means the subsample size $C$ need not be very
large to control the subsample error toward a small accuracy level $\epsilon$.

###  3.3 Theoretical Justification of Algorithm 3

Algorithm 3 can be justified for providing the pilot probability estimate. We
show that its estimated probability $\\{\hat{p}_{n}\\}_{n=1}^{N}$ remains
close to the target (2.2) in the sense of (2.3). Its proof is derived in
Section E in the Appendix.

######  Proposition 9.

Set $\Delta_{0}=\mathrm{intdim}(\widehat{\mathscr{C}}_{XX})$,
$\gamma_{0}=\sqrt{(\Delta_{0}+1)\log(120\Delta_{0})}$ and
$G_{R}=(1-q_{R,R+1})^{-1}$. For large enough subsample size $C$, with
probability at least $0.9$, the pilot sampling probability
$\\{\hat{p}_{n}\\}_{n=1}^{N}$ computed by Algorithm 3 with $\alpha=0.5$
satisfy

$\min\Big{\\{}\frac{p^{\mathrm{Exact}}}{\hat{p}_{n}},\,\frac{\hat{p}_{n}}{p^{\mathrm{Exact}}}\Big{\\}}\,\geq\,\frac{1-\gamma_{1}C^{-1/2}-\gamma_{2}C^{-1}}{1+\gamma_{1}C^{-1/2}+\gamma_{2}C^{-1}},$

where
$\gamma_{1}=35(1+G_{R})(R+\Delta_{R})+8\hat{\sigma}_{1}^{2}\gamma_{0}/\hat{\sigma}_{R}^{2}$
and
$\gamma_{2}=32K_{R}\hat{\sigma}_{1}^{6}G_{R}^{2}/(g_{R}^{2}\hat{\sigma}_{R}^{2})$.

One implication of Proposition 9 is that, when the subsample size $C$ is large
enough, it holds with high probability that $\hat{p}_{n}\geq
p^{\mathrm{Exact}}/\beta$ for
$\beta=\frac{1+\gamma_{1}C^{-1/2}+\gamma_{2}C^{-1}}{1-\gamma_{1}C^{-1/2}-\gamma_{2}C^{-1}}$.
This corresponds to our definition of the nearly exact FunPrinSS probability
(2.3). Besides, the value of $\beta(\geq 1)$ will converge to one as $C$
increases to infinity.

###  3.4 Choice of $R$

For functional principal component analysis, the subspace dimension $R$ is
usually chosen based on the fraction of the variance explained (FVE) by the
leading $R$ eigenfuctions. That is, $R$ is chosen such that the FVE is greater
than a pre-specified threshold value (e.g. $90\%$). Given the full sample
estimate $\widehat{\mathscr{P}}_{R}$ , its fraction of variance explained can
be computed by

$\widehat{\text{FVE}}=\frac{\sum_{r=1}^{R}\hat{\sigma}_{r}^{2}}{\sum_{r=1}^{\infty}\hat{\sigma}_{r}^{2}}=\frac{\sum_{n=1}^{N}\|\widehat{\mathscr{P}}_{R}x_{n}\|^{2}}{\sum_{n=1}^{N}\|x_{n}\|^{2}};$

and correspondingly, the fraction of variance explained by the subsampled
estimate $\widetilde{\mathscr{P}}_{R}$ is
$\widetilde{\text{FVE}}=\sum_{n=1}^{N}\|\widetilde{\mathscr{P}}_{R}x_{n}\|^{2}/\sum_{n=1}^{N}\|x_{n}\|^{2}$.
As a direct consequence of Theorem 8, we can bound the difference between the
$\widehat{\text{FVE}}$ of the full sample estimator and the
$\widetilde{\text{FVE}}$ of the subsampled estimator.

######  Corollary 10.

Under the conditions of Theorem 8, it holds that

$\big{|}\widehat{\mathrm{FVE}}-\widetilde{\mathrm{FVE}}\big{|}\leq\epsilon+K_{R}\cdot\hat{\sigma}_{1}^{4}\epsilon^{2}/g_{R}^{2},$
(3.6)

with probability at least
$1-12(R+\Delta_{R})\exp\big{(}-\frac{C\epsilon^{2}/2}{V+L\epsilon/3}\big{)}.$

###### Proof.

With Theorem 8, the result (3.6) is straightforward by observing

$\displaystyle\Big{|}\frac{\sum_{n=1}^{N}\|\widehat{\mathscr{P}}_{R}x_{n}\|^{2}}{\sum_{n=1}^{N}\|x_{n}\|^{2}}-\frac{\sum_{n=1}^{N}\|\widetilde{\mathscr{P}}_{R}x_{n}\|^{2}}{\sum_{n=1}^{N}\|x_{n}\|^{2}}\Big{|}$
$\displaystyle=$ $\displaystyle\Big{|}\frac{\sum_{n=1}^{N}\langle
x_{n},(\widehat{\mathscr{P}}_{R}-\widetilde{\mathscr{P}}_{R})x_{n}\rangle}{\sum_{n=1}^{N}\|x_{n}\|^{2}}\Big{|}\leq\epsilon+K_{R}\cdot\hat{\sigma}_{1}^{4}\epsilon^{2}/g_{R}^{2}.$

The last inequality employs $|\langle
x_{n},(\widehat{\mathscr{P}}_{R}-\widetilde{\mathscr{P}}_{R})x_{n}\rangle|\leq\|\widehat{\mathscr{P}}_{R}-\widetilde{\mathscr{P}}_{R}\|\cdot\|x_{n}\|^{2}$.
∎

Corollary 10 indicates $\widetilde{\text{FVE}}$ is a reliable estimate of
$\widehat{\text{FVE}}$. When applying the randomized algorithm, we can use the
criterion $\widetilde{\text{FVE}}$ to select $R$. After obtaining the
subsampled estimate $\widetilde{\mathscr{P}}_{R}$, computing
$\widetilde{\text{FVE}}$ requires one additional scan over the whole dataset
and have a computational complexity $\mathscr{O}(NRL)$. This complexity is no
more than that of Algorithm 3, see Remark 2.

##  4 Extenstion: Randomized Functional Linear Regression

A considerable number of functional data models rely on functional principal
components for dimension reduction. One of these models is functional linear
regression (FLR) with scalar response. In this section, we extend the results
from Section 2 and Section 3 to FLR and develop its theoretical guarantees.

###  4.1 Randomized FLR

FLR is a generalization of the classical linear regression to the functional
setting. For a functional predictor $x_{n}\in\mathscr{H}_{X}$ and a scalar
response $Y_{n}\in\mathbb{R}$, FLR assumes the relation $Y_{n}=\alpha+\langle
x_{n},\Psi\rangle+\epsilon_{n}$ for $n=1,\cdots,N$. In this model,
$\alpha\in\mathbb{R}$ is the intercept term, $\Psi\in\mathscr{H}_{X}$ is the
regression function and $\epsilon_{n}$ is some random error with zero mean.
When both $x_{n}$ and $\Psi$ belong to $L_{2}(T)$ for some compact interval
$T$, the regression model is written as
$Y_{n}=\alpha+\int_{T}x_{n}(t)\Psi(t)\,\mathrm{d}t+\epsilon_{n}$. Given any
estimator $\widehat{\Psi}$, the intercept is directly available via
$\hat{\alpha}=(1/N)\sum_{n=1}^{N}(Y_{n}-\langle x_{n},\widehat{\Psi}\rangle)$.
In the following, we focus on the estimation of the regression function $\Psi$
by assuming that $Y_{n}$ and $x_{n}$ have been centered with zero mean, such
that $\hat{\alpha}=0$.

Considerable amount of works has been proposed for estimating $\Psi$ (e.g. Yao
et al., 2005; Yuan and Cai, 2010). One of the commonly used estimators for
$\Psi$ is based on the truncated inverse of the empirical covariance operator
$\widehat{\mathscr{C}}_{XX}=\frac{1}{N}\sum_{i=1}^{N}x_{i}\otimes x_{i}$; see
Hall et al. (2007); Cardot, Mas and Sarda (2007). Suppose the covariance
operator admits the eigenvalue-eigenfunction decomposition
$\widehat{\mathscr{C}}_{XX}=\sum_{r=1}^{\infty}\hat{\sigma}_{r}^{2}\cdot\hat{\theta}_{r}\otimes\hat{\theta}_{r}$.
Its rank-$R$ truncated inverse is defined as
$\widehat{\mathscr{C}}_{XX}^{+}=\sum_{r=1}^{R}(1/\hat{\sigma}_{r}^{2})\cdot\hat{\theta}_{r}\otimes\hat{\theta}_{r}$.
Based on this, the regression function is estimated by
$\widehat{\Psi}=\widehat{\mathscr{C}}_{XX}^{+}\hat{z}$ with
$\hat{z}=(1/N)\sum_{n=1}^{N}Y_{n}x_{n}$.

Input: Dataset $\\{(x_{n},Y_{n})\\}_{n=1}^{N}$; sampling probability
$\\{p_{n}\\}_{n=1}^{N}$; subsample size $C$.

Output: the regression function $\widetilde{\Psi}$

1: for $c=1,\cdots,C$ do

2: Sample $(\tilde{x}_{c},\tilde{Y}_{c})$ from $\\{(x_{n},Y_{n})\\}_{n=1}^{N}$
and get $\tilde{p}_{c}$ with
$\mathbb{P}\big{(}(\tilde{x}_{c},\tilde{Y}_{c},\tilde{p}_{c})=(x_{n},Y_{n},p_{n})\big{)}=p_{n}.$

3: end for

4: Compute
$\widetilde{\mathscr{C}}_{XX}=\frac{1}{C}\sum_{c=1}^{C}\frac{1}{N\tilde{p}_{c}}\tilde{x}_{c}\otimes\tilde{x}_{c}$
and
$\tilde{z}=\frac{1}{C}\sum_{c=1}^{C}\frac{1}{N\tilde{p}_{c}}\tilde{Y}_{c}\tilde{x}_{c}$.

5: Compute $\widetilde{\Psi}=\widetilde{\mathscr{C}}^{+}_{XX}\tilde{z}$

Algorithm 4 Base Randomized FLR

To reduce the cost of computing $\widehat{\Psi}$ from the full sample, we
adapt the randomized algorithm to estimate it based on a subset of the data.
The algorithm will sample a subset
$\\{(\tilde{x}_{c},\tilde{Y}_{c})\\}_{c=1}^{C}$ of observation pairs from the
full sample $\\{(x_{n},Y_{n})\\}_{n=1}^{N}$. The sampling is taken with
replacement and according to some probability distribution
$\\{p_{n}\\}_{n=1}^{N}$. Given the subsampled pairs, we can compute two
quantities,
$\widetilde{\mathscr{C}}_{XX}=\frac{1}{C}\sum_{c=1}^{C}\frac{1}{N\tilde{p}_{c}}\tilde{x}_{c}\otimes\tilde{x}_{c}$
and
$\tilde{z}=\frac{1}{C}\sum_{c=1}^{C}\frac{1}{N\tilde{p}_{c}}\tilde{Y}_{c}\tilde{x}_{c}$,
which are unbiased estimators of $\widehat{\mathscr{C}}$ and $\hat{z}$,
respectively. Finally, the full sample regression function is estimated by its
subsampled counterpart
$\widetilde{\Psi}=\widetilde{\mathscr{C}}^{+}_{XX}\tilde{z}$. This procedure
is listed in Algorithm 4.

In Algorithm 4, notice the truncated inverse
$\widetilde{\mathscr{C}}^{+}_{XX}$ is computed from the eigenvalue-
eigenfunction decomposition of the subsampled $\widetilde{\mathscr{C}}_{XX}$.
It is evident that the accuracy of the randomized estimator $\widetilde{\Psi}$
relies on the precision of the principal subspace estimated by
$\widetilde{\mathscr{C}}_{XX}$. Our proposed FunPrinSS probability (2.2) is a
proper choice for the randomized FLR to capture the principal subspace
information.

In practice, the FunPrinSS probability (2.2) is unknown without full sample
computation. We can carry out a two-step procedure similarly to that of
Section 2.3. Firstly, the FunPrinSS probability is estimated via Algorithm 3.
Secondly, the estimated probability is supplied to Algorithm 4 to obtain
$\widetilde{\Psi}$. As noted in Proposition 9, the probability obtained in the
first step is close to the target FunPrinSS probability with high probability
when the subsample size $C$ is large enough.

Application of the randomized FLR algorithm also requires the specification of
the principal subspace dimension $R$. Practitioners can choose $R$ such that
the fraction of variance explained (FVE) by the principal subspace is close to
$100\%$. In this way, most of the variability of the predictor $x$ is captured
and the remaining variability can be regarded as being negligible.

######  Remark 11.

When the functional observation is recorded as a high dimensional vector of
length $L$, the computational complexity of such two-step procedure can be
similarly analyzed as in Remark 2. In this case, the computational complexity
of FLR with full sample is $\mathscr{O}(NL\min\\{L,N\\}+L(R+N))$. The proposed
two-step procedure for randomized FLR requires
$\mathscr{O}(CL\min\\{C,L\\}+(C+NR)L)$. Compared with the full sample
computation, the proposed algorithm has much smaller complexity when $C\ll N$
and $R\ll L$.

###  4.2 Analysis of Randomized FLR with FunPrinSS

We continue to theoretically analyze the subsampled regression estimator
$\widetilde{\Psi}$ obtained from Algorithm 4 with the nearly exact sampling
probability (2.3). Toward this target, we express the full sample estimator
$\widehat{\Psi}$ and the subsampled estimator $\widetilde{\Psi}$ in compact
expressions. These expressions resemble the commonly used ones for the
classical linear regression. First of all, define an operator
$\mathscr{T}:\mathscr{H}_{X}\to\mathbb{R}^{N}$ such that an element
$u\in\mathscr{H}_{X}$ is mapped to

$\mathscr{T}u=(\langle x_{1},u\rangle,\cdots,\langle x_{N},u\rangle)^{T}.$
(4.1)

Let $\mathbf{Y}=(y_{1},\cdots,y_{N})^{T}$ and
$\boldsymbol{\epsilon}=(\epsilon_{1},\cdots,\epsilon_{N})^{T}$ be
$N$-dimensional vectors of the responses and the residuals, respectively.
Then, the regression model for all $N$ observations can be expressed as
$\mathbf{Y}=\mathscr{T}\Psi+\boldsymbol{\epsilon}$. The operator $\mathscr{T}$
plays a similar role as the design matrix for classical linear regression.

Recall from the end of Section 1, the Euclidean space $\mathbb{R}^{N}$ is
equipped with the normalized inner product $\langle\cdot,\cdot\rangle_{N}$. We
denote $\mathscr{T}^{*}:\mathbb{R}^{N}\to\mathscr{H}_{X}$ as the adjoint
operator of $\mathscr{T}$. For any $\mathbf{a}\in\mathbb{R}^{N}$ and
$u\in\mathscr{H}_{X}$, it holds that

$\langle\mathscr{T}^{*}\mathbf{a},u\rangle\stackrel{{\scriptstyle(i)}}{{=}}\langle\mathbf{a},\mathscr{T}u\rangle_{N}\stackrel{{\scriptstyle(ii)}}{{=}}\frac{1}{N}\sum_{n=1}^{N}a_{n}\langle
x_{n},u\rangle=\big{\langle}\sum_{n=1}^{N}a_{n}x_{n}/N,\ u\big{\rangle}.$

In the above, (i) uses the definition of adjoint operator, and (ii) employs
(4.1). From the above equation, we find that the adjoint operator
$\mathscr{T}^{*}$ maps any $\mathbf{a}\in\mathbb{R}^{N}$ to
$(1/N)\sum_{n=1}^{N}a_{n}x_{n}$. Moreover, for any $u\in\mathscr{H}_{X}$, we
can find that

$\mathscr{T}^{*}\mathscr{T}u=\mathscr{T}^{*}\begin{pmatrix}\langle
x_{1},u\rangle\\\ \vdots\\\ \langle
x_{N},u\rangle\end{pmatrix}=\frac{1}{N}\sum_{n=1}^{N}\langle x_{n},u\rangle
x_{n}=\widehat{\mathscr{C}}_{XX}u.$

The above equation implies that
$\widehat{\mathscr{C}}_{XX}=\mathscr{T}^{*}\mathscr{T}$. Similarly, we can
verify $\hat{z}=\mathscr{T}^{*}\mathbf{Y}$. Combining these results, the full
sample estimator of the regression function can be equivalently expressed as
$\widehat{\Psi}=\widehat{\mathscr{C}}_{XX}^{+}\hat{z}=(\mathscr{T}^{*}\mathscr{T})^{+}(\mathscr{T}^{*}\mathbf{Y})$.

To find a closed form expression for the subsampled $\widetilde{\Psi}$, we
represent the subsampling process by an operator
$\mathscr{D}:\mathbb{R}^{N}\to\mathbb{R}^{C}$. As a mapping between Euclidean
spaces, $\mathscr{D}$ can be represented by a matrix
$\mathbf{D}\in\mathbb{R}^{C\times N}$. The $(c,n)$-th entry of the matrix
$\mathbf{D}$ is $D_{cn}=1/\sqrt{Np_{n}}$ if $x_{n}$ is selected as the $c$-th
subsample; all the other entries are zero. The adjoint operator of
$\mathscr{D}$ is denoted as $\mathscr{D}^{*}:\mathbb{R}^{C}\to\mathbb{R}^{N}$.
Note the Euclidean space $\mathbb{R}^{C}$ is also equipped with the normalized
inner product $\langle\cdot,\cdot\rangle_{C}$, such that for
$\mathbf{a}^{\prime},\mathbf{b}^{\prime}\in\mathbb{R}^{C}$ we have
$\langle\mathbf{a}^{\prime},\mathbf{b}^{\prime}\rangle_{C}=(1/C)\sum_{c=1}^{C}a^{\prime}_{c}b^{\prime}_{c}$.
Now for any $\mathbf{a}^{\prime}\in\mathbb{R}^{C}$ and
$\mathbf{b}\in\mathbb{R}^{N}$, we have

$\langle\mathbf{a}^{\prime},\mathbf{D}\mathbf{b}\rangle_{C}=\frac{1}{C}\langle\mathbf{a}^{\prime},\mathbf{D}\mathbf{b}\rangle_{2}=\frac{1}{C}\langle\mathbf{D}^{T}\mathbf{a}^{\prime},\mathbf{b}\rangle_{2}=\langle(N/C)\mathbf{D}^{T}\mathbf{a}^{\prime},\mathbf{b}\rangle_{N},$

where $\langle\cdot,\cdot\rangle_{2}$ is the classical (unnormalized) inner
product for Euclidean space, which is the summation of element-wise products.
From the above equation, we can see that the matrix representation of
$\mathscr{D}^{*}$ is $(N/C)\mathbf{D}^{T}$, which is the transpose of
$\mathbf{D}$ and multiplied by a ratio factor $N/C$. With some computation, we
can verify that for any $u\in\mathscr{H}_{X}$

$\mathscr{T}^{*}\mathscr{D}^{*}\mathscr{D}\mathscr{T}u=(N/C)\mathscr{T}^{*}\mathbf{D}^{T}\mathbf{D}\begin{pmatrix}\langle
x_{1},u\rangle\\\ \vdots\\\ \langle
x_{N},u\rangle\end{pmatrix}=\frac{1}{C}\sum_{c=1}^{C}\frac{1}{N\tilde{p}_{c}}\langle\tilde{x}_{c},u\rangle\tilde{x}_{c}=\widetilde{\mathscr{C}}_{XX}u.$

It follows that
$\widetilde{\mathscr{C}}_{XX}=\mathscr{T}^{*}\mathscr{D}^{*}\mathscr{D}\mathscr{T}$,
and $\tilde{z}=\mathscr{T}^{*}\mathscr{D}^{*}\mathscr{D}\mathbf{Y}$ can be
verified due to a similar reasoning. This leads to a compact expression for
the randomized regression estimator
$\widetilde{\Psi}=\widetilde{\mathscr{C}}^{+}_{XX}\tilde{z}=(\mathscr{T}^{*}\mathscr{D}^{*}\mathscr{D}\mathscr{T})^{+}(\mathscr{T}^{*}\mathscr{D}^{*}\mathscr{D}\mathbf{Y})$.
The above results are summarized in Lemma 12 below, which gives an alternative
expression of the full sample estimator and the subsampled estimator obtained
via Algorithm 4.

######  Lemma 12.

The full sample regression estimator can be expressed as

$\widehat{\Psi}=(\mathscr{T}^{*}\mathscr{T})^{+}(\mathscr{T}^{*}\mathbf{Y}),$

and the subsampled estimator can be expressed as

$\widetilde{\Psi}=(\mathscr{T}^{*}\mathscr{D}^{*}\mathscr{D}\mathscr{T})^{+}(\mathscr{T}^{*}\mathscr{D}^{*}\mathscr{D}\mathbf{Y}).$

Now we show that the subsampled estimator provides a good approximation of the
full sample estimator under suitable criterion. Our aim is to characterize the
closeness between the full sample estimator $\widehat{\Psi}$ and the
subsampled estimator $\widetilde{\Psi}$. In the literature of functional
linear regression, strong consistency of the estimated regression function is
not generally possible. The estimators are only guaranteed to converge in some
weak norm such as the prediction error on a new sample (see Cardot, Mas and
Sarda, 2007; Yuan and Cai, 2010). We choose to measure the difference
$\widetilde{\Psi}-\widehat{\Psi}$ based on the full sample prediction error
$\|\mathscr{T}\widetilde{\Psi}-\mathscr{T}\widehat{\Psi}\|_{N}^{2}=\frac{1}{N}\sum_{n=1}^{N}\langle
x_{n},\widetilde{\Psi}-\widehat{\Psi}\rangle^{2}$. When only a fraction of the
whole dataset is employed to estimate $\widetilde{\Psi}$, this metric
partially gauges the quality of $\widetilde{\Psi}$ by the samples out of the
subsampled (training) dataset. Our theoretical target boils down to the
analysis of the error $\mathscr{T}\widetilde{\Psi}-\mathscr{T}\widehat{\Psi}$.
Lemma 13 below provides a decomposition of the subsampling error, the proof of
which can be found in Section D in the Appendix. Recall that $x_{n}$ can be
represented as
$x_{n}=\sum_{r=1}^{\infty}\hat{\sigma}_{r}\hat{\xi}_{nr}\hat{\theta}_{r}$ due
to the Karhunen-Loève expansion (2.1). The following decomposition involves an
operator $\mathscr{T}^{+}:\mathbb{R}^{N}\to\mathscr{H}_{X}$, which maps
$\mathbf{a}\in\mathbb{R}^{N}$ to
$\mathscr{T}^{+}\mathbf{a}=\frac{1}{N}\sum_{n=1}^{N}a_{n}x_{n}^{+}\in\mathscr{H}_{X}$,
with $x_{n}^{+}=\sum_{r=1}^{R}\hat{\xi}_{nr}\hat{\theta}_{r}/\hat{\sigma}_{r}$
corresponding to each $x_{n}$.

######  Lemma 13.

The subsampling error in full sample prediction has the following
decomposition

$\displaystyle\mathscr{T}\widetilde{\Psi}-\mathscr{T}\widehat{\Psi}=\underbrace{\mathscr{T}(\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R})\mathscr{T}^{+}\mathbf{Y}}_{\mathrm{I}}$
$\displaystyle+\underbrace{\mathscr{T}(\widetilde{\mathscr{C}}^{+}_{XX}-\widehat{\mathscr{C}}^{+}_{XX})\mathscr{T}^{*}(\mathscr{D}^{*}\mathscr{D})\mathbf{Y}^{\perp}}_{\mathrm{II}}$
$\displaystyle\qquad+\underbrace{\mathscr{T}\widehat{\mathscr{C}}^{+}_{XX}\mathscr{T}^{*}(\mathscr{D}^{*}\mathscr{D})\mathbf{Y}^{\perp}}_{\mathrm{III}}.$
(4.2)

In the above, $\mathbf{Y}^{\perp}=\mathbf{Y}-\mathscr{T}\widehat{\Psi}$ is the
residual of the full sample estimator. In addition, we have that

$\mathscr{T}\widehat{\mathscr{C}}_{XX}^{+}\mathscr{T}^{*}\mathbf{Y}^{\perp}=\mathbf{0}.$
(4.3)

Lemma 13 shows that the subsampling error
$\mathscr{T}\widetilde{\Psi}-\mathscr{T}\widehat{\Psi}$ can be decomposed into
three interpretable terms. The first term (I) accounts for the subsampling
estimation error of the principal subspace
$\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}$. The second term (II)
arises due to the difference between the truncated inverses
$\widetilde{\mathscr{C}}^{+}_{XX}-\widehat{\mathscr{C}}^{+}_{XX}$. The third
term (III) can be interpreted as follows. Notice we have
$\mathscr{T}\widehat{\mathscr{C}}_{XX}^{+}\mathscr{T}^{*}\mathbf{Y}^{\perp}=\mathbf{0}$
by the last conclusion of Lemma 13. This relation holds due to the similar
reason that the residual is orthogonal to the column space of the design
matrix for classical linear regression. In fact, (4.3) holds because our
$\mathbf{Y}^{\perp}$ is orthogonal to the space spanned by the first $R$ score
vectors , i.e.,
$\text{span}(\boldsymbol{\xi}_{1},\cdots,\boldsymbol{\xi}_{R})$ with
$\boldsymbol{\xi}_{r}=(\hat{\xi}_{1r},\cdots,\hat{\xi}_{Nr})^{T}$. This
orthogonality is violated due to subsampling, which explains the appearance of
the last term (III) in (4.2). With Lemma 13 and the operator perturbation
expansion techniques in Section 3.1, we can establish the following result for
the randomized FLR.

######  Theorem 14.

Under Condition 3, set
$Z_{1}=\hat{\sigma}_{1}^{6}K_{R}/(g_{R}^{2}\hat{\sigma}_{R}^{2})$ and
$Z_{2}=\beta(R+\Delta_{R})$. In addition, let
$V=2+(2+G_{R}^{2})Z_{2}^{2}/\beta$ and
$L=\sqrt{2R}+(\sqrt{2}+G_{R}/\sqrt{2})Z_{2}$. Suppose $C$ subsamples are
obtained according to the probability $\\{p_{n}\\}$ satisfying $p_{n}\geq
p_{n}^{\mathrm{Exact}}/\beta$ with some $\beta\geq 1$. Then with probability
at least

$\displaystyle
1-16(R+\Delta_{R})\exp\Big{(}-\frac{C\epsilon^{2}}{V+L\epsilon/3}\Big{)}-\frac{3Z_{2}}{\epsilon^{2}C},$
(4.4)

it holds that

$\displaystyle\|\mathscr{T}\widetilde{\Psi}-\mathscr{T}\widehat{\Psi}\|_{N}$
$\displaystyle\leq\underbrace{\big{(}\epsilon+\hat{\sigma}_{R}Z_{1}/\hat{\sigma}_{1}\epsilon^{2}\big{)}\|\mathbf{Y}\|_{N}}_{a}+\underbrace{4\big{(}\epsilon+Z_{1}\epsilon^{2}\big{)}\|\mathbf{Y}^{\perp}\|_{N}}_{b}.$
(4.5)

In the above, $\epsilon$ satisfy $\epsilon\geq\sqrt{V/C}+L/(3C)$ and
$\epsilon\leq\min\\{1,g_{R}/3\\}$.

Theorem 14 directly implies that, given a subsample size $C$, we can control
$\|\mathscr{T}\widetilde{\Psi}-\mathscr{T}\widehat{\Psi}\|_{N}$ within a
precision level $\epsilon$ in (4.5) with high probability. For a large enough
$C$ and with probability at least $0.9$, the bound (4.5) holds for the
precision level $\epsilon$ of magnitude order
$\mathscr{O}\big{(}(R+\Delta_{R})\log^{1/2}[320(R+\Delta_{R})]/C^{1/2}\big{)}$.

The right hand side of (4.5) contains two parts. The first term (a) bounds the
subspace difference (I) in (4.2), and it is an error term relative to
$\|\mathbf{Y}\|_{N}$. The second term (b) bounds (II) and (III) in (4.2)
together, and its magnitude is relative to $\|\mathbf{Y}^{\perp}\|_{N}$.

To help get a better understanding of the two terms on the right hand side of
(4.5), we can relate them to similar quantities in finite dimensional
problems. For the randomized classical linear regression (Drineas et al.,
2011), the error bound can be simplified to one term containing
$\|\mathbf{Y}^{\perp}\|_{N}$ . This is because their corresponding predictors
$x_{n}$ span a finite dimensional space. Their work agrees with our result
when $\mathrm{dim}(\mathscr{H}_{X})$ is finite and relatively small, and we
directly take $R$ to be $\mathrm{dim}(\mathscr{H}_{X})$ without truncation.
Indeed, with large enough subsample size $C$, the finite dimensional space can
be easily recovered, and
$\widehat{\mathscr{P}}_{R}=\widetilde{\mathscr{P}}_{R}$ is guaranteed with
high probability. This implies that the term (I) will vanish in (4.2) for
linear regression, and consequently the first term (a) will disappear in
(4.5). The result (4.5) would become an error bound whose magnitude is only
related to the size of the residual, $\|\mathbf{Y}^{\perp}\|_{N}$.

When $\mathrm{dim}(\mathscr{H}_{X})$ is finite and
$R<\mathrm{dim}(\mathscr{H}_{X})$, our randomized FLR corresponds to
randomized multivariate principal component regression (PCR). In the
literature, additive error bound relative to $\|\mathbf{Y}\|_{2}$ has been
established for the randomized PCR (see. Theorem 5 of Mor-Yosef and Avron,
2019). Our theoretical bound (4.5) reveals different sources of uncertainties
whose sizes are proportional to $\|\mathbf{Y}^{\perp}\|_{N}$ or
$\|\mathbf{Y}\|_{N}$. This is a consequence of the subsample error
decomposition in Lemma 13.

Notice the precision parameter $\epsilon$ has an upper bound $g_{R}/3$ in the
theorem. This is because we need to control the difference between the
truncated inverses
$\widetilde{\mathscr{C}}^{+}_{XX}-\widehat{\mathscr{C}}^{+}_{XX}$. The
difference is generally unbounded unless we require
$\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|$ to be smaller
than the eigengap $g_{R}$. Thereby, the theorem is derived conditional on the
event
$\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|\leq(1/3)g_{R}$,
which holds at least with the stated probability (4.4).

##  5 Simulation Study

This section conducts a simulation study for assessing the randomized
functional principal component analysis in Algorithm 2 and the randomized
functional linear regression in Algorithm 4. Specifically, for these
randomized algorithms, we will consider three sampling probabilities: 1) the
uniform sampling probability (UNIF), $p_{n}=1/N$; 2) the importance sampling
probability (IMPO, He and Yan, 2020),
$p_{n}=\|x_{n}\|^{2}/\sum_{m=1}^{N}\|x_{m}\|^{2}$; 3) the proposed functional
principal subspace sampling probability (FunPrinSS) estimated by Algorithm 3
with $\alpha=0.5$. These sampling probabilities will be supplied to the
randomized algorithms, and the accuracy of the subsampled estimation will be
compared.

###  5.1 Data Generation

In each replication of the simulation study, we will generate one synthetic
dataset $\\{x_{n}\\}_{n=1}^{N}$ as follows. Each observation $x_{n}$ is a
continuous function defined over the compact interval $T=[0,1]$. In
particular, it is generated according to
$x_{n}(t)=\sum_{r=1}^{50}\sigma_{r}\xi_{nr}\theta_{r}(t)$, where the
$\theta_{r}$ is set as the Fourier function which is
$\theta_{r}(t)=\sqrt{2}\sin(2\pi rt)$ when $r$ is odd, and
$\theta_{r}(t)=\sqrt{2}\cos(2\pi rt)$ when $r$ is even. The scores
$\xi_{nr}$’s are independent random variables with zero mean and unit
variance. For the random function $x_{n}$, its covariance operator
$\mathscr{C}_{XX}$ is determined by the integral
$(\mathscr{C}_{XX}u)(t^{\prime})=\int_{0}^{1}k(t,t^{\prime})u(t)\,\mathrm{d}t$
for any $u\in L_{2}([0,1])$, with integral kernel
$k(t,t^{\prime})=\sum_{j=1}^{p}\sigma_{r}^{2}\theta_{r}(t)\theta_{r}(t^{\prime})$.
We can see the $r$-th eigenvalue and eigenfunction of $\mathscr{C}_{XX}$ are
$\sigma_{r}^{2}$ and $\theta_{r}$, respectively.

The settings to specify the eigenvalues $\\{\sigma_{r}^{2}\\}_{r=1}^{50}$ and
the distribution of $\xi_{nr}$ is described in Table 1. The eigenvalues
$\\{\sigma^{2}_{r}\\}_{r=1}^{50}$ are specified as two cases: Exponential
Decay (ED) and Polynomial Decay (PD). The decline rate of the PD eigenvalues
is slower than the ED eigenvalues. Hence, the eigengap for the PD setting is
smaller, and the estimation of the subspace spanned by the first R
eigenfunctions is more difficult.

Similar to the simulation setting of Ma, Mahoney and Yu (2015), the random
scores, $\xi_{nr}$’s, are independently drawn from one of the three
distributions: Nearly Uniform (NU) distribution, for which $\xi_{nr}$’s are
drawn independently from standard normal distribution $N(0,1)$; Moderately
Nonuniform (MN) distribution, where the scores $\xi_{nr}$’s are independently
generated by t-distribution with three degree of freedom and scaled to unit
variance; and Very Nonuniform (VN) distribution, where $\xi_{nr}$ follows
t-distribution with one degree of freedom and scaled to unit variance. Notice
among these three distributions, the score distribution for NU has the
lightest tail, while has the heaviest tail for VN. As a result, the functions
in a synthetic dataset from NU will behave more homogeneously. On the other
hand, the synthetic datasets from MN or VN are likely to contain some
observations with extreme magnitude. Ma, Mahoney and Yu (2015) note
observations in the VN setting could have drastically different leverages,
adversely affecting the naive uniform sampling.

Table 1: Simulation Parameter Specification.

Parameter | Type | Description
---|---|---
Eigenvalue $\\{\sigma^{2}_{r}\\}_{r=1}^{50}$ | Exponential Decay (ED) |  $\sigma^{2}_{r}=c\cdot\kappa^{r}$ with $c=2^{51}$ and $\kappa=0.5$
Polynomial Decay (PD) |  $\sigma^{2}_{r}=c\cdot r^{-\kappa}$ with $c=100$ and $\kappa=1.5$
Score $\xi_{nr}$’s | Nearly Uniform (NU) | standard normal distribution $N(0,1)$
Moderately Nonuniform (MN) | t-distribution with three degree of freedom and scaled to unit variance
| Very Nonuniform (VN) | t-distribution with one degree of freedom and scaled to unit variance

By combining two different choices of eigenvalues (ED or PD) and three
different choices of score distributions (NU, MN, or VN), we have totally six
distinct ways to generate a full sample dataset $\\{x_{n}\\}_{n=1}^{N}$ of
size $N=10,000$. In each replicate of the simulation, one dataset will be
generated in a specified way. The procedure is repeated $1000$ times to report
the average loss and associated standard errors based on the subsample size
$C\in\\{100,300,500,700,1000,3000,5000,7000\\}$. We will apply the randomized
algorithms to estimate the leading functional principal components in Section
5.2 and solve the functional linear regression in Section 5.3. We present here
the results with the ED eigenvalue setting. See Section F in the Appendix for
the results with the PD eigenvalue setting.

Figure 2: Randomized covariance operator and FPCA estimation for the ED and NU
setting. The vertical axes represent comparison metrics in $\log_{10}$ scale.
The blue dotted, red dashed and black solid lines correspond to the results of
the UNIF, IMPO and FunPrinSS sampling, respectively. Figure 3: Randomized
covariance operator and FPCA estimation for the ED and MN setting.  Figure 4:
Randomized covariance operator and FPCA estimation for the ED and VN setting.

###  5.2 Randomized FPCA

To assess the randomized FPCA algorithm, we supply Algorithm 2 with the three
sampling probabilities discussed at the beginning of this section. In
particular, we will estimate the first $R=5$ functional principal components
($\theta_{1},\cdots,\theta_{5}$) and their spanned subspace. A variety of
criteria are used to assess the algorithm performance. The functional
principal subspace estimation error is assessed via the Hilbert-Schmidt norm
and the operator norm
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|_{HS}$ and
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|$ for the functional
principal subspace projection operator. The estimation error for each
principal component function will be computed by
$\|\tilde{\theta}_{r}-\text{sign}(\langle\tilde{\theta}_{r},\hat{\theta}_{r}\rangle)\hat{\theta}_{r}\|$.
As a byproduct, the subsampled covariance operator errors
$\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|_{HS}$ and
$\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|$ are also
reported. Notice these criteria are assessing the differences between the
subsampled estimates and the full sample estimates, as developed in our
theoretical results.

The simulation results with exponentially decaying (ED) eigenvalues are
presented in Figures 2–4. The three figures correspond to the cases with NU,
MN and VN distributions, respectively. In each figure, the first two panels in
the first row demonstrate the estimation error of
$\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|_{HS}$ and
$\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|$. The third and
fourth panels in the first row report the error of functional subspace
projection operator
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|_{HS}$ and
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|$. The estimation
error for the first four eigenfunctions are reported in the second row. In
addition, for each panel, the horizontal axes represent the subsample size
$C$, while the vertical axes are the comparison metrics in the logarithm scale
with base 10. The blue dotted, red dashed and black solid lines correspond to
the results of UNIF, IMPO and FunPrinSS probability, respectively.

From the figures, we can see that as the subsample size increases, all errors
tend to decrease monotonically. There are more points worth emphasizing.
Firstly, UNIF has the worst performance in most of the cases. It almost fails
completely for the VN datasets in Figure 4, where the some principal component
scores have extreme magnitude. Secondly, the sampling strategy IMPO performs
best in estimating $\widehat{\mathscr{C}}_{XX}$ and $\hat{\theta}_{1}$, which
agrees with the result of He and Yan (2020). It also confirms our observations
from Figure 1 in Section 1 that this sampling strategy over-emphasizes the
eigenfunction $\hat{\theta}_{1}$ with the most dominant eigenvalues. As for
the estimation of $\hat{\theta}_{3}$ and $\hat{\theta}_{4}$, the performance
of IMPO is no longer the best. Especially for $\hat{\theta}_{4}$ with a small
eigenvalue, its accuracy can even be worse than UNIF, as shown in the last
penal of Figure 2. Thirdly, our proposed FunPrinSS has the smallest error for
estimating $\widehat{\mathscr{P}}_{R}$, which provides empirical support for
our theoretical results of Theorem 8. It also leads in the accuracy of
estimating $\hat{\theta}_{2},\cdots,\hat{\theta}_{4}$, and the accuracy
remains appealing on the VN dataset with extreme observations in Figure 4. The
above shows the result for the ED eigenvalue setting. The randomized FPCA
results have similar interpretation for the PD eigenvalue setting. See Section
F in the Appendix.

###  5.3 Randomized FLR

Figure 5: Randomized FLR for predictors $x_{n}$ generated with the
exponentially decaying eigenvalues. The blue dotted, red dashed and black
solid lines correspond to the results of the UNIF, IMPO and FunPrinSS
sampling, respectively.

To examine the performance of the proposed randomized FLR algorithm, we
generate data pairs $\\{(x_{n},Y_{n})\\}_{n=1}^{N}$ with $N=10,000$ as
follows. The functional predictor $x_{n}$ is simulated by the procedure in
Section 5.1. The additional scalar response $Y_{n}$ is set by $Y_{n}=\langle
x_{n},\Psi\rangle+\epsilon_{n}$ with noise $\epsilon_{n}\sim N(0,1)$. Besides,
the regression function is specified as
$\Psi(t)=\sum_{r=1}^{50}\theta_{r}(t)$, where $\theta_{r}$’s are the Fourier
bases specified in Section 5.1. Recall by combining different choices of
eigenvalues (ED or PD), and different choices of score distributions (NU, MN,
or VN), six types of simulation settings are considered.

As in Section 5.2, we consider the same sampling probabilities for Algorithm 4
to obtain the estimator $\widetilde{\Psi}$. We consider two metrics to compare
these sampling schemes. The first is the prediction error
$\|\mathscr{T}\widehat{\Psi}-\mathscr{T}\widetilde{\Psi}\|_{N}^{2}=(1/N)\sum_{n=1}^{N}\langle
x_{n},\widehat{\Psi}-\widetilde{\Psi}\rangle^{2}$ discussed in Section 4.2.
The second is direct estimation error
$\|\widehat{\Psi}-\widetilde{\Psi}\|^{2}=\int_{0}^{1}[\widehat{\Psi}(t)-\widetilde{\Psi}(t)]^{2}\,\mathrm{d}t$
in $L_{2}$ norm. The latter is a stronger norm on the difference between the
full sample empirical estimate $\widehat{\Psi}$ and subsampled estimator
$\widetilde{\Psi}$.

The simulation is carried out with 1000 replicates and the operator inverse is
truncated at $R=5$. Figure 5 presents the prediction error and the estimation
error for the ED eigenvalue setting and different score distributions (NU, MN,
VN). The first row depicts the prediction error, while the second row shows
the estimation error. The three columns correspond to the NU, MN and VN score
distribution, respectively. For each panel, the horizontal axis is the
subsample size $C$, and the vertical axis plots the error in the logarithm
scale with base 10. In addition, the blue dotted, red dashed and black solid
lines correspond to the results of UNIF, IMPO and FunPrinSS sampling
probabilities, respectively. It is remarkable that the proposed FunPrinSS has
the smallest errors in all cases and its performance is stable in the VN
datasets with extreme observations. Between different versions of sampling
probabilities, their performance difference enlarges as the tails of the score
distributions grow heavier. In some cases (e.g. when the scores has NU
distribution), the non-informative uniform sampling (UNIF) is slightly better
than the the IMPO sampling. However, UNIF has the worst performance for the VN
score distribution. Figure 14 in Appendix F presents the simulation results
for the PD eigenvalue setting. It also conveys a similar message.

##  6 Real Data Analysis

Figure 6: Five stella spectra from LAMOST with 0.002 offsets in the vertical
direction. The horizontal axis is wavelength (Å) and the vertical axis is
normalized flux.

We now illustrate the effectiveness of our randomized algorithm on a stellar
spectrum dataset. The dataset contains observations from the Large sky Area
Multi-Object fiber Spectroscopic Telescope (LAMOST, Zhao et al., 2012). The
LAMOST project began carrying out its spectroscopic survey of millions of
stars and galaxies in 2012. The current LAMOST data release includes more than
nine million spectra. Aided by this rapid progress in the availability of
spectra data, we can conduct scientific research on many advanced
astrophysical topics, such as exploring galaxy star generation, tracing the
structure and evolution history of our Milky way galaxy and searching for
signatures of dark matter distribution and sub-structures in the Milky Way
halo.

Table 2: Randomized covariance operator and FPCA estimation with $R=5$ for the
stellar spectra dataset. The first column is the comparison metric and the
scale of the reported numbers. In Column 2–10 are the Monte Carlo averages and
the standard errors in the parentheses.

| $C=1,000$ | $C=5,000$ | $C=10,000$
---|---|---|---
| UINF | IMPO | FunPrinSS | UINF | IMPO | FunPrinSS | UINF | IMPO | FunPrinSS
$\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|_{HS}$ $(\times 10^{-9})$ | 20.75 | 4.62 | 12.55 | 9.58 | 2.05 | 5.64 | 7.09 | 1.44 | 3.91
| (0.45) | (0.04) | (0.24) | (0.17) | (0.02) | (0.11) | (0.12) | (0.01) | (0.08)
$\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\|$ $(\times 10^{-9})$ | 19.72 | 3.52 | 12.30 | 9.05 | 1.56 | 5.54 | 6.72 | 1.10 | 3.84
| (0.46) | (0.03) | (0.25) | (0.18) | (0.01) | (0.11) | (0.12) | (0.01) | (0.08)
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|_{HS}$ $(\times 10^{-1})$ | 11.65 | 6.45 | 2.31 | 7.11 | 2.87 | 1.03 | 5.35 | 2.05 | 0.73
| (0.11) | (0.10) | (0.04) | (0.13) | (0.05) | (0.01) | (0.12) | (0.03) | (0.01)
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|$ $(\times 10^{-1})$ | 7.72 | 4.34 | 1.52 | 4.76 | 1.93 | 0.68 | 3.59 | 1.37 | 0.48
| (0.08) | (0.07) | (0.03) | (0.10) | (0.04) | (0.01) | (0.08) | (0.03) | (0.01)
$\|\tilde{\theta}_{1}-\text{sign}(\langle\tilde{\theta}_{1},\hat{\theta}_{1}\rangle)\hat{\theta}_{1}\|$ $(\times 10^{-2})$ | 1.79 | 0.77 | 0.91 | 0.89 | 0.34 | 0.40 | 0.64 | 0.24 | 0.28
| (0.04) | (0.01) | (0.01) | (0.01) | (0.00) | (0.00) | (0.01) | (0.00) | (0.00)
$\|\tilde{\theta}_{2}-\text{sign}(\langle\tilde{\theta}_{2},\hat{\theta}_{2}\rangle)\hat{\theta}_{2}\|$ $(\times 10^{-2})$ | 21.70 | 7.94 | 5.15 | 7.70 | 3.54 | 2.31 | 5.35 | 2.49 | 1.62
| (0.81) | (0.09) | (0.07) | (0.10) | (0.04) | (0.03) | (0.06) | (0.03) | (0.02)
$\|\tilde{\theta}_{3}-\text{sign}(\langle\tilde{\theta}_{3},\hat{\theta}_{3}\rangle)\hat{\theta}_{3}\|$ $(\times 10^{-2})$ | 39.39 | 12.37 | 6.46 | 16.69 | 4.53 | 2.88 | 9.55 | 3.20 | 1.99
| (1.40) | (0.36) | (0.07) | (0.86) | (0.05) | (0.03) | (0.52) | (0.03) | (0.02)
$\|\tilde{\theta}_{4}-\text{sign}(\langle\tilde{\theta}_{4},\hat{\theta}_{4}\rangle)\hat{\theta}_{4}\|$ $(\times 10^{-2})$ | 99.42 | 34.68 | 9.91 | 53.67 | 14.53 | 4.29 | 37.89 | 9.93 | 3.03
| (0.92) | (0.83) | (0.20) | (1.10) | (0.35) | (0.07) | (0.91) | (0.23) | (0.05)
$\|\tilde{\theta}_{5}-\text{sign}(\langle\tilde{\theta}_{5},\hat{\theta}_{5}\rangle)\hat{\theta}_{5}\|$ $(\times 10^{-2})$ | 90.89 | 55.06 | 18.11 | 61.87 | 24.30 | 7.93 | 47.71 | 17.08 | 5.62
| (1.10) | (0.99) | (0.29) | (1.28) | (0.46) | (0.11) | (1.15) | (0.29) | (0.08)

From the the current LAMOST data release, we randomly select a dataset
containing $N=110,000$ stellar spectra with $r$-band signal-to-noise ratio
larger than 35. Each spectrum is normalized to have unit norm and get denoised
by a spline wavelet method. In addition, the average spectrum is subtracted
from each observation such that the dataset has zero mean. We treat this as
our full dataset and apply the subsampling algorithms over this dataset.
Figure 6 presents five example spectra from the pre-processed dataset. In the
figure, the horizontal axis is wavelength, and the vertical axis is flux which
measures the brightness of an object. Each spectrum observation can be viewed
as functional data, where the flux value $x_{n}(\lambda)$ varies as a function
of the wavelength $\lambda$. In Figure 6, several important spectral features
are marked as vertical dashed lines. These spectral features are adopted by
Liu et al. (2015) for spectral classification.

Characterizing the diversity of stars and extracting stellar atmospheric
parameters are of key astrophysical importance for astronomical observational
studies. We will apply the randomized FPCA in Algorithm 2 and the randomized
FLR in Algorithm 4 to this dataset.

###  6.1 Results for Randomized FPCA

Dimensional reduction has been an effective tools for astronomical spectrum
data analysis (Bermejo, Ramos and Prieto, 2013; McGurk, Kimball and Ivezic,
2010). We apply the randomized FPCA in Algorithm 2 with the three sampling
strategies (UNIF, IMPO and FunPrinSS) in Section 5. We will estimate the
leading $R=5$ functional principal components, together with the covariance
operator. The parameter $R=5$ is chosen such that FVE is approximately $98\%$.
The improvement of FVE with $R=6$ is marginal (see Figure 7).

Based on the whole selected 110,000 spectra, the full sample estimators are
computed with brute force. Besides, we apply the randomized algorithm with the
subsample size $C$ varying among $1,000$, $5,000$ and $10,000$. The randomized
estimates are compared against the full sample results by the metrics
introduced in Section 5.2. For each subsample size, the procedure is repeated
$1,000$ times to compute the average estimation loss. Table 2 presents the
resulting average accuracy and related standard errors. The method performance
in Table 2 can be interpreted similarly as that from Section 5.2. The IMPO
sampling strategy delivers the smallest error in estimating $\hat{\theta}_{1}$
and $\widehat{\mathscr{C}}_{XX}$. However, FunPrinSS has the best performance
in estimating the projection operator $\widehat{\mathscr{P}}_{R}$ and the
other functional principal components.

Figure 7: The boxplot of the FVE of the subsampled functional principal
component estimates under 1000 replicates. The three panels correspond to the
subsample size $C=1000$, $5000$ and $10000$, respectively. In each panel,
different sampling probabilities (UNIF, IMPO, FunPrinSS) are compared over
various $R$. The FVEs of the full sample estimator are shown as dashed
horizontal lines.  Figure 8:  The green solid lines (in the first row)
represent the leading five eigenfunctions of full sample covariance operator.
The blue, red and black solid lines (in the second row) represent the averaged
absolute difference between the full sample eigenfunctions and the subsampled
eigenfunctions obtained by UNIF, IMPO and FunPrinSS sampling, respectively.

In Table 2, for the randomized FPCA with the FunPrinSS probability, we can see
the Monte Carlo average error
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|\approx 0.152$,
$0.068$, $0.048$ for $C=1000$, $5000$ and $10000$, respectively. Neglecting
the second order term ($\epsilon^{2}$) in (3.5) of Theorem 8, we can find the
expected subsample error
$\mathbb{E}\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|$ is upper
bounded by a first-order term with magnitude order
$\mathscr{O}\big{(}(R+\Delta_{R})\log^{1/2}(R+\Delta_{R})/C^{1/2}\big{)}$. For
$C=1000$, $5000$ and $10000$, the values of this theoretical order are
$\mathscr{O}\big{(}0.469)$, $\mathscr{O}\big{(}0.210)$ and
$\mathscr{O}\big{(}0.148)$, respectively. The ratio between the Monte Carlo
average error and the theoretical first-order approximate magnitude is stable
at a constant around $0.32$.

In Figure 8, we plot the leading five full sample eigenfunctions in the first
row of panels, and the averaged absolute difference between the full sample
and subsampled eigenfunctions in the second row of panels. The absolute
difference $|\hat{\theta}_{r}(\lambda)-\tilde{\theta}_{r}(\lambda)|$ at each
wavelength $\lambda$ is averaged over the $1000$ replicates when the subsample
size is $C=1,000$. From the figure, we can compare the ability of capturing
the spectral features ($H_{\beta}$, $Na$, $Mg$ and $H_{\alpha}$ marked by the
vertical dashed lines) by various sampling probabilities. As expected, the
UNIF sampling (colored by blue) has the largest error for characterizing the
eigenfunctions and spectral features. The IMPO probability (red) is
competitive in recovering spectral lines for the first eigenfunctions but not
so for the rest of eigenfunctions. For the other eigenfunctions, the error of
estimating the spectral features by the FunPrinSS probability (black) is the
smallest.

We examine the accuracy of the FVE of the subsampled estimator (see
$\widetilde{\mathrm{FVE}}$ in Section 3.4). Figure 7 draws the boxplots for
the FVE of the subsampled functional principal component estimates under 1000
replicates. The three panels correspond to the subsample size $C=1000$, $5000$
and $10000$, respectively. In each panel, the sampling probabilities (UNIF,
IMPO, FunPrinSS) are compared over various $R$. The FVEs of the full sample
estimator are shown as dashed horizontal lines. From the figure, we can see
the FVE of the subsampled estimator serves as an estimate of the FVE of the
full sample estimator. In particular, the FVE of the subsampled estimator with
the FunPrinSS probability is very accurate with small variance.

###  6.2 Results for Randomized FLR

Figure 9: Subsampling error in estimation and in prediction for the randomized
FLR for the LAMOST spectra dataset. The blue dotted, red dashed and black
solid lines correspond to the results of the UNIF, IMPO and FunPrinSS
sampling, respectively. The gray solid curves in the first row show the
magnitude
$\mathscr{O}\big{(}(R+\Delta_{R})(\|\mathbf{Y}\|_{N}+\|\mathbf{Y}^{\perp}\|_{N})\log^{1/2}(R+\Delta_{R})/C^{1/2}\big{)}$
of the first-order upper bound for the expected subsample error. Figure 10:
The comparison of $\widetilde{\Psi}$ obtained by the randomized FLR with
different sampling probabilities. The first row (green lines) presents the
full sample regression functions $\widehat{\Psi}$ and the second row shows the
average absolute difference between the full sample and the subsampled
estimates. The blue, red and black solid lines (in the second row) correspond
to the subsampled estimator $\widetilde{\Psi}$’s obtained by UNIF, IMPO and
FunPrinSS, respectively. The left, middle and right column have the results
for the $\log\,g$, $T_{\text{eff}}$ and [Fe/H] response, respectively.

The fundamental stellar atmospheric parameters such as effective temperature
($T_{\text{eff}}$), surface gravity (log $g$), and metallicity ([Fe/H]) are
useful probes for revealing the kinematic and chemical properties of stars. We
conduct three functional linear regression tasks to compare the different
sampling schemes (UNIF, IMPO and FunPrinSS). For each task, the response
$Y_{n}$ is selected as one of the stellar atmospheric parameters: the
effective temperature ($T_{\text{eff}}$), the surface gravity ($\log g$), and
the metallicity ([Fe/H]). In particular, the regression model is
$Y_{n}=\int_{\Lambda}x_{n}(\lambda)\Psi(\lambda)\,\mathrm{d}\lambda+\epsilon_{n}$.
Each spectrum $x_{n}$ is integrated with the regression function $\Psi$ over
the observed wavelength range. Both the response and the spectrum predictor
are centered with zero mean, so we do not need to estimate the intercept term.

The accuracy metrics in Section 5.3 are also employed here for comparison
between the full sample and subsampled estimator. The full sample estimator
$\widehat{\Psi}$ is obtained via brute force computation on the full dataset.
In Figure 9, the subsampling prediction error and estimation error are shown
in the first and second row, respectively. The three columns of Figure 9
corresponds to the results with the response $T_{\text{eff}}$, $\log g$, and
[Fe/H], respectively. In each penal, the vertical axis is the error in the
scale of $\log_{10}$; and the horizontal axis is the subsample size $C$. These
plots suggest that the proposed FunPrinSS method outperforms the other two
sampling schemes under all circumstances. Theorem 14 indicates the expected
subsample error
$\mathbb{E}\|\mathscr{T}\widetilde{\Psi}-\mathscr{T}\widehat{\Psi}\|_{N}$ can
be approximately upper bounded by a first-order term with magnitude order
$\mathscr{O}\big{(}(R+\Delta_{R})(\|\mathbf{Y}\|_{N}+\|\mathbf{Y}^{\perp}\|_{N})\log^{1/2}(R+\Delta_{R})/C^{1/2}\big{)}$.
The values of the theoretical order are plotted as the gray solid curves in
the first row of Figure 9. The vertical axis is in the scale of $\log_{10}$,
and the theoretical order predicts the Monte Carlo performance up to an
offset. Note the theoretical order underpredicts the actual Monte Carlo
average error for small $C$, as the higher order error terms in Theorem 14
have been discarded.

Figure 10 plots the full sample regression function estimates $\widehat{\Psi}$
in the first row. The second row shows the averaged absolute difference
between the full sample estimates $\widehat{\Psi}$ and the subsampled
estimates $\widetilde{\Psi}$ over $1000$ replicates with $C=1,000$. The work
of Li et al. (2014) has identified a few spectral features in determining the
stellar atmosphere parameters. Some of these key spectral features are also
marked by the vertical dashed lines. The proposed FunPrinSS has generally
smaller error in identifying these spectral features (e.g. $H_{\alpha}$
$\lambda 6563$) and provides a better estimation of the feature strength.

## Acknowledgment

The authors thank the editor, the associate editor, and an anonymous reviewer
for constructive comments that helped significantly improve this work. The
authors also thank Professor Jianhua Z. Huang, Professor Xiangyu Chang and
Professor Yixuan Qiu for comments on the initial draft of this work. Shiyuan
He’s research was supported by National Natural Science Foundation of China
(No.11801561). This research was supported by Public Computing Cloud, Renmin
University of China.

## References

  * Bermejo, Ramos and Prieto (2013) [author] Bermejo, J MuñozJ. M., Ramos, A AsensioA. A. and Prieto, C AllendeC. A. (2013). A PCA approach to stellar effective temperatures. Astronomy & Astrophysics 553 A95.
  * Cardot, Degras and Josserand (2013) [author] Cardot, HervéH., Degras, DavidD. and Josserand, EtienneE. (2013). Confidence bands for Horvitz–Thompson estimators using sampled noisy functional data. Bernoulli 19 2067–2097.
  * Cardot, Goga and Lardin (2013) [author] Cardot, HervéH., Goga, CameliaC. and Lardin, PaulineP. (2013). Uniform convergence and asymptotic confidence bands for model-assisted estimators of the mean of sampled functional data. Electronic journal of statistics 7 562–596.
  * Cardot and Josserand (2011) [author] Cardot, HervéH. and Josserand, EtienneE. (2011). Horvitz–Thompson estimators for functional data: Asymptotic confidence bands and optimal allocation for stratified sampling. Biometrika 98 107–118.
  * Cardot, Mas and Sarda (2007) [author] Cardot, HervéH., Mas, AndréA. and Sarda, PascalP. (2007). CLT in functional linear regression models. Probability Theory and Related Fields 138 325–361.
  * Cardot et al. (2010) [author] Cardot, HervéH., Chaouch, MohamedM., Goga, CameliaC. and Labruère, CatherineC. (2010). Properties of design-based functional principal components analysis. Journal of statistical planning and inference 140 75–91.
  * Connolly et al. (1994) [author] Connolly, Andy JA. J., Szalay, ASA., Bershady, MAM., Kinney, ALA. and Calzetti, DD. (1994). Spectral classification of galaxies: an orthogonal approach. arXiv preprint astro-ph/9411044.
  * Degras (2014) [author] Degras, DavidD. (2014). Rotation sampling for functional data. Statistica Sinica 1075–1095.
  * Delaigle et al. (2010) [author] Delaigle, AuroreA., Hall, PeterP. et al. (2010). Defining probability density for a distribution of random functions. The Annals of Statistics 38 1171–1193.
  * Dicker et al. (2017) [author] Dicker, Lee HL. H., Foster, Dean PD. P., Hsu, DanielD. et al. (2017). Kernel ridge vs. principal component regression: Minimax bounds and the qualification of regularization operators. Electronic Journal of Statistics 11 1022–1047.
  * Drineas, Kannan and Mahoney (2006a) [author] Drineas, PetrosP., Kannan, RaviR. and Mahoney, Michael WM. W. (2006a). Fast Monte Carlo algorithms for matrices I: Approximating matrix multiplication. SIAM Journal on Computing 36 132–157.
  * Drineas, Kannan and Mahoney (2006b) [author] Drineas, PetrosP., Kannan, RaviR. and Mahoney, Michael WM. W. (2006b). Fast Monte Carlo algorithms for matrices II: Computing a low-rank approximation to a matrix. SIAM Journal on computing 36 158–183.
  * Drineas, Kannan and Mahoney (2006c) [author] Drineas, PetrosP., Kannan, RaviR. and Mahoney, Michael WM. W. (2006c). Fast Monte Carlo algorithms for matrices III: Computing a compressed approximate matrix decomposition. SIAM Journal on Computing 36 184–206.
  * Drineas, Mahoney and Muthukrishnan (2006) [author] Drineas, PetrosP., Mahoney, Michael W.M. W. and Muthukrishnan, S.S. (2006). Subspace Sampling and Relative-Error Matrix Approximation: Column-Row-Based Methods. In Algorithms - ESA 2006, 14th Annual European Symposium, Zurich, Switzerland, September 11-13, 2006, Proceedings.
  * Drineas and Mahoney (2018) [author] Drineas, PetrosP. and Mahoney, Michael WM. W. (2018). Lectures on randomized numerical linear algebra. The Mathematics of Data 25 1\.
  * Drineas et al. (2011) [author] Drineas, PetrosP., Mahoney, Michael WM. W., Muthukrishnan, SS. and Sarlós, TamásT. (2011). Faster least squares approximation. Numerische mathematik 117 219–249.
  * Drineas et al. (2012) [author] Drineas, PetrosP., Magdon-Ismail, MalikM., Mahoney, Michael WM. W. and Woodruff, David PD. P. (2012). Fast approximation of matrix coherence and statistical leverage. Journal of Machine Learning Research 13 3475–3506.
  * Eisenstein et al. (2011) [author] Eisenstein, Daniel JD. J., Weinberg, David HD. H., Agol, EricE., Aihara, HiroakiH., Prieto, Carlos AllendeC. A., Anderson, Scott FS. F., Arns, James AJ. A., Aubourg, ÉricÉ., Bailey, StephenS., Balbinot, EduardoE. et al. (2011). SDSS-III: Massive spectroscopic surveys of the distant universe, the Milky Way, and extra-solar planetary systems. The Astronomical Journal 142 72\.
  * Hadjipantelis and Müller (2018) [author] Hadjipantelis, Pantelis ZenonP. Z. and Müller, Hans-GeorgH.-G. (2018). Functional data analysis for big data: A case study on california temperature trends. In Handbook of Big Data Analytics 457–483. Springer.
  * Halko, Martinsson and Tropp (2011) [author] Halko, NathanN., Martinsson, PergunnarP. and Tropp, Joel AJ. A. (2011). Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions. Siam Review 53 217–288.
  * Hall et al. (2007) [author] Hall, PeterP., Horowitz, Joel LJ. L. et al. (2007). Methodology and convergence rates for functional linear regression. The Annals of Statistics 35 70–91.
  * He and Yan (2020) [author] He, ShiyuanS. and Yan, XiaomengX. (2020). Randomized estimation of functional covariance operator via subsampling. Stat 9 e311.
  * Horváth and Kokoszka (2012) [author] Horváth, LajosL. and Kokoszka, PiotrP. (2012). Inference for functional data with applications 200\. Springer Science & Business Media.
  * Horváth, Kokoszka and Reeder (2013) [author] Horváth, LajosL., Kokoszka, PiotrP. and Reeder, RonR. (2013). Estimation of the mean of functional time series and a two-sample problem. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 75 103–122.
  * Hsing and Eubank (2015) [author] Hsing, TailenT. and Eubank, RandallR. (2015). Theoretical foundations of functional data analysis, with an introduction to linear operators. John Wiley & Sons.
  * James, Hastie and Sugar (2000) [author] James, Gareth MG. M., Hastie, Trevor JT. J. and Sugar, Catherine AC. A. (2000). Principal component models for sparse functional data. Biometrika 87 587–602.
  * Kato (2012) [author] Kato, KengoK. (2012). Estimation in functional linear quantile regression. Annals of Statistics 40 3108–3136.
  * Kokoszka and Reimherr (2017) [author] Kokoszka, PiotrP. and Reimherr, MatthewM. (2017). Introduction to functional data analysis. CRC Press.
  * Koltchinskii et al. (2016) [author] Koltchinskii, VladimirV., Lounici, KarimK. et al. (2016). Asymptotics and concentration bounds for bilinear forms of spectral projectors of sample covariance. In Annales de l’Institut Henri Poincaré, Probabilités et Statistiques 52 1976–2013. Institut Henri Poincaré.
  * Lardin-Puech, Cardot and Goga (2014) [author] Lardin-Puech, PaulineP., Cardot, HervéH. and Goga, CameliaC. (2014). Analysing large datasets of functional data: a survey sampling point of view. Journal de la Société Française de Statistique 155 70–94.
  * Li et al. (2014) [author] Li, XiangruX., Wu, QM JonathanQ. J., Luo, AliA., Zhao, YonghengY., Lu, YuY., Zuo, FangF., Yang, TanT. and Wang, YongjunY. (2014). SDSS/SEGUE spectral feature analysis for stellar atmospheric parameter estimation. The Astrophysical Journal 790 105\.
  * Liu et al. (2015) [author] Liu, ChaoC., Cui, Wen-YuanW.-Y., Zhang, BoB., Wan, Jun-ChenJ.-C., Deng, Li-CaiL.-C., Hou, Yong-HuiY.-H., Wang, Yue-FeiY.-F., Yang, MingM. and Zhang, YongY. (2015). Spectral classification of stars based on LAMOST spectra. Research in Astronomy and Astrophysics 15 1137\.
  * Liu et al. (2018) [author] Liu, PengP., Di, LipingL., Du, QianQ. and Wang, LizheL. (2018). Remote Sensing Big Data: Theory, Methods and Applications. Remote Sensing 10 711\.
  * Ma, Mahoney and Yu (2015) [author] Ma, PingP., Mahoney, Michael WM. W. and Yu, BinB. (2015). A statistical perspective on algorithmic leveraging. The Journal of Machine Learning Research 16 861–911.
  * McGurk, Kimball and Ivezic (2010) [author] McGurk, Rosalie C.R. C., Kimball, Amy E.A. E. and Ivezic, ZeljkoZ. (2010). Principal Component Analysis of SDSS Stellar Spectra. Astron. J. 139 1261–1268. 10.1088/0004-6256/139/3/1261
  * Minsker (2017) [author] Minsker, StanislavS. (2017). On some extensions of Bernstein’s inequality for self-adjoint operators. Statistics & Probability Letters 127 111–119.
  * Mor-Yosef and Avron (2019) [author] Mor-Yosef, LironL. and Avron, HaimH. (2019). Sketching for principal component regression. SIAM Journal on Matrix Analysis and Applications 40 454–485.
  * Peng and Paul (2009) [author] Peng, JieJ. and Paul, DebashisD. (2009). A geometric approach to maximum likelihood estimation of the functional principal components from sparse longitudinal data. Journal of Computational and Graphical Statistics 18 995–1015.
  * Pilanci and Wainwright (2015) [author] Pilanci, MertM. and Wainwright, Martin JM. J. (2015). Randomized sketches of convex programs with sharp guarantees. IEEE Transactions on Information Theory 61 5096–5115.
  * Pilanci and Wainwright (2017) [author] Pilanci, MertM. and Wainwright, Martin JM. J. (2017). Newton sketch: A near linear-time optimization algorithm with linear-quadratic convergence. SIAM Journal on Optimization 27 205–245.
  * Ramsay (2004) [author] Ramsay, James OJ. O. (2004). Functional data analysis. Encyclopedia of Statistical Sciences 4\.
  * Raskutti and Mahoney (2016) [author] Raskutti, GarveshG. and Mahoney, Michael WM. W. (2016). A statistical perspective on randomized sketching for ordinary least-squares. Journal of Machine Learning Research 17 7508–7538.
  * Tian (2010) [author] Tian, Tian SivaT. S. (2010). Functional Data Analysis in Brain Imaging Studies. Frontiers in Psychology 1 35–35.
  * Tropp et al. (2015) [author] Tropp, Joel AJ. A. et al. (2015). An introduction to matrix concentration inequalities. Foundations and Trends® in Machine Learning 8 1–230.
  * Turkbrowne (2013) [author] Turkbrowne, Nicholas BN. B. (2013). Functional interactions as big data in the human brain. Science 342 580–584.
  * Wang (2019) [author] Wang, HaiYingH. (2019). More Efficient Estimation for Logistic Regression with Optimal Subsamples. Journal of Machine Learning Research 20 1–59.
  * Wang, Gittens and Mahoney (2017) [author] Wang, ShusenS., Gittens, AlexA. and Mahoney, Michael WM. W. (2017). Sketched ridge regression: Optimization perspective, statistical perspective, and model averaging. The Journal of Machine Learning Research 18 8039–8088.
  * Wang, Zhu and Ma (2018) [author] Wang, HaiYingH., Zhu, RongR. and Ma, PingP. (2018). Optimal subsampling for large sample logistic regression. Journal of the American Statistical Association 113 829–844.
  * Woodruff (2014) [author] Woodruff, David PD. P. (2014). Sketching as a tool for numerical linear algebra. arXiv preprint arXiv:1411.4357.
  * Yang et al. (2017) [author] Yang, YunY., Pilanci, MertM., Wainwright, Martin JM. J. et al. (2017). Randomized sketches for kernels: Fast and optimal nonparametric regression. The Annals of Statistics 45 991–1023.
  * Yao, Müller and Wang (2005) [author] Yao, FangF., Müller, Hans-GeorgH.-G. and Wang, Jane-LingJ.-L. (2005). Functional data analysis for sparse longitudinal data. Journal of the American Statistical Association 100 577–590.
  * Yao et al. (2005) [author] Yao, FangF., Müller, Hans-GeorgH.-G., Wang, Jane-LingJ.-L. et al. (2005). Functional linear regression analysis for longitudinal data. The Annals of Statistics 33 2873–2903.
  * Yuan and Cai (2010) [author] Yuan, MingM. and Cai, T TonyT. T. (2010). A reproducing kernel Hilbert space approach to functional linear regression. The Annals of Statistics 38 3412–3444.
  * Zhang, Ning and Ruppert (2021) [author] Zhang, TaoT., Ning, YangY. and Ruppert, DavidD. (2021). Optimal sampling for generalized linear models under measurement constraints. Journal of Computational and Graphical Statistics 30 106–114.
  * Zhao et al. (2012) [author] Zhao, GangG., Zhao, Yong-HengY.-H., Chu, Yao-QuanY.-Q., Jing, Yi-PengY.-P. and Deng, Li-CaiL.-C. (2012). LAMOST spectral survey—An overview. Research in Astronomy and Astrophysics 12 723\.

##  Appendix A Comparison with sampling probabilities in the multivariate
context

There are some underlying connections between the FunPrinSS probability (2.2)
for functional data and some sampling probabilities for multivariate
statistics in the literature. In this section, we discuss some of the
similarities and differences.

The FunPrinSS probability (2.2) is closely related to the leverage sampling
probability (Drineas et al., 2012) for classical linear regression. Let
$\mathbf{X}=\left(\mathbf{x}_{1},\ldots,\mathbf{x}_{N}\right)^{T}$ be a
$N\times p$ fixed design matrix with $N$ observations and $p$ variables (with
$N>p$). Suppose the singular value decomposition (SVD) of $\mathbf{X}$ is
$\mathbf{X}=\mathbf{U}\mathbf{D}\mathbf{V}^{T}$, where $\mathbf{D}$ is a
diagonal matrix with singular values, and the columns of $\mathbf{U}$ and
$\mathbf{V}$ contain the left and right singular vectors, respectively. Denote
$U_{nj}$ as the $(n,j)$-th element of $\mathbf{U}$. The leverage sampling
probability $p_{n}$ is proportional to the statistical leverage $h_{nn}$ of
$n$-th observation, i.e. $p_{n}\propto
h_{nn}=\mathbf{x}_{n}^{T}\left(\mathbf{X}^{T}\mathbf{X}\right)^{-}\mathbf{x}_{n}=\sum_{j=1}^{p}U_{nj}^{2}.$
When the samples $\mathbf{x}_{n}$’s have zero mean, $U_{nj}$ can be
interpreted as multivariate principal component scores, and the statistical
leverage $h_{nn}$ is the summation of all squared scores in the
$p$-dimensional space. In contrast, our FunPrinSS probability (2.2) makes a
rank-$R$ truncation to the Karhunen-Loève expansion (2.1) of the infinite
dimensional functional data, and takes into account the effect of the
remaining scores via the projected residual
$(I-\widehat{\mathscr{P}}_{R})x_{n}$.

The FunPrinSS probability is also connected to the subspace sampling of
Drineas, Mahoney and Muthukrishnan (2006). To perform the CUR decomposition
for a matrix $\mathbf{X}\in\mathbb{R}^{N\times p}$, they compute the row
sampling probabilities as follows. A matrix $\mathbf{C}\in\mathbb{R}^{N\times
q}$ is formed by sampling columns from $\mathbf{X}$, and its SVD
$\mathbf{C}=\mathbf{U}_{C}\mathbf{D}_{C}\mathbf{V}_{C}^{T}$ is computed. Then,
the sampling probability for $n$-th row of $\mathbf{X}$ is computed from the
$n$-th row of $\mathbf{U}_{C}$, the $n$-th row of
$\mathbf{X}^{\perp}_{C}=(\mathbf{I}-\mathbf{U}_{C}\mathbf{U}_{C}^{T})\mathbf{X}$
and their product

$p_{n}=\frac{(1/3)\|(\mathbf{U}_{C})_{(n)}\|_{2}^{2}}{\sum_{m=1}^{N}\|(\mathbf{U}_{C})_{(m)}\|_{2}^{2}}+\frac{(1/3)\|(\mathbf{U}_{C})_{(n)}\|_{2}\|(\mathbf{X}^{\perp}_{C})_{(n)}\|_{2}}{\sum_{m=1}^{N}\|(\mathbf{U}_{C})_{(m)}\|_{2}\|(\mathbf{X}^{\perp}_{C})_{(m)}\|_{2}}+\frac{(1/3)\|(\mathbf{X}^{\perp}_{C})_{(n)}\|_{2}^{2}}{\sum_{m=1}^{N}\|(\mathbf{X}^{\perp}_{C})_{(m)}\|_{2}^{2}}.$

In the above, $(\mathbf{A})_{(n)}$ represents the $n$-th row of a matrix
$\mathbf{A}$ and $\|\cdot\|_{2}$ is the classical unnormalized Euclidean norm.
When we view the rows of $\mathbf{X}$ as samples (observations) of a zero mean
random vector and the columns of $\mathbf{X}$ as variables, the matrix
$\mathbf{U}_{C}$ serves as an estimate of the principal component scores, and
$(\mathbf{I}-\mathbf{U}_{C}\mathbf{U}_{C}^{T})\mathbf{X}$ is the residual
matrix.

##  Appendix B Concentration Inequality for Compact Operators

Theorem 7.3.1 of Tropp et al. (2015) developed intrinsic matrix Bernstein
concentration inequality. Together with the techniques in Section 3.2 of
Minsker (2017), it can be easily extended to compact operators. Similar
results can also be found in Lemma 5 of Dicker et al. (2017). In the below,
$\mathscr{H}_{1}$ and $\mathscr{H}_{2}$ are two Hilbert spaces.

######  Lemma 15.

For a finite sequence of random operators $\mathscr{Z}_{i}$ mapping from
$\mathscr{H}_{1}$ to $\mathscr{H}_{2}$, they satisfies
$\mathbb{E}\mathscr{Z}_{i}=0$ and $\|\mathscr{Z}_{i}\|\leq L$. Suppose their
summation is $\mathscr{S}=\sum_{i}\mathscr{Z}_{i}$ and $\mathscr{S}^{*}$ is
its adjoint. Let $\mathscr{V}_{1}$ and $\mathscr{V}_{2}$ be semidefinite upper
bounds for $\mathrm{Var}_{1}(\mathscr{S})$ and
$\mathrm{Var}_{2}(\mathscr{S})$, respectively. That is,
$\mathscr{V}_{1}\succeq\mathrm{Var}_{1}(\mathscr{S})=\mathbb{E}(\mathscr{S}\mathscr{S}^{*})=\sum_{i}\mathbb{E}(\mathscr{Z}_{i}\mathscr{Z}_{i}^{*})$,
and
$\mathscr{V}_{2}\succeq\mathrm{Var}_{2}(\mathscr{S})=\mathbb{E}(\mathscr{S}^{*}\mathscr{S})=\sum_{i}\mathbb{E}(\mathscr{Z}_{i}^{*}\mathscr{Z}_{i})$.
Define an intrinsic bound
$d=\mathrm{intdim}(\mathscr{V}_{1})+\mathrm{intdim}(\mathscr{V}_{2})$ and a
variance bound $v=\max\\{\|\mathscr{V}_{1}\|,\|\mathscr{V}_{2}\|\\}$, then for
$t\geq\sqrt{v}+L/3$,

$\mathbb{P}(\|\mathscr{S}\|\geq t)\leq
4d\exp\left(-\frac{t^{2}/2}{v+Lt/3}\right).$

The next lemma provides upper bound for self-adjoint operator.

######  Lemma 16.

For a finite sequence of random self-adjoint operators $\mathscr{Z}_{i}$
mapping from $\mathscr{H}_{1}$ to $\mathscr{H}_{1}$, they satisfies
$\mathbb{E}\mathscr{Z}_{i}=0$ and $\|\mathscr{Z}_{i}\|\leq L$. Suppose their
summation is $\mathscr{S}=\sum_{i}\mathscr{Z}_{i}$, and $\mathscr{V}$ is the
semidefinite upper bound for $\sum_{i}\mathbb{E}(\mathscr{Z}_{i}^{2})$. Define
an intrinsic bound $d=\mathrm{intdim}(\mathscr{V})$ and a variance bound
$v=\|\mathscr{V}\|$, then for $t\geq\sqrt{v}+L/3$,

$\mathbb{P}(\|\mathscr{S}\|\geq t)\leq
4d\exp\left(-\frac{t^{2}/2}{v+Lt/3}\right).$

##  Appendix C Proof for Section 3

###  C.1 Proof of Lemma 6

Define $\widehat{\mathscr{P}}_{r}=\hat{\theta}_{r}\otimes\hat{\theta}_{r}$.
With the definition of the resolvent
$\mathscr{R}_{\widehat{\mathscr{C}}_{XX}}(\eta)=(\widehat{\mathscr{C}}_{XX}-\eta
I)^{-1}=-\sum_{r=1}^{\infty}\frac{1}{\eta-\hat{\sigma}_{r}^{2}}\widehat{\mathscr{P}}_{r}$,
the linear error term can be expanded as

$\displaystyle L_{R}(\mathscr{E})$ $\displaystyle=\frac{1}{2\pi
i}\oint_{\Gamma_{R}}\mathscr{R}_{\widehat{\mathscr{C}}_{XX}}(\eta)\mathscr{E}\mathscr{R}_{\widehat{\mathscr{C}}_{XX}}(\eta)\mathrm{d}\eta=\sum_{r=1}^{R}\sum_{s=R+1}^{\infty}\frac{1}{\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}}\left(\widehat{\mathscr{P}}_{r}\mathscr{E}\widehat{\mathscr{P}}_{s}+\widehat{\mathscr{P}}_{s}\mathscr{E}\widehat{\mathscr{P}}_{r}\right).$

by Cauchy integration formula. The above results can also be deduced from
Theorem 5.1.4 of Hsing and Eubank (2015). Notice that

$\displaystyle\widehat{\mathscr{P}}_{r}(x_{n}\otimes
x_{n})\widehat{\mathscr{P}}_{s}$
$\displaystyle=(\hat{\theta}_{r}\otimes\hat{\theta}_{r})(x_{n}\otimes
x_{n})(\hat{\theta}_{s}\otimes\hat{\theta}_{s})$
$\displaystyle=\langle\hat{\theta}_{r},x_{n}\rangle\cdot\langle\hat{\theta}_{s},x_{n}\rangle\cdot(\hat{\theta}_{s}\otimes\hat{\theta}_{r})=\hat{\sigma}_{r}\hat{\sigma}_{s}\hat{\xi}_{ns}\hat{\xi}_{nr}(\hat{\theta}_{s}\otimes\hat{\theta}_{r}),$

and
$\widehat{\mathscr{P}}_{r}\widehat{\mathscr{C}}_{XX}\widehat{\mathscr{P}}_{s}=0$
for any $r\neq s$. From these, we know for $r\neq s$ that

$\displaystyle\widehat{\mathscr{P}}_{r}\mathscr{E}\widehat{\mathscr{P}}_{s}=\widehat{\mathscr{P}}_{r}(\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX})\widehat{\mathscr{P}}_{s}=\frac{\hat{\sigma}_{r}\hat{\sigma}_{s}}{CN}\sum_{c=1}^{C}\frac{\tilde{\xi}_{cr}\tilde{\xi}_{cs}}{\tilde{p}_{c}}(\hat{\theta}_{s}\otimes\hat{\theta}_{r}).$

In the above, $\tilde{\xi}_{cs}$ is the $s$-th score of the subsampled
$\tilde{x}_{c}$. It holds
$\tilde{\xi}_{cs}=\langle\tilde{x}_{c},\hat{\theta}_{s}\rangle/\hat{\sigma}_{s}$
when $\hat{\sigma}_{s}>0$; and $\tilde{\xi}_{cs}=1$ when $\hat{\sigma}_{s}=0$.
Notice it is computed from the full sample eigenfunction $\hat{\theta}_{s}$
and the subsampled $\tilde{x}_{c}$. As a result, the linear error term can be
expressed as
$L_{R}(\mathscr{E})=\frac{1}{C}\sum_{c=1}^{C}\frac{\mathscr{Z}_{c}}{N\tilde{p}_{c}}$,
with

$\mathscr{Z}_{c}=\sum_{r=1}^{R}\sum_{s=R+1}^{\infty}\frac{\hat{\sigma}_{r}\hat{\sigma}_{s}}{\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}}\times\tilde{\xi}_{cr}\tilde{\xi}_{cs}\times\big{[}\hat{\theta}_{r}\otimes\hat{\theta}_{s}+\hat{\theta}_{s}\otimes\hat{\theta}_{r}\big{]}.$

###  C.2 Proof of Theorem 8

By Lemma 5, we have
$\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|\leq\|L_{R}(\mathscr{E})\|+\|S_{R}(\mathscr{E})\|\leq\|L_{R}(\mathscr{E})\|+K_{R}\|\mathscr{E}\|^{2}/g_{R}^{2}$.
It follows that

$\mathbb{P}\big{(}\|\widetilde{\mathscr{P}}_{R}-\widehat{\mathscr{P}}_{R}\|\geq\epsilon+\hat{\sigma}_{1}^{4}K_{R}\epsilon^{2}/g_{R}^{2}\big{)}\leq\mathbb{P}\big{(}\|L_{R}(\mathscr{E})\|\geq\epsilon\big{)}+\mathbb{P}\big{(}\|\mathscr{E}\|\geq\hat{\sigma}_{1}^{2}\epsilon\big{)}.$

The two probabilities on the right hand side are controlled by the next two
lemmas. The conclusion of Theorem 8 follows by combing the next two lemmas and
noting that $\Delta_{0}\leq R+\Delta_{R}$.

######  Lemma 17.

Under the conditions of Theorem 8, it holds that

$\displaystyle\mathbb{P}\big{(}\|L_{R}(\mathscr{E})\|\geq\epsilon\big{)}\leq
8R\times\exp\Big{(}-\frac{C\epsilon^{2}/2}{\frac{\beta
G_{R}^{2}}{2}\big{(}R+\Delta_{R}\big{)}^{2}+[G_{R}\beta\big{(}R+\Delta_{R}\big{)}/\sqrt{2}]\times\epsilon/3}\Big{)}$

for $\epsilon\cdot C\geq\sqrt{\beta
C/2}G_{R}(R+\Delta_{R})+[G_{R}\beta\big{(}R+\Delta_{R}\big{)}/\sqrt{2}]/3$.

######  Lemma 18.

Denote $\Delta_{0}=\mathrm{intdim}(\widehat{\mathscr{C}}_{XX})$ as the
intrinsic dimension of $\widehat{\mathscr{C}}_{XX}$. Under the conditions of
Theorem 8, we have

$\mathbb{P}\left(\left\|\widetilde{\mathscr{C}}_{XX}-\widehat{\mathscr{C}}_{XX}\right\|\geq\epsilon\hat{\sigma}_{1}^{2}\right)\leq
4\Delta_{0}\exp\Big{(}-\frac{C\epsilon^{2}/2}{\beta(R+\Delta_{R})+[\beta(R+\Delta_{R})+1]\epsilon/3}\Big{)},$

for $\epsilon$ satisfying $\epsilon\cdot
C>\sqrt{C\beta(R+\Delta_{R})}+[\beta(R+\Delta_{R})+1]/3$

###  C.3 Proof of Lemma 17

Recall from Lemma 6, we have the relation
$L_{R}(\mathscr{E})=\frac{1}{C}\sum_{c=1}^{C}\frac{\mathscr{Z}_{c}}{N\tilde{p}_{c}}$
with

$\frac{\mathscr{Z}_{c}}{N\tilde{p}_{c}}=\sum_{r=1}^{R}\sum_{s=R+1}^{\infty}\frac{\hat{\sigma}_{r}\hat{\sigma}_{s}}{\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}}\times\frac{\tilde{\xi}_{cr}\tilde{\xi}_{cs}}{\tilde{p}_{c}N}\times\big{[}\hat{\theta}_{r}\otimes\hat{\theta}_{s}+\hat{\theta}_{s}\otimes\hat{\theta}_{r}\big{]}.$

We will apply Lemma 16 to the summation
$\mathscr{S}=\sum_{c=1}^{C}\frac{\mathscr{Z}_{c}}{N\tilde{p}_{c}}$. We first
bound the operator norm of each summand by

$\displaystyle\|\frac{\mathscr{Z}_{c}}{N\tilde{p}_{c}}\|^{2}$
$\displaystyle\leq\|\frac{\mathscr{Z}_{c}}{N\tilde{p}_{c}}\|_{HS}^{2}=2\sum_{r=1}^{R}\sum_{s=R+1}^{\infty}\frac{\hat{\sigma}_{r}^{2}\hat{\sigma}_{s}^{2}}{[\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}]^{2}}\times\frac{\tilde{\xi}_{cr}^{2}\tilde{\xi}_{cs}^{2}}{\tilde{p}_{c}^{2}N^{2}}$
$\displaystyle=2\sum_{r=1}^{R}\sum_{s=R+1}^{\infty}\frac{\hat{\sigma}_{r}^{2}\hat{\sigma}_{R}^{2}}{[\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}]^{2}}\times\frac{\tilde{\xi}_{cr}^{2}(\hat{\sigma}_{s}^{2}\tilde{\xi}_{cs}^{2}/\hat{\sigma}_{R}^{2})}{\tilde{p}_{c}^{2}N^{2}}$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}\frac{2G_{R}^{2}}{\tilde{p}_{c}^{2}N^{2}}\Big{(}\sum_{r=1}^{R}\tilde{\xi}_{cr}^{2}\Big{)}\Big{(}\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}\tilde{\xi}_{cs}^{2}/\hat{\sigma}_{R}^{2}\Big{)}$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\leq}}\frac{G_{R}^{2}}{2\tilde{p}_{c}^{2}N^{2}}\Big{(}\sum_{r=1}^{R}\tilde{\xi}_{cr}^{2}+\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}\tilde{\xi}_{cs}^{2}/\hat{\sigma}_{R}^{2}\Big{)}^{2}.$

In the above, (i) used that
$\frac{\hat{\sigma}_{r}^{2}}{\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}}\leq
G_{R}$ and
$\hat{\sigma}_{r}^{2}-\hat{\sigma}_{s}^{2}\geq\hat{\sigma}_{R}^{2}-\hat{\sigma}_{R+1}^{2}$
for any $r\leq R$ and $s>R$. Inequality (ii) used that
$ab\leq(\frac{a+b}{2})^{2}$ for real numbers $a$ and $b$. For the sampling
probability $\\{p_{n}\\}_{n=1}^{N}$ satisfying $p_{n}\geq
p_{n}^{\text{Exact}}/\beta$, we have

$p_{n}\geq\frac{1}{\beta}\Big{(}\sum_{r=1}^{R}\hat{\xi}_{nr}^{2}+\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}\hat{\xi}_{ns}^{2}/\hat{\sigma}_{R}^{2}\Big{)}/\sum_{m=1}^{N}\Big{(}\sum_{r=1}^{R}\hat{\xi}_{mr}^{2}+\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}\hat{\xi}_{ms}^{2}/\hat{\sigma}_{R}^{2}\Big{)},$

as $\sum_{n=1}^{N}\hat{\xi}_{nr}^{2}/N=1$, we have

$\displaystyle p_{n}$ $\displaystyle\geq\frac{1}{\beta\cdot
N}\Big{(}\sum_{r=1}^{R}\hat{\xi}_{nr}^{2}+\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}\hat{\xi}_{ns}^{2}/\hat{\sigma}_{R}^{2}\Big{)}/\big{(}R+\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}/\hat{\sigma}_{R}^{2}\big{)}$
$\displaystyle\geq\frac{1}{\beta\cdot
N\big{(}R+\Delta_{R}\big{)}}\Big{(}\sum_{r=1}^{R}\hat{\xi}_{nr}^{2}+\sum_{s=R+1}^{\infty}\hat{\sigma}_{s}^{2}\hat{\xi}_{ns}^{2}/\hat{\sigma}_{R}^{2}\Big{)}.$
(C.1)

As a result,
$\|\frac{\mathscr{Z}_{c}}{N\tilde{p}_{c}}\|^{2}\leq\frac{G_{R}^{2}\beta^{2}}{2}\big{(}R+\Delta_{R}\big{)}^{2}$,
which is

$\|\frac{\mathscr{Z}_{c}}{N\tilde{p}_{c}}\|\leq
L=G_{R}\beta\big{(}R+\Delta_{R}\big{)}/\sqrt{2}.$ (C.2)

Next, we bound the variance of the summation
$\mathscr{S}=\sum_{c=1}^{C}\frac{\mathscr{Z}_{c}}{N\tilde{p}_{c}}$. Note that