Sparsity learning via structured functional factor augmentation

Hanteng Ma
School of Statistics and Data Science, Shanghai University of Finance and Economics
Ziliang Shen
School of Statistics and Data Science, Shanghai University of Finance and Economics
Xingdong Feng
School of Statistics and Data Science, Shanghai University of Finance and Economics
Xin Liu
School of Statistics and Data Science, Shanghai University of Finance and Economics The authors gratefully acknowledge please remember to list all relevant funding sources in the unblinded version

Abstract

As one of the most powerful tools for examining the association between functional covariates and a response, the functional regression model has been widely adopted in various interdisciplinary studies. Usually, a limited number of functional covariates are assumed in a functional linear regression model. Nevertheless, correlations may exist between functional covariates in high-dimensional functional linear regression models, which brings significant statistical challenges to statistical inference and functional variable selection. In this article, a novel functional factor augmentation structure (fFAS) is proposed for multivariate functional series, and a multivariate functional factor augmentation selection model (fFASM) is further proposed to deal with issues arising from variable selection of correlated functional covariates. Theoretical justifications for the proposed fFAS are provided, and statistical inference results of the proposed fFASM are established. Numerical investigations support the superb performance of the novel fFASM model in terms of estimation accuracy and selection consistency.

Keywords: correlated functional covariates, functional factor augmentation structure, functional variable selection, factor augmentation regression.

1 Introduction

Functional data, usually referred to as a sequential collection of instances over time with serial dependence, have been widely observed in various scenarios from different scientific disciplines, such as earth, medical, and social sciences (Peng et al., 2005; Centofanti et al., 2021). Instead of employing conventional time series modeling techniques, an underlying trajectory is usually assumed for a functional process and sampled at a set of time points in some interval, so that statistical modeling and inference are feasible, including estimation of the mean process and its intrinsic covariance structure, trajectory recovery, and prediction of a functional; see, for example, Ramsay & Silverman (2005); Yao et al. (2005); Ferraty (2006); Hall & Horowitz (2007); Hörmann & Kokoszka (2010); Cuevas (2014); Hörmann et al. (2015); Aneiros, Horová, Hušková & Vieu (2022); Petersen (2024), among others.

Conventionally, functional data analysis has focused on a single or fixed number of separate functional series, while recent studies have explored associations between multiple functional series, and the most widely adopted model may be the functional linear model, which assumes to (linearly) link a response with a bunch of covariates and at least one functional covariate by using unknown functional coefficient curves. Such a model offers a statistical tool to infer linear association between the response and functional covariates. Several studies have attempted along this path, employing parametric, non-parametric, and semi-parametric models; see, for instance, Cardot et al. (2003); Yao et al. (2005); Li & Hsing (2007); Müller & Yao (2008); Yuan & Cai (2010); Chen et al. (2011); Kong et al. (2016); Chen et al. (2022), among others.

However, two main issues arise in multivariate functional data analysis in practice. One is that functional series may be correlated with each other, especially in high dimensions or less dense observation cases, while the current functional linear model analysis may typically assume independence between multiple functional covariates. To address this problem, some efforts are paid, for example, to assume pairwise covariance structures for functional covariates explicitly (Chiou & Müller, 2016). However, these methods turn out to be less robust or accurate due to lack of flexibility and computational efficiency. Alternatively, one can describe such associations by using a factor model with lower ranks (Bai, 2003; Fan et al., 2020), where all functional covariates are assumed to share a finite number of certain latent functional factors (Castellanos et al., 2015), without assuming any explicit correlation structure for multivariate functional covariates. A few studies have touched the so-called functional factor models, e.g., Hays et al. (2012); Chen et al. (2021); Yang & Ling (2024), while these methods provide limited exploration of common functional factors shared by functional covariates. To the best of our knowledge, very limited research has focused on how to capture common associations for multivariate functional series efficiently and effectively with inferential justifications from a statistical perspective.

Another issue is how to select useful functional covariates in multivariate functional linear regression models when correlations between functional covariates exist, which is also frequently needed in practice, such as Internet of Things (IoT) data analysis (Gonzalez-Vidal et al., 2019), antibiotic resistance prediction (Jiménez et al., 2020), nonsmall cell lung cancer (NSCLC) clinical study (Fang et al., 2015), and stock market forecasting (Nti et al., 2019; Htun et al., 2023). A common way to achieve functional variable selection is to impose a penalty on the corresponding functional coefficient curves in a group-wise manner, using popular penalties such as the lasso (Tibshirani, 1996), SCAD (Fan & Li, 2001) and MCP (Zhang, 2010). A few studies have touched on functional variable selection in certain scenarios; see, for example, Matsui & Konishi (2011); Kong et al. (2016); Aneiros, Horová, Hušková & Vieu (2022), among others. However, when correlation exists between functional covariates, such a strategy turns out to be less accurate in functional coefficient estimation and fails to capture the truly useful ones, as demonstrated in simulation studies in this article. Consequently, it remains statistically challenging how to select useful functional covariates with selection consistency in multivariate functional linear models with correlated functional covariates in high dimensions.

Inspired by the challenges above, we firstly propose a novel functional factor augmentation structure (fFAS) to capture associations for correlated multivariate functional data, and further propose a functional factor augmentation selection model (fFASM) to achieve selection of functional covariates in high dimensions when correlations exist between functional covariates using the penalized method. Not only is the correlation addressed without assuming an explicit covariance structure, but theoretical properties of the estimated functional factors are also established. Further, pertaining to the correlated functional covariates, the proposed fFASM method successfully captures the useful functional covariates simultaneously in the context of functional factor models. Numerical studies on both simulated and real datasets support the superior performance of the proposed fFASM. The main contributions of the proposed method are two-fold as follows.

•

We propose a feasible time-independent functional factor augmentation structure (fFAS) for functional data, and establish theoretical justifications for statistical inference, revealing valuable insights into current literature. Also, the impact from truncation inherent in functional data analysis to the proposed fFAS is discussed, and solutions are provided to enhance robustness and applicability of our model. A key result is how the difference may be addressed and controlled by comparing the true factor model and the estimated one when using truncated expansion.
•

A multivariate functional linear regression model with correlated functional covariates is proposed, by employing the proposed functional structural factor. As correlations often exist in multiple functional covariates in high dimensions but may be difficult to estimate, the proposed fFASM method decomposes the functional covariates into two parts with low correlations with each other, and simultaneously expands the dimensions of parameters to be estimated. By this way our approach improves the performance of functional variable selection.

The rest of the paper is organized as follows. Section 2 introduces a novel fFAS for functional processes with its statistical properties, and the fFASM is further proposed in Section 3 for correlated functional covariates and functional variable selection in detail with theoretical justifications. Section 4 employs simulated data to examine the proposed method in various scenarios, and Section 5 presents its applications on two sets of real-world datasets. Section 6 concludes the article with discussions.

Notation: $\bm{I}_{n}$ denotes the $n\times n$ identity matrix; $\bm{0}$ refers to the $n\times m$ zero matrix; $\bm{0}_{n}$ and $\bm{1}_{n}$ represent the all-zero and all-one vectors in $\mathbb{R}^{n}$ , respectively. For a vector $\bm{a}$ , $\|\bm{a}\|_{2}$ denotes its Euclidean norm. For a vector $\bm{v}\in\mathbb{R}^{p}$ and $S\subseteq[p]$ , denote $\bm{v}_{S}=\left(\bm{v}_{i}\right)_{i\in S}$ as its sub-vector. For a matrix $\bm{M}\in\mathbb{R}^{n\times m},I\subseteq[n]$ and $J\subseteq[m]$ , define $\bm{M}_{IJ}=\left(\bm{M}_{ij}\right)_{i\in I,j\in J},\bm{M}_{I.}=\left(\bm{M}_{ij}\right)_{i\in I,j\in[m]}$ and $\bm{M}_{.J}=\left(\bm{M}_{ij}\right)_{i\in[n],j]}$ , and its matrix entry-wise max norm is denoted as $\|\bm{M}\|_{\mathop{\rm max}}=\mathop{\rm max}_{i,j}\left|M_{ij}\right|$ , and $\|\bm{M}\|_{F}$ and $\|\bm{M}\|_{p}$ as its Frobenius and induced $p$ -norms, respectively. Denote $\lambda_{\min}(\bm{M})$ as the minimum eigenvalue of $\bm{M}$ if it is symmetric. Let $\nabla$ and $\nabla^{2}$ be the gradient and Hessian operators. For $f:\mathbb{R}^{p}\rightarrow\mathbb{R}$ and $I,J\subseteq[p]$ , define $\nabla_{I}f(x)=(\nabla f(x))_{I}$ and $\nabla_{IJ}^{2}f(x)=\left(\nabla^{2}f(x)\right)_{IJ}$ . $N(\boldsymbol{\mu},\boldsymbol{\Sigma})$ refers to the normal distribution with mean $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$ .

2 A novel functional factor augmentation structure

2.1 Functional observations and its expansion

To start with, suppose that $G$ functional processes are observed and collected sequentially in a sample $\{W_{i}^{(g)}(t_{ij}^{(g)}),t_{ij}^{(g)}\in\mathcal{T};~{}g=1,\ldots,G,~{}i=1,\ldots,n,~{}j=1,\ldots,n_{i}^{(g)}\}$ for the $g$ -th functional process from the $i$ -th subject at time $t_{ij}$ in a certain time interval $\mathcal{T}$ , and we use $W_{i}^{(g)}(t_{j})$ for abbreviation without confusion. Usually in the analysis of a single functional process, $W_{i}^{(g)}(t_{j})$ is assumed to consist of the underlying process $X_{i}^{(g)}(t_{j})$ and an independent noise $\varepsilon_{i}^{(g)}(t_{j})$ , i.e., $W_{i}^{(g)}(t_{j})=X_{i}^{(g)}(t_{j})+\varepsilon_{i}^{(g)}(t_{j})$ , where $X_{i}^{(g)}(\cdot)$ and $\varepsilon_{i}^{(g)}(\cdot)$ are assumed to be identically and independently distributed (i.i.d.) with a mean function and a covariance structure, respectively, and $\mathbb{E}(\varepsilon_{i}^{(g)}(t))=0,\mathbb{E}(X_{i}^{(g)}(t))=\mu(t)$ . For convenience, we centralize these functional processes and still use the notation $X_{i}^{(g)}(\cdot)$ to represent them. To recover the functional trajectory $X_{i}^{(g)}(t)$ , it is usually smoothed by assuming an expansion over a set of pre-specified orthogonal basis functions $\{\phi_{0j}^{(g)}(\cdot),j=1,\ldots,m_{g}\}$ as

X_{i}^{(g)}(t)=\sum_{j=1}^{m_{g}}a_{0,ij}^{(g)}\phi_{0j}^{(g)}(t)+e_{0i}^{(g)}(t)=\bm{a}_{0,i}^{(g)}{}^{\top}\bm{\phi}_{0}^{(g)}(t)+e_{0i}^{(g)}(t),

(1)

where $\bm{a}_{0,i}^{(g)}=(a_{0,i1}^{(g)},\ldots,a_{0,im_{g}}^{(g)})^{\top}$ are the time-invariant coefficients, $\bm{\phi}_{0}^{(g)}(t)=(\phi_{01}^{(g)}(t),\ldots,\phi_{0m_{g}}^{(g)}(t))^{\top}$ , and $e_{0i}^{(g)}(\cdot)$ is the residual orthogonal to $\phi_{0j}^{(g)}(\cdot)$ . As in practice, $\bm{\phi}_{0}^{(g)}(\cdot)$ is usually unknown, an identifiable estimator $\widehat{\bm{a}}_{0,i}^{(g)}$ can be obtained under some conditions, shown in Lemma 1.

Lemma 1

For $X_{i}^{(g)}(\cdot)$ which has the following properties,

•

$\mathbb{E}(X_{i}^{(g)}(t))=\mathbb{E}(a_{0,ij}^{(g)})=\mathbb{E}(a_{0,ij}^{(g)}a_{0,ik}^{(g)})=\mathbb{E}(e_{0,i}^{(g)}(t))=\mathbb{E}(a_{0,ij}^{(g)}e_{0i}^{(g)}(t))=0$ , where $k\neq j$ ,
•

$\mathbb{E}(X_{i}^{(g_{1})}(s)e_{0i}^{(g_{2})}(t))=0$ for $g_{1}\neq g_{2}$ , $g_{1},g_{2}=1,\ldots,G$ , $s,t\in\mathcal{T}$ ,
•

$\operatorname{Var}(a_{0,i1}^{(g)})>\operatorname{Var}(a_{0,i2}^{(g)})>\ldots>\operatorname{Var}(a_{0,im_{g}}^{(g)})$ ,
•

all the eigenvalues of $\operatorname{Cov}(e_{0,i}^{(g)}(s),e_{0,i}^{(g)}(t))$ are less than $\lambda_{\min}(\operatorname{Cov}(\bm{a}_{0,i}^{(g)}))$ , where $\lambda_{\min}(\operatorname{Cov}(\bm{a}_{0,i}^{(g)}))$ is the smallest eigenvalue of $\operatorname{Cov}(\bm{a}_{0,i}^{(g)})$ ,

where the $j$ -th eigenfunction $\gamma_{j}^{(g)}(\cdot)$ of $\operatorname{Cov}(X_{i}^{(g)}(s),X_{i}^{(g)}(t))$ is $\phi_{0j}^{(g)}(\cdot)$ , so that the $j$ -th functional score $a_{ij}^{(g)}\coloneqq\int_{\mathcal{T}}X_{i}^{(g)}(t)\gamma_{j}^{(g)}(t)dt=a_{0,ij}^{(g)}$ .

Remark 1

Lemma 1 states $\bm{a}_{i}^{(g)}=({a}_{i1}^{(g)},\ldots,{a}_{i{m_{g}}}^{(g)})^{\top}$ is a reasonable approximation of $\bm{a}_{0,i}^{(g)}$ . All eigenvalues of $\operatorname{Cov}(e_{0,i}^{(g)}(s),e_{0,i}^{(g)}(t))$ are required to be less than $\lambda_{\min}(\operatorname{Cov}(\bm{a}_{0,i}^{(g)}))$ , indicating that segregation can be conducted based on contribution to the variance of $X_{i}^{(g)}(t)$ .

The whole recovery process can be achieved by the popular functional principal component analysis (fPCA) with the Karhunen-Loéve (KL) expansion (Yao et al., 2005).

2.2 A functional factor augmentation structure

As correlations may exist between multivariate functional processes, we propose a functional factor augmentation structure (fFAS) to address the issue. Consider a simplified scenario with only two correlated functional processes, $X_{i}^{(1)}(t)$ and $X_{i}^{(2)}(t)$ , generated by

	$\displaystyle X_{i}^{(1)}(t)$	$\displaystyle=\bm{a}_{0,i}^{(1)}{}^{\top}\cdot\bm{\phi}_{0}^{(1)}(t)+e_{0i}^{(1)}(t),$		(2)
	$\displaystyle X_{i}^{(2)}(t)$	$\displaystyle=\bm{a}_{0,i}^{(2)}{}^{\top}\cdot\bm{\phi}_{0}^{(2)}(t)+e_{0i}^{(2)}(t),$

where $e_{0i}^{(1)}(t)$ and $e_{0i}^{(2)}(t)$ are independent of each other. Assume $\bm{a}_{0,i}^{(g_{1})}$ and $e_{0i}^{(g_{2})}(t)$ are uncorrelated for $g_{1},g_{2}=1,\ldots,G$ , where $G=2$ , and hence the correlation between $X_{i}^{(1)}(t)$ and $X_{i}^{(2)}(t)$ arises only from that between $\bm{a}_{0,i}^{(1)}$ and $\bm{a}_{0,i}^{(2)}$ . To capture such a correlation, it is assumed that each $\bm{a}_{0,i}^{(g)}$ shares $K$ common underlying factors using a linear combination

\displaystyle\left(\begin{array}[]{ll}\bm{a}_{0,i}^{(1)}\\ \bm{a}_{0,i}^{(2)}\end{array}\right)=\left(\begin{array}[]{ll}B^{(1)}\\ B^{(2)}\end{array}\right)\bm{f}_{i}+\left(\begin{array}[]{ll}\bm{u}_{i}^{(1)}\\ \bm{u}_{i}^{(2)}\end{array}\right)=\bm{B}\bm{f}_{i}+\bm{u}_{i},

(9)

where $\bm{B}=\left(B^{(1)}{}^{\top},B^{(2)}{}^{\top}\right)^{\top}$ is a $p\times K$ factor loading matrix with $p=\sum_{g=1}^{G}m_{g}$ , $\bm{f}_{i}$ is a $K\times 1$ vectorized latent factors, and $\bm{u}_{i}$ is an idiosyncratic component independent of $\bm{f}_{i}$ which carries a weak correlation. Note that the covariance $C_{\bm{a}_{0}}:=\operatorname{Cov}(\bm{a}_{0,i})=\operatorname{Cov}(\bm{B}\bm{f}_{i}+\bm{u}_{i})=\bm{B}\cdot\operatorname{Cov}(\bm{f}_{i})\cdot\bm{B}^{\top}+\operatorname{Cov}(\bm{u}_{i})$ . For model identifiability, it is assumed that $\operatorname{Cov}(\bm{f}_{i})=\bm{I}_{K}$ and $\operatorname{Cov}(\bm{u}_{i})=\omega\Lambda_{u}$ with a $p\times p$ matrix $\Lambda_{u}$ , where $\|\Lambda_{u}\|_{\mathop{\rm max}}\leq C_{u}$ , a constant depending on the distribution of $\bm{u}_{i}$ . By the spectral decomposition, $\bm{B}=\left(\sqrt{\lambda_{1}^{(B)}}\bm{\xi}_{1}^{(B)},\ldots,\sqrt{\lambda_{K}^{(B)}}\bm{\xi}_{K}^{(B)}\right)$ with the eigenvalues $\lambda_{1}^{(B)}>\ldots>\lambda_{K}^{(B)}$ of $\bm{BB}^{\top}$ and the corresponding eigenvectors $\{\bm{\xi}_{k}^{(B)},k=1,\ldots,K\}$ . In this way, $X_{i}^{(1)}(t)$ and $X_{i}^{(2)}(t)$ in (2) can be expanded as

	$\displaystyle X_{i}^{(1)}(t)$	$\displaystyle={\bm{\phi}_{0}^{(1)}(t)}{}^{\top}(B^{(1)}\bm{f}_{i}+\bm{u}_{i}^{(1)})+e_{0i}^{(1)}(t)={{\tilde{\bm{\phi}_{0}}}{}^{(1)}(t)}{}^{\top}\bm{f}_{i}+{\bm{\phi}_{0}^{(1)}(t)}{}^{\top}\bm{u}_{i}^{(1)}+e_{0i}^{(1)}(t),$		(10)
	$\displaystyle X_{i}^{(2)}(t)$	$\displaystyle={\bm{\phi}_{0}^{(2)}(t)}{}^{\top}(B^{(2)}\bm{f}_{i}+\bm{u}_{i}^{(2)})+e_{0i}^{(2)}(t)={{\tilde{\bm{\phi}_{0}}}{}^{(2)}(t)}{}^{\top}\bm{f}_{i}+{\bm{\phi}_{0}^{(2)}(t)}{}^{\top}\bm{u}_{i}^{(2)}+e_{0i}^{(2)}(t),$

where ${{\tilde{\bm{\phi}_{0}}}{}^{(g)}(t)}{}^{\top}={\bm{\phi}_{0}^{(g)}(t)}{}^{\top}B^{(g)}$ . This indicates $X_{i}^{(g)}(t)$ can be decomposed into two correlated parts, namely, the functional factor part $\boldsymbol{f}_{i}$ and the weakly correlated part ${\bm{\phi}_{0}^{(g)}(t)}{}^{\top}\bm{u}_{i}^{(g)}$ , plus an independent error term. We provide an example to illustrate such a fFAS when two functional covariates share a linear structure.

Example 1

Suppose $X_{i}^{(1)}{(\cdot)}$ and $X_{i}^{(2)}{(\cdot)}$ are associated with a linear structure as

\displaystyle X_{i}^{(2)}{(t)}=\mathbb{E}(X_{i}^{(2)}(t)|X_{i}^{(1)}{(t)})+\epsilon_{i}(t)=\int_{\mathcal{T}}\beta(s,t)X_{i}^{(1)}(s)ds+\epsilon_{i}{(t)},

with $\beta(s,t)=\sum_{k=1}^{\infty}\sum_{m=1}^{\infty}\frac{\mathbb{E}(a_{0,im}^{(1)}a_{0,ik}^{(2)})}{\mathbb{E}(a_{0,im}^{(1)}{}^{2})}\phi_{0k}^{(2)}(t)\phi_{0m}^{(1)}(s)$ . Then it is obtained that

	$\displaystyle\int_{\mathcal{T}}\beta(s,t)X_{i}^{(1)}(s)ds$	$\displaystyle=\int_{\mathcal{T}}\beta(s,t)\sum_{j=1}^{\infty}a_{0,ij}^{(1)}\phi_{0j}^{(1)}(s)ds$
		$\displaystyle=\int_{\mathcal{T}}\sum_{k=1}^{\infty}\sum_{m=1}^{\infty}\frac{\mathbb{E}(a_{0,im}^{(1)}a_{0,ik}^{(2)})}{\mathbb{E}(a_{0,im}^{(1)}{}^{2})}\phi_{0k}^{(2)}(t)\phi_{0m}^{(1)}(s)\sum_{j=1}^{\infty}a_{0,ij}^{(1)}\phi_{0j}^{(1)}(s)ds$
		$\displaystyle=\sum_{k=1}^{\infty}\left(\sum_{j=1}^{\infty}\frac{\mathbb{E}(a_{0,ij}^{(1)}a_{0,ik}^{(2)})}{\mathbb{E}(a_{0,ij}^{(1)}{}^{2})}a_{0,ij}^{(1)}\right)\phi_{0k}^{(2)}(t),$

so that $a_{0,ik}^{(2)}\approx\sum_{j=1}^{\infty}\frac{\mathbb{E}(a_{0,ij}^{(1)}a_{0,ik}^{(2)})}{\mathbb{E}(a_{0,ij}^{(1)}{}^{2})}a_{0,ij}^{(1)}=\sum_{j=1}^{m_{1}}\frac{\mathbb{E}(a_{0,ij}^{(1)}a_{0,ik}^{(2)})}{\mathbb{E}(a_{0,ij}^{(1)}{}^{2})}a_{0,ij}^{(1)}+\sum_{j=m_{1}+1}^{\infty}\frac{\mathbb{E}(a_{0,ij}^{(1)}a_{0,ik}^{(2)})}{\mathbb{E}(a_{0,ij}^{(1)}{}^{2})}a_{0,ij}^{(1)}$ . With a truncation of the first $m_{1}$ and $m_{2}$ scores in $X_{i}^{(1)}(t)$ and $X_{i}^{(2)}(t)$ , and combining them as ${\bm{a}_{0,i}}=({\bm{a}_{0,i}^{(1)}}^{\top},{\bm{a}_{0,i}^{(2)}}^{\top})^{\top}$ , and a matrix $\bm{E}$ with elements $E_{kj}=\frac{\mathbb{E}(a_{0,ij}^{(1)}a_{0,ik}^{(2)})}{\mathbb{E}(a_{0,ij}^{(1)}{}^{2})}$ and a diagonal matrix $\Lambda$ with elements $\Lambda_{jj}=\mathbb{E}(a_{0,ij}^{(1)}{}^{2})$ , it is easily obtained that

\left(\begin{array}[]{ll}\bm{a}_{0,i}^{(1)}\\ \bm{a}_{0,i}^{(2)}\end{array}\right)\approx\left(\begin{array}[]{cc}\bm{I}_{m_{1}}\\ \bm{E}\end{array}\right)\bm{a}_{0,i}^{(1)}=\left(\left(\begin{array}[]{cc}\bm{I}_{m_{1}}\\ \bm{E}\end{array}\right)\Lambda^{\frac{1}{2}}\right)\left(\Lambda^{-\frac{1}{2}}\bm{a}_{0,i}^{(1)}\right)=\left(\left(\begin{array}[]{cc}\bm{I}_{m_{1}}\\ \bm{E}\end{array}\right)\Lambda^{\frac{1}{2}}\bm{P}\right)(\bm{P}^{\top}\Lambda^{-\frac{1}{2}}\bm{a}_{0,i}^{(1)}),

where $\bm{P}$ is an orthogonal matrix, so that $\left(\left(\begin{array}[]{cc}\bm{I}_{m_{1}}\\ \bm{E}\end{array}\right)\Lambda^{\frac{1}{2}}\bm{P}\right)^{\top}\left(\left(\begin{array}[]{cc}\bm{I}_{m_{1}}\\ \bm{E}\end{array}\right)\Lambda^{\frac{1}{2}}\bm{P}\right)$ is diagonal, and $\bm{f}_{i}=\bm{P}^{\top}\Lambda^{-\frac{1}{2}}\bm{a}_{0,i}^{(1)}$ with $\operatorname{Cov}(\bm{f}_{i})=\bm{I}_{m_{1}}$ . Note that no matter what values of $m_{1}$ and $m_{2}$ are in practice, one can always obtain such a functional factor augmentation structure.

Next, we consider a more general case where the structure $\bm{a}_{0,i}^{(g)}$ contains correlations for $g=1,\ldots,G$ , where $\mathbb{E}(a_{0,ij}^{(g)}a_{0,ik}^{(g)})\neq 0~{}(\text{for all }0<j<k<m_{g})$ . Then, Lemma 2 shows the relation between $\bm{a}_{i}^{(g)}$ and $\bm{a}_{0,i}^{(g)}$ in this circumstance that $\bm{a}_{i}^{(g)}$ is still an approximation of $\bm{a}_{0,i}^{(g)}$ by imposing an orthogonal rotation induced by the basis functions in K-L expansion.

Lemma 2

Under the conditions in Lemma 1 without $\mathbb{E}(a_{0,ij}^{(g)}a_{0,ik}^{(g)})=0$ , there exists an orthogonal matrix $P^{(g)}$ , such that $\bm{\phi}_{0}^{(g)}(\cdot)=P^{(g)}({\gamma}_{1}^{(g)}(\cdot),\ldots,{\gamma}_{m_{g}}^{(g)}(\cdot))^{\top}$ , and

\displaystyle\left(\begin{array}[]{ll}a_{0,i1}^{(g)}\\ a_{0,i2}^{(g)}\\ \vdots\\ a_{0,im_{g}}^{(g)}\end{array}\right)=P^{(g)}\left(\begin{array}[]{ll}a_{i1}^{(g)}\\ a_{i2}^{(g)}\\ \vdots\\ a_{im_{g}}^{(g)}\end{array}\right).

Furthermore, denote $\bm{a}_{i}^{\top}=(\bm{a}_{i}^{(1)}{}^{\top},\ldots,\bm{a}_{i}^{(G)}{}^{\top})$ and $\bm{a}_{0,i}^{\top}=(\bm{a}_{0,i}^{(1)}{}^{\top},\ldots,\bm{a}_{0,i}^{(G)}{}^{\top})$ . Then

\displaystyle\operatorname{Cov}(\bm{a}_{i})=\begin{pmatrix}P^{(1)}&\bm{0}&\cdots&\bm{0}\\ \bm{0}&P^{(2)}&\cdots&\bm{0}\\ \vdots&\vdots&\ddots&\vdots\\ \bm{0}&\bm{0}&\cdots&P^{(G)}\end{pmatrix}^{\top}\operatorname{Cov}(\bm{a}_{0,i})\begin{pmatrix}P^{(1)}&\bm{0}&\cdots&\bm{0}\\ \bm{0}&P^{(2)}&\cdots&\bm{0}\\ \vdots&\vdots&\ddots&\vdots\\ \bm{0}&\bm{0}&\cdots&P^{(G)}\end{pmatrix},

where

\displaystyle{\operatorname{Cov}(\bm{a}_{i})}=\begin{pmatrix}\Sigma_{m_{1}}&\Sigma_{m_{1}m_{2}}&\cdots&\Sigma_{m_{1}m_{G}}\\ \Sigma_{m_{2}m_{1}}&\Sigma_{m_{2}}&\cdots&\Sigma_{m_{2}m_{G}}\\ \vdots&\vdots&\ddots&\vdots\\ \Sigma_{m_{G}m_{1}}&\Sigma_{m_{G}m_{2}}&\cdots&\Sigma_{m_{G}}\end{pmatrix},

and

\displaystyle{\operatorname{Cov}(\bm{a}_{0,i})}=\begin{pmatrix}\Sigma_{0,m_{1}}&\Sigma_{0,m_{1}m_{2}}&\cdots&\Sigma_{0,m_{1}m_{G}}\\ \Sigma_{0,m_{2}m_{1}}&\Sigma_{0,m_{2}}&\cdots&\Sigma_{0,m_{2}m_{G}}\\ \vdots&\vdots&\ddots&\vdots\\ \Sigma_{0,m_{G}m_{1}}&\Sigma_{0,m_{G}m_{2}}&\cdots&\Sigma_{0,m_{G}}\end{pmatrix}.

Remark 2

Lemma 2 indicates that if there is a fFAS on $\bm{a}_{0,i}$ , then there will be a fFAS on $\bm{a}_{i}$ by the fact that

\displaystyle\bm{a}_{i}=\begin{pmatrix}P^{(1)}&\bm{0}&\cdots&\bm{0}\\ \bm{0}&P^{(2)}&\cdots&\bm{0}\\ \vdots&\vdots&\ddots&\vdots\\ \bm{0}&\bm{0}&\cdots&P^{(G)}\end{pmatrix}^{\top}\bm{a}_{0,i}=\bm{P}^{\top}\bm{a}_{0,i}=\bm{P}^{\top}\bm{B}\bm{f}_{i}+\bm{P}^{\top}\bm{u}_{i},

(11)

and treating $\bm{P}^{\top}\bm{B}$ and $\bm{P}^{\top}\bm{u}_{i}$ as the updated loading matrix and the idiosyncratic component in (9), respectively, where $\bm{P}$ is still an orthogonal matrix. More specifically, if $\bm{a}_{0,i}$ has a fFAS induced by $\bm{f}_{i}$ , $\bm{a}_{i}$ also has such a fFAS with the same factor $\bm{f}_{i}$ . Consequently, even if the functional covariates are correlated with each other, we can still use the functional scores $\bm{a}_{i}$ to estimate $\bm{a}_{0,i}$ as if they were uncorrelated.

To further get the estimates of the loading matrix ${\bm{B}}$ and the functional factors $\bm{F}=\left(\bm{f}_{1},\ldots,\bm{f}_{n}\right)^{\top}\in\mathbb{R}^{n\times K}$ after obtaining the score matrix $\bm{A}=(\bm{a}_{1},\ldots,\bm{a}_{n})^{\top}$ with $\bm{a}_{i}$ in (11), we employ the principal component analysis (Bai, 2003) on the covariance of $\bm{A}$ , and further obtain that $\widehat{\bm{U}}={\bm{A}}-\widehat{\bm{F}}\widehat{\bm{B}}^{\top}$ , where $\widehat{\bm{U}}=\left(\widehat{\bm{u}}_{1},\ldots,\widehat{\bm{u}}_{n}\right)^{\top}\in\mathbb{R}^{n\times p}$ . More specifically, the columns of $\widehat{\bm{F}}/\sqrt{n}$ are the eigenvectors of $\bm{AA}^{\top}$ corresponding to the top ${K}$ eigenvalues, and $\widehat{\bm{B}}=n^{-1}\bm{A}^{\top}\widehat{\bm{F}}$ . This is similar to the case $\widehat{\bm{B}}=\left(\sqrt{\lambda_{1}}\bm{\xi}_{1},\ldots,\sqrt{\lambda_{{K}}}\bm{\xi}_{{K}}\right)$ and $\widehat{\bm{F}}=\bm{A}\widehat{\bm{{B}}}\operatorname{diag}\left(\lambda_{1}^{-1}\ldots,\lambda_{{K}}^{-1}\right)$ , where $\left\{\lambda_{j}\right\}_{j=1}^{{K}}$ and $\left\{\boldsymbol{\xi}_{j}\right\}_{i=1}^{{K}}$ are top ${{K}}$ eigenvalues in descending order and their associated eigenvectors of the sample covariance matrix.

A practical issue is how to determine the number of factors $K$ . Given that latent factors, loadings, and idiosyncratic components are all unobservable, a conditional sparsity criterion is then adopted (Ahn & Horenstein, 2013; Fan et al., 2020). Specifically, let $\lambda_{k}\left({\bm{A}}^{\top}{\bm{A}}\right)$ denote the $k$ -th largest eigenvalue of ${\bm{A}}^{\top}{\bm{A}}$ , $K_{\mathop{\rm max}}$ be a prespecified upper bound for $K$ , and $C_{n}$ be a constant dependent on $n$ and $p=\sum_{g=1}^{G}m_{g}$ , and then $K$ is determined by

\widehat{K}=\underset{k\leq K_{\mathop{\rm max}}}{\operatorname{argmin}}\frac{\lambda_{k+1}\left({\bm{A}}^{\top}{\bm{A}}\right)+C_{n}}{\lambda_{k}\left({\bm{A}}^{\top}{\bm{A}}\right)+C_{n}}

(12)

for a given $C_{n}$ and $K_{\mathop{\rm max}}$ . An alternative criterion may be the information criteria proposed by Bai & Ng (2002) and Fan et al. (2013), referred to as PC and IC, respectively.

2.3 Properties of functional factor augmentation structure

In this subsection, we evaluate the estimation error of the proposed functional factors $\bm{f}_{i}$ . Denote the first $K$ eigenvectors and eigenvalues of $\bm{B}\bm{B}^{\top}$ as $\bm{\xi}^{(B)}_{k}$ , and $\Lambda_{B}=\operatorname{diag}\{\lambda_{1}^{(B)},\ldots,\lambda_{K}^{(B)}\}$ , respectively, where $\lambda_{k}^{(B)}$ , $k=1,\dots,K$ are sorted in a descending order. The largest eigenvalue of $\operatorname{Cov}(\bm{u}_{i})$ is $\lambda_{\mathop{\rm max}}^{\bm{u}}$ , where $\operatorname{Cov}(\bm{u}_{i})=\omega\Lambda_{u}$ , the covariance matrix of $\bm{a}_{i}$ , $\operatorname{Cov}(\bm{a}_{i})$ , equals to $\mathbb{E}(\bm{a}_{i}\bm{a}_{i}^{\top})=\mathbb{E}(\bm{A}^{\top}\bm{A}/n)=\bm{B}\bm{B}^{\top}+\mathbb{E}(\bm{U}^{\top}\bm{U}/n)$ from (11). The covariance can be further expanded using a perturbation

\displaystyle C_{\bm{a}}(\omega)=\bm{B}\bm{B}^{\top}+\Lambda_{\bm{u}}\omega+O(\omega^{2}),

where $\Lambda_{\bm{u}}$ is a perturbing matrix with $\omega$ such that

	$\displaystyle\lambda_{k}^{(C)}(\omega)$	$\displaystyle=\lambda_{k}^{(B)}+(\bm{\xi}^{(B)}_{k}{}^{\top}\Lambda_{\bm{u}}\bm{\xi}^{(B)}_{k})\omega+O(\omega^{2}),$
	$\displaystyle\bm{\xi}_{k}^{(C)}(\omega)$	$\displaystyle=\bm{\xi}^{(B)}_{k}-(\bm{B}\bm{B}^{\top}-\lambda_{k}^{(B)}I)^{+}\Lambda_{\bm{u}}\bm{\xi}_{k}^{(B)}\omega+O(\omega^{2}),$

with $\lambda_{k}^{(C)}(\omega)$ and $\bm{\xi}^{(C)}_{k}$ , $k=1,\dots,K$ being the first $K$ eigenvalues and eigenvectors of $C_{\bm{a}}(\omega)$ (Shi, 1997), respectively. Based on $C_{\bm{a}}(\omega)$ we have ${\bm{B}(\omega)}=\left(\sqrt{\lambda_{1}^{(C)}(\omega)}\bm{\xi}_{1}^{(C)}(\omega),\ldots,\sqrt{\lambda_{K}^{(C)}(\omega)}\bm{\xi}_{K}^{(C)}(\omega)\right)$ and ${\bm{F}}(\omega)=\bm{A}\bm{B}(\omega)\left[\operatorname{diag}\left(\lambda_{1}^{(C)}\ldots,\lambda_{K}^{(C)}\right)\right]^{-1}$ . Accordingly, the following lemma shows the estimation error of the proposed functional structure factors for a given $\bm{B}$ .

Lemma 3

For the fixed truncation numbers $m_{1},\ldots,m_{G}$ , if $\underset{1\leq i<j\leq K}{\min}(\lambda_{i}^{(B)}-\lambda_{j}^{(B)})>c_{1}$ and $\lambda_{\mathop{\rm max}}^{\bm{u}}<c_{2}$ for some constant $c_{1},c_{2}>0$ , then

	$\displaystyle\\|{\bm{f}}_{i}^{\top}(\omega)-\bm{f}_{i}^{\top}\\|_{2}^{2}$	$\displaystyle=O_{P}(\omega^{2}\\|\bm{a}_{i}\\|_{2}^{2}+\\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\\|_{2}^{2}),~{}\quad{\mbox{and}}$
	$\displaystyle\\|{\bm{u}}_{i}^{\top}(\omega)-\bm{u}_{i}^{\top}\\|_{2}^{2}$	$\displaystyle=O_{P}(\omega^{2}\\|\bm{a}_{i}\\|_{2}^{2}+\\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\bm{B}^{\top}\\|_{2}^{2})$

with $\bm{f}_{i}^{\top}(\omega)$ and $\bm{u}_{i}^{\top}(\omega)$ calculated based on $C_{\bm{a}}(\omega)$ .

Remark 3

In terms of errors induced by detecting functional structure factors, there are mainly two sources. One is from the estimation error of eigenvalues and eigenvectors, and the other is from the projection of $\bm{u}_{i}$ in the column space of $\bm{B}$ .

Furthermore, in real applications when only repeated observations $\{(t_{ij}^{(g)},W_{i}^{(g)}(t_{ij}^{(g)})),t_{ij}^{(g)}\in\mathcal{T},j=1,\ldots,n_{i}^{(g)}\}$ are available for each subject, a common practice is to estimate a smooth trajectory $\{\widehat{X}_{i}^{(g)}(\cdot),i=1,\ldots,n,g=1,\ldots,G\}$ and hence obtain $\widehat{\bm{A}}$ , $\widehat{\bm{B}}$ and $\widehat{\bm{F}}$ . To further consider the estimation error for $\widehat{\bm{F}}$ , more assumptions are required to describe the properties of the scores estimated by fPCA. For each functional covariate, denote the covariance function $\text{Cov}_{X}^{(g)}(s,t)=\operatorname{Cov}(X_{i}^{(g)}(s),X_{i}^{(g)}(t)),g=1,\ldots,G$ , and $\tau_{k}^{(g)}$ is the k-th eigenvalue of covariance functions $\text{Cov}_{X}^{(g)}(s,t)$ . Assume for $g=1,\ldots,G$ , we have $m_{g}\leq s_{n}$ .

Assumption 1

(A1)

(The decay rates of eigenvalues) For $g=1,\ldots,G$ , $\tau_{1}^{(g)}<\infty$ , $\tau_{k}^{(g)}-\tau_{k+1}^{(g)}\geq Ck^{-a-1}$ for $k\geq 1$ .
(A2)

(The common truncation parameter $s_{n}$ cannot be too large) $(s_{n}^{2a+2}+s_{n}^{a+4})/n=o(1).$

This implies that $\tau_{k}^{(g)}\geq Ck^{-a}$ . As the covariance functions $\text{Cov}_{X}^{(g)}(s,t),g=1,\ldots,G$ are bounded, one has $a>1$ .

Assumption 2

We defer to the Supplementary Material for Conditions (B1)–(B4) on the underlying processes $X_{i}^{(g)}(t)$ . These conditions specify how data are sampled and smoothed.

Lemma 4

For a fixed number $G$ of trajectory types, under Assumptions 1 and 2,

\displaystyle\{(\widehat{\bm{A}}^{\top}\widehat{\bm{A}})_{ij}-E(\bm{A}^{\top}\bm{A})_{ij}\}/n=O_{p}(1/\sqrt{n}),\quad\text{for any}~{}i,j\leq p.

Further, $\|\widehat{\bm{a}}_{i}-{\bm{a}}_{i}\|_{2}^{2}\leq O_{P}(1/n)$ .

Lemma 4 indicates the convergence rate of estimating $\mathbb{E}(\bm{A}^{\top}\bm{A})_{ij}/n$ using $(\widehat{\bm{A}}^{\top}\widehat{\bm{A}})_{ij}$ is the same as using $({\bm{A}}^{\top}{\bm{A}})_{ij}$ . Combined with Lemma 3, the estimation error of functional structure factors is formally established.

Theorem 1

With $\widehat{\bm{F}}=\widehat{\bm{A}}\widehat{\bm{{B}}}\operatorname{diag}\left(\widehat{\lambda}_{1}^{-1}\ldots,\widehat{\lambda}_{{K}}^{-1}\right)$ and $\widehat{\bm{U}}=\widehat{\bm{A}}-\widehat{\bm{F}}\widehat{\bm{B}}^{\top}$ , we have

\displaystyle\left\{\begin{array}[]{l}\|\widehat{\bm{f}}_{i}^{\top}-\bm{f}_{i}^{\top}\|_{2}^{2}=O_{P}((\omega^{2}+1/n)\|\bm{a}_{i}\|_{2}^{2}+\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\|_{2}^{2}),\\ \|\widehat{\bm{u}}_{i}^{\top}-\bm{u}_{i}^{\top}\|_{2}^{2}=O_{P}((\omega^{2}+1/n)\|\bm{a}_{i}\|_{2}^{2}+\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\bm{B}^{\top}\|_{2}^{2}).\end{array}\right.

2.4 Truncation analysis and the relationship between $K$ and $m_{g}$

A practical issue when modeling functional data with sample instances is how to determine the number of basis functions using truncation, i.e., $m_{g}$ for $X_{i}^{(g)}(t)$ , and how this truncation would affect the estimate of functional factors $\bm{f}_{i}$ in (9). Usually, $m_{g}$ would be determined by the cumulative contribution to the variance from the corresponding functional scores, which tends to get an overestimated value of $m_{g}$ , namely, $\widehat{m}_{g}>m_{g}$ . To illustrate, assume $\widehat{m}_{1}>m_{1}$ and $\widehat{m}_{g}=m_{g}$ for $g=2,\ldots,G$ , and define $\tilde{\bm{a}}_{i}=(\bm{a}_{i}^{(1)\top},\tilde{\bm{a}}_{i}^{(1)\top},\bm{a}_{i}^{(2)\top},\ldots,\bm{a}_{i}^{(G)\top})^{\top}$ where $\tilde{\bm{a}}_{i}^{(g)\top}$ is a redundant variable for the factor model. In this case,

\displaystyle\operatorname{Cov}((\bm{a}_{i}^{(g)\top},\tilde{\bm{a}}_{i}^{(g)\top})^{\top})=\begin{pmatrix}\Sigma_{m_{g}}&\bm{0}\\ \bm{0}&\tilde{\Lambda}_{0}^{(g)}\end{pmatrix},

and

\displaystyle\operatorname{Cov}(\tilde{\bm{a}}_{i})=\begin{pmatrix}\Lambda_{m_{1}}&\bm{0}&\Sigma_{m_{1}m_{2}}&\cdots&\Sigma_{m_{1}m_{G}}\\ \bm{0}&\tilde{\Lambda}_{0}^{(1)}&\bm{0}&\cdots&\bm{0}\\ \Sigma_{m_{2}m_{1}}&\bm{0}&\Lambda_{m_{2}}&\cdots&\Sigma_{m_{2}m_{G}}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \Sigma_{m_{G}m_{1}}&\bm{0}&\Sigma_{m_{G}m_{2}}&\cdots&\Lambda_{m_{G}}\end{pmatrix}.

If $\xi^{\star}=(\xi_{1}^{\star\top},\xi_{2}^{\star\top})^{\top}$ is one of the first $K$ eigenvectors of $\operatorname{Cov}({\bm{a}}_{i})$ , $(\xi_{1}^{\star\top},\bm{0},\xi_{2}^{\star\top})^{\top}$ is an eigenvector of $\operatorname{Cov}(\tilde{\bm{a}}_{i})$ with the same eigenvalue, and the position of $\bm{0}$ in $(\xi_{1}^{\star\top},\bm{0},\xi_{2}^{\star\top})^{\top}$ corresponds to that of $\tilde{\Lambda}_{0}^{(1)}$ in $\operatorname{Cov}(\tilde{\bm{a}}_{i})$ . As $\widehat{\bm{F}}=\bm{A}\widehat{\bm{B}}\operatorname{diag}\left(\lambda_{1}^{-1}\ldots,\lambda_{K}^{-1}\right)$ , $\tilde{\bm{a}}_{i}^{(g)}$ will not affect $\widehat{\bm{F}}$ since the row of $\bm{B}$ corresponding to $\tilde{\bm{a}}_{i}^{(g)}$ is $\bm{0}$ when ${\widehat{K}}=K$ . Actually, in this overestimating situation, we can treat

\displaystyle\tilde{\bm{a}}_{i}=\begin{pmatrix}B_{1}\\ \bm{0}\\ B_{2}\end{pmatrix}\bm{f}_{i}+\begin{pmatrix}\bm{u}_{i}^{(1)}\\ \tilde{\bm{a}}_{i}^{(1)}\\ \bm{u}_{i}^{(C)}\end{pmatrix}=B^{\star}\bm{f}_{i}+\tilde{\bm{u}}_{i},

as the real model with $\bm{u}_{i}=(\bm{u}_{i}^{(1)}{}^{\top},\bm{u}_{i}^{(C)}{}^{\top})^{\top}$ and $\bm{B}=\left(B_{1}^{\top},B_{2}^{\top}\right)^{\top}$ .

Another case is that $m_{g}$ may be underestimated for some specific $g_{0}\leq G$ , especially when $\omega$ and $K$ are relatively small. In this case, a majority amount of variance of $X_{i}(t)$ is concentrated in ${{\tilde{\bm{\phi}_{0}}}{}^{(g_{0})}(t)}{}^{\top}\bm{f}_{i}$ by (10). When $\tilde{\bm{\phi}}_{0}^{(g_{0})}(t)$ are linearly independent for $g=1,\ldots,G$ , $\widehat{m}_{g_{0}}$ is determined as $K$ by the variance contribution criteria. However, when they are correlated, an underestimation of $m_{g_{0}}$ may not significantly affect the estimation of $\widehat{\bm{f}}_{i}$ . To illustrate, consider that the covariance matrix of $\bm{a}_{i}$ with $G=2$ as

\displaystyle{\operatorname{Cov}(\bm{a}_{i})}=\begin{pmatrix}\Lambda_{m_{1}}&\Sigma_{m_{1}m_{2}}\\ \Sigma_{m_{2}m_{1}}&\Lambda_{m_{2}}\end{pmatrix},

with $K<m_{1},m_{2}$ , and for convenience, we assume $\Lambda_{u}=\bm{I}_{p}$ . Since $B^{(g)}$ is an $m_{g}\times K$ matrix, $B^{(g)}B^{(g)}{}^{\top}$ has at most $K$ non-zero eigenvalues, denoted as $\nu_{j}^{(g)},j=1,\ldots,K$ . By plugging in $\Lambda_{m_{g}}=P^{(g)}{}^{\top}B^{(g)}B^{(g)}{}^{\top}P^{(g)}+\omega\bm{I}_{m_{g}}=\operatorname{diag}\{\nu_{1}^{(g)}+\omega,\ldots,\nu_{K}^{(g)}+\omega,\omega,\ldots,\omega\}$ , ${\operatorname{Var}(\bm{a}_{i})}$ is written as

	$\displaystyle{\operatorname{Cov}(\bm{a}_{i})}$	$\displaystyle=\begin{pmatrix}\Lambda_{m_{1(K)}}&0&\Sigma_{m_{1(K)}m_{2(K)}}&\Sigma_{m_{1(K)}m_{2(C)}}\\ 0&\Lambda_{m_{1(C)}}&\Sigma_{m_{1(C)}m_{2(K)}}&\Sigma_{m_{1(C)}m_{2(C)}}\\ \Sigma_{m_{2(K)}m_{1(K)}}&\Sigma_{m_{2(K)}m_{1(C)}}&\Lambda_{m_{2(K)}}&0\\ \Sigma_{m_{2(C)}m_{1(K)}}&\Sigma_{m_{2(C)}m_{1(C)}}&0&\Lambda_{m_{2(C)}}\end{pmatrix}$
		$\displaystyle\approx\begin{pmatrix}\Lambda_{m_{1(K)}}&0&\Sigma_{m_{1(K)}m_{2(K)}}&0\\ 0&0&0&0\\ \Sigma_{m_{2(K)}m_{1(K)}}&0&\Lambda_{m_{2(K)}}&0\\ 0&0&0&0\end{pmatrix},$

with $\Lambda_{m_{g}(K)}=\operatorname{diag}\{\tau_{1}^{(g)}+\omega,\ldots,\tau_{K}^{(g)}+\omega\}$ and $\Lambda_{m_{g}(C)}=\operatorname{diag}\{\omega,\ldots,\omega\}$ . Consequently, when $\omega\approx 0$ , underestimating $m_{g}$ will not significantly affect $\widehat{\bm{f}}_{i}$ as long as $\widehat{m}_{g}\geq K$ , which is easy to achieve in practice.

3 The functional factor augmentation selection model

3.1 The functional linear regression model

In this section, we address a multivariate functional linear regression model with correlated functional covariates, using the proposed fFAS in Section 2. To start with, consider a functional linear regression model where a scalar response $Y$ with $\mathbb{E}(Y)=\mu_{Y}$ is generated by a group of $G$ functional covariates $\{X^{(1)}(t),X^{(2)}(t),\ldots,X^{(G)}(t);t\in\mathcal{T}\}$ as

Y=\mu_{Y}+\sum_{g=1}^{G}\int_{\mathcal{T}}\beta^{*}{}^{(g)}(t)\left(X^{(g)}(t)-\mu^{(g)}(t)\right)dt+\epsilon,

(13)

where $\beta^{*}{}^{(g)}(\cdot)$ is the square-integrable regression parameter function, and $\epsilon$ is a random noise with a zero mean and a constant variance $\sigma^{2}$ . We note $\beta_{0}^{*}=\mu_{Y}$ as the intercept term. Thus, with given i.i.d samples $\{y_{i},X_{i}^{(g)}(t),g=1,\ldots,G;t\in\mathcal{T},i=1,\ldots,n$ } which have the same distribution as $\{Y,X^{(g)}(t),g=1,\ldots,G;t\in\mathcal{T}\}$ , the sample functional linear regression model is described as

\displaystyle y_{i}

\displaystyle=\beta_{0}^{*}+\sum_{g=1}^{G}\int_{\mathcal{T}}\beta^{*}{}^{(g)}(t)\left(X_{i}^{(g)}(t)-\mu^{(g)}(t)\right)dt+\epsilon_{i}.

(14)

To detect active functional covariates with correlation, we develop a functional factor augmentation selection model (fFASM) as follows. Without loss of generalizability, we use $X_{i}^{(g)}(\cdot)$ and $y_{i}$ as the centered functional covariates and scalar response variable, respectively, and accordingly (14) is equivalent to

y_{i}=\sum_{g=1}^{G}\int_{\mathcal{T}}\beta^{*}{}^{(g)}(t)X_{i}^{(g)}(t)dt+\epsilon_{i},

(15)

and can be further expanded by K-L expansion as

\displaystyle y_{i}

\displaystyle=\sum_{g=1}^{G}\sum_{j=1}^{m_{g}}a_{ij}^{(g)}\eta_{j}^{*}{}^{(g)}+\sum_{g=1}^{G}\int_{\mathcal{T}}\beta^{*}{}^{(g)}(t)e_{i}^{(g)}(t)dt+\epsilon_{i}=\bm{H}^{*}{}^{\top}\cdot\bm{a}_{i}\ +\tilde{\epsilon}_{i},

where $\bm{H}^{*}=\left(\bm{\eta}^{*}{}^{(1)}{}^{\top},\ldots,\bm{\eta}^{*}{}^{(G)}{}^{\top}\right)^{\top}\in\mathbb{R}^{p},\eta_{j}^{*}{}^{(g)}=\int_{\mathcal{T}}\beta^{*}{}^{(g)}(t){\gamma}_{j}^{(g)}(t)dt,j=1,\ldots,m_{g}$ , and $\tilde{\epsilon}_{i}=\epsilon_{i}+\sum_{g=1}^{G}\int_{\mathcal{T}}\beta^{*}{}^{(g)}(t)e_{i}^{(g)}(t)dt$ . For $g=1,\ldots,G$ , $\int_{\mathcal{T}}\beta^{*}{}^{(g)}(t)e_{i}^{(g)}(t)dt$ has a zero mean. Consequently, to select useful functional covariates $X_{i}^{(g)}(t)$ is to find $\beta^{*}{}^{(g)}(\cdot)$ such that $\beta^{*}{}^{(g)}(\cdot)\neq 0$ , which is further assumed as $\bm{\eta}^{*}{}^{(g)}\neq\bm{0}$ . Plugging in the proposed fFAS, one easily obtains

\displaystyle y_{i}=\bm{H}^{*}{}^{\top}\left(\bm{B}\bm{f}_{i}+\bm{u}_{i}\right)+\tilde{\epsilon}_{i},

or equivalently,

\displaystyle\bm{y}=\bm{F}\bm{B}^{\top}\bm{H}^{*}{}+\bm{U}\bm{H}^{*}{}+\tilde{\bm{\epsilon}},

(16)

and with $\widehat{\bm{F}}$ , $\widehat{\bm{B}}$ and $\widehat{\bm{U}}$ obtained by the fFAS, (16) is further equivalent to

\displaystyle\left(\bm{I}_{n}-\bm{P}_{\widehat{\bm{F}}}\right)\bm{y}

\displaystyle=\left(\bm{I}_{n}-\bm{P}_{\widehat{\bm{F}}}\right)\widehat{\bm{U}}\bm{H}^{*}{}+\left(\bm{I}_{n}-\bm{P}_{\widehat{\bm{F}}}\right)\left(\tilde{\bm{e}}+\tilde{\bm{\epsilon}}\right),

(17)

where $\bm{P}_{\widehat{\bm{F}}}={\widehat{\bm{F}}}({\widehat{\bm{F}}}^{\top}{\widehat{\bm{F}}})^{-1}{\widehat{\bm{F}}}^{\top}$ is the orthogonal projection matrix onto the column space $C({\widehat{\bm{F}}})$ , and $\tilde{\bm{e}}=(\tilde{e}_{1},\ldots,\tilde{e}_{n})^{\top}=\left(\bm{A}-\widehat{\bm{A}}+\mathbb{E}(\widehat{\bm{A}})\right)\bm{H}^{*}{}$ . Consequently to select useful functional covariates, the penalized loss function

\displaystyle L_{n}\left(\bm{y},\widehat{\bm{U}}\bm{H},\widehat{\bm{F}}\right)=\frac{1}{n}\left\|(\bm{I}-P_{\widehat{\bm{F}}})(\bm{y}-\widehat{\bm{U}}\bm{H})\right\|_{2}^{2}+\sum_{g=1}^{G}\sum_{k=1}^{m_{g}}J_{\lambda}\left(\eta_{k}^{(g)}\right)

is minimized with respect to $\bm{H}$ , where $J_{\lambda}=\lambda\cdot J(\cdot)$ is a penalty controlled by the parameter $\lambda$ and $J(\cdot)$ can be set as the popular penalties such as lasso, SCAD or MCP. Note that $\lambda$ can be selected using cross-validation. Hence, we successfully transform the problem from model selection with highly correlated functional covariates to model selection with weakly correlated or uncorrelated ones by lifting the space to higher dimension. As $\widehat{\bm{H}}=\mathop{\rm argmin}_{\bm{H}}L_{n}(\bm{y},\widehat{\bm{U}}\bm{H},\widehat{\bm{F}})$ is obtained, ${\beta}^{(g)}(t)$ is estimated as

\displaystyle\widehat{\beta}^{(g)}(t)=\sum_{j=1}^{\widehat{m}_{g}}\widehat{\eta}_{j}^{(g)}\widehat{\gamma}_{j}^{(g)}(t).

Accordingly, when $\widehat{\beta}^{(g)}(t)\neq 0$ , ${X_{i}}^{(g)}(t)$ is selected as a useful functional covariate. Note that the group selection method, such as the GM strategy by Aneiros, Novo & Vieu (2022), is not adopted in our case.

Also, this procedure may work even in the generalized linear context. Honestly,

$\displaystyle\bm{y}$	$\displaystyle=\left(\widehat{\bm{A}}-\mathbb{E}(\widehat{\bm{A}})\right)\bm{H}^{}+\left(\bm{A}-\widehat{\bm{A}}+\mathbb{E}(\widehat{\bm{A}})\right)\bm{H}^{}{}+\tilde{\bm{\epsilon}}$
	$\displaystyle=\widehat{\bm{F}}\widehat{\bm{B}}^{\top}\bm{H}^{}{}+\widehat{\bm{U}}\bm{H}^{}{}+\left(\bm{A}-\widehat{\bm{A}}+\mathbb{E}(\widehat{\bm{A}})\right)\bm{H}^{*}{}+\tilde{\bm{\epsilon}}$
	$\displaystyle=\widehat{\bm{F}}\widehat{\bm{B}}^{\top}\bm{H}^{}{}+\widehat{\bm{U}}\bm{H}^{}{}+\tilde{\bm{e}}+\tilde{\bm{\epsilon}},$	(18)

indicating that the explanatory variables are switched from $\bm{A}$ to $\widehat{\bm{A}}$ . In practice, when using the sample mean of $\widehat{\bm{A}}$ to substitute $\mathbb{E}(\widehat{\bm{A}})$ , $\mathbb{E}(\widehat{\bm{a}}_{i})\rightarrow 0$ under some regularization conditions (Kong et al., 2016). After centralizing $\widehat{\bm{A}}$ , by $\left(\widehat{\bm{A}}-\mathbb{E}(\widehat{\bm{A}})\right)\bm{H}^{*}{}=\widehat{\bm{F}}\widehat{\bm{B}}^{\top}\bm{H}^{*}{}+\widehat{\bm{U}}\bm{H}^{*}{}=\widehat{\bm{F}}\bm{\gamma}^{*}+\widehat{\bm{U}}\bm{H}^{*}{}$ , the unknown parameters are transformed into $(\bm{H},\bm{\gamma})$ , so that their corresponding covariates $\widehat{\bm{U}}$ and $\widehat{\bm{F}}$ are weakly correlated by introducing $\bm{\gamma}$ . Note that in linear cases, $\bm{\gamma}^{*}$ is further eliminated by using the projection matrix in (17). Consequently, for generalized linear models and samples without centralization, the loss function is updated as

	$\displaystyle L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)=$	$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[-y_{i}\left(\beta_{0}+\widehat{\bm{f}}_{i}^{\top}\bm{\gamma}+\widehat{\bm{u}}_{i}^{\top}\bm{H}\right)+b\left(\beta_{0}+\widehat{\bm{f}}_{i}^{\top}\bm{\gamma}+\widehat{\bm{u}}_{i}^{\top}\bm{H}\right)\right]$
		$\displaystyle+\sum_{g=1}^{G}\sum_{k=1}^{m_{g}}J_{\lambda}\left(\eta_{k}^{(g)}\right)+J_{\lambda}(\beta_{0}),$		(19)

where $\widehat{\bm{w}}_{i}=\left(1,\widehat{\bm{u}}_{i}^{\top},\widehat{\bm{f}}_{i}^{\top}\right)^{\top}$ , $\widehat{\bm{W}}=\left(\widehat{\bm{w}}_{1},\ldots,\widehat{\bm{w}}_{n}\right)$ and $\bm{\theta}=\left(\beta_{0},\bm{H}^{\top},\bm{\gamma}^{\top}\right)^{\top}$ , $\beta_{0}$ is the intercept term and $b(\cdot)$ is a known function, where $b(z)=z^{2}/2$ in linear models. The estimate is given by $\left(\widehat{\beta_{0}},\widehat{\bm{H}}^{\top},\widehat{\bm{\gamma}}^{\top}\right)^{\top}=\mathop{\rm argmin}_{\bm{\theta}}L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)$ .

3.2 Theoretical justifications for functional variable selection

In this section, the proposed method is theoretically investigated under the general linear model context in (3.1), and hence the linear model will be covered as a special case. To start with, some notations and assumptions are introduced. Recall $p=\sum_{g=1}^{G}m_{g}$ , and define $p_{1}=1+p$ , $\bm{H}_{1}=(\beta_{0},\bm{H}^{\top})^{\top}$ , $\bm{H}_{1}^{*}=(\beta_{0}^{*},\bm{H}^{*}{}^{\top})^{\top}$ , $\bm{\theta}^{*}=(\beta_{0}^{*},\bm{H}^{*}{}^{\top},(\bm{B}^{\top}\bm{H}^{*})^{\top})^{\top}$ , $S=\operatorname{supp}(\bm{\theta}^{*})$ , $S_{1}=\operatorname{supp}(\bm{H}_{1}^{*})$ , and $S_{2}=[p_{1}+K]\backslash S$ . Suppose $\widehat{\bm{F}}$ and $\widehat{\bm{U}}$ are obtained given $\bm{A}$ . We write

\displaystyle X_{i}^{(g)}(t)

\displaystyle=X_{f,i}^{(g)}(t)+X_{u,i}^{(g)}(t)+e_{0i}^{(g)}(t),

(20)

with $X_{f,i}^{(g)}(t)=\bm{\gamma}^{(g)}(t){}^{\top}B^{(g)}\bm{f}_{i}$ and $X_{u,i}^{(g)}(t)={\bm{\gamma}^{(g)}(t)}{}^{\top}\bm{u}_{i}^{(g)}$ .

Assumption 3

(Smoothness). $b(z)\in C^{3}(\mathbb{R})$ , i.e., for some constants $M_{2}$ and $M_{3}$ , $0\leq b^{\prime\prime}(z)\leq M_{2}$ and $\left|b^{\prime\prime\prime}(z)\right|$ $\leq M_{3},\forall z$ .

Assumption 4

(Restricted strong convexity and irrepresentable condition). For ${\bm{W}}=\left({\bm{w}}_{1},\ldots,{\bm{w}}_{n}\right)$ where ${\bm{w}}_{i}=\left(1,{\bm{u}}_{i}^{\top},{\bm{f}}_{i}^{\top}\right)^{\top}$ , there exist $\kappa_{2}>\kappa_{\infty}>0$ and $\tau\in(0,0.5)$ such that

	$\displaystyle\text{(Convexity)}\left\\|\left[\nabla_{SS}^{2}L_{n}\left(\bm{y},\bm{W}\boldsymbol{\theta}^{*}\right)\right]^{-1}\right\\|_{\ell}\leq\frac{1}{4\kappa_{\ell}},\quad\text{ for }\ell=2\text{ and }\infty,$		(21)
	$\displaystyle\text{(Irrepresentable condition)}\left\\|\nabla_{S_{2}S}^{2}L_{n}\left(\bm{y},\bm{W}\boldsymbol{\theta}^{}\right)\left[\nabla_{SS}^{2}L_{n}\left(\bm{y},\bm{W}\boldsymbol{\theta}^{}\right)\right]^{-1}\right\\|_{\infty}\leq 1-2\tau.$		(22)

Assumption 5

(Estimation of factor model). $\|\gamma_{j}^{(g)}(\cdot)\|_{\infty}\leq M_{\gamma}$ for $j=1,\ldots,m_{g},~{}g=1,\ldots,G$ , and $\|X_{f,i}^{(g)}(\cdot)\|_{\infty}\leq\frac{1}{2}M_{\gamma}M_{0}\sqrt{\lambda_{1}^{(B)}}$ and $\|X_{u,i}^{(g)}(\cdot)\|_{\infty}\leq\frac{1}{2}M_{\gamma}M_{0}$ for some constant $M_{0}>2$ . In addition, there exists a $K\times K$ nonsingular matrix $\bm{V}_{0}$ , and $\bm{V}=\left(\begin{array}[]{cc}\bm{I}_{p_{1}}&\bm{0}_{{p_{1}}\times K}\\ \bm{0}_{K\times{p_{1}}}&\bm{V}_{0}\end{array}\right)$ such that for $\overline{\bm{W}}=\widehat{\bm{W}}\bm{V}$ , we have $\|\overline{\bm{W}}-\bm{W}\|_{\mathop{\rm max}}\leq\frac{M_{0}}{2}$ and $\mathop{\rm max}_{j\in[p_{1}+K]}$ $\left(\frac{1}{n}\sum_{i=1}^{n}\left|\bar{w}_{ij}-w_{ij}\right|^{2}\right)^{1/2}\leq\frac{2\kappa_{\infty}\tau}{3M_{0}M_{2}|S|}$ , with $w_{ij}$ and $\bar{w}_{ij}$ being the $(i,j)$ -th element of $\mathbf{W}$ and $\overline{\mathbf{W}}$ .

Assumption 3 holds for a large family of generalized linear models. For example, linear model has $b(z)=\frac{1}{2}z^{2}$ , $M_{2}=1$ , and $M_{3}=0$ ; logistic model has $b(z)=\log\left(1+e^{z}\right)$ and finite $M_{2},M_{3}$ . Assumption 4 is easily satisfied with a small matrix and holds with high probability as long as $\mathbb{E}\left[\nabla^{2}L_{n}\left(\bm{y},\bm{W}\boldsymbol{\theta}^{*}\right)\right]$ satisfies similar conditions by standard concentration inequalities (Merlevède et al., 2011). Assumption 5 indicates the norm $\|\bm{W}\|_{\mathop{\rm max}}$ when $X_{i}^{(g)}(\cdot)$ and $\varepsilon_{i}^{(g)}(\cdot)$ is bounded, and

	$\displaystyle\\|\widehat{\bm{f}}_{i}^{\top}-\bm{f}_{i}^{\top}\\|_{2}^{2}$	$\displaystyle=O(\omega^{2}+\\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\\|_{2}^{2}+1/n),$
	$\displaystyle\\|\widehat{\bm{u}}_{i}^{\top}-\bm{u}_{i}^{\top}\\|_{2}^{2}$	$\displaystyle=O(\omega^{2}+\\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\bm{B}^{\top}\\|_{2}^{2}+1/n),$

will be satisfied with a high probability when $n$ is large and $\omega$ is small. Furthermore, with $J(\cdot)=|\cdot|$ and $\bm{\widehat{H}}$ obtained by minimizing (3.1), the following result is established.

Theorem 2

Suppose Assumptions 3-5 holds. Define $M=M_{0}^{3}M_{3}|S|^{3/2}$ , and

\varepsilon^{*}=\mathop{\rm max}_{j\in[p_{1}+K]}\left|\frac{1}{n}\sum_{i=1}^{n}\bar{w}_{ij}\left[-y_{i}+b^{\prime}\left({\bm{a}}_{i}^{T}\bm{H}^{*}+\beta_{0}^{*}\right)\right]\right|.

If $\frac{7\varepsilon^{*}}{\tau}<\lambda<\frac{\kappa_{2}\kappa_{\infty}\tau}{12M\sqrt{|S|}}$ , then $\operatorname{supp}(\widehat{\bm{H}})\subseteq\operatorname{supp}\left(\bm{H}^{*}\right)$ and

\left\|\widehat{\bm{H}}-\bm{H}^{*}\right\|_{\infty}\leq\frac{6\lambda}{5\kappa_{\infty}},\quad\left\|\widehat{\bm{H}}-\bm{H}^{*}\right\|_{2}\leq\frac{4\lambda\sqrt{|S|}}{\kappa_{2}},\quad\left\|\widehat{\bm{H}}-\bm{H}^{*}\right\|_{1}\leq\frac{6\lambda|S|}{5\kappa_{\infty}}.

In addition, if $\varepsilon^{*}<\frac{\kappa_{2}\kappa_{\infty}\tau^{2}}{12CM\sqrt{|S|}}$ and $\min\left\{\left|{\bm{H}_{1}^{*}}_{j}\right|:{\bm{H}_{1}^{*}}_{j}\neq 0,j\in[p_{1}]\right\}>\frac{6C\varepsilon^{*}}{5\kappa_{\infty}\tau}$ hold for some $C>7$ , where ${\bm{H}_{1}^{*}}_{j}$ means the j-th element of ${\bm{H}_{1}^{*}}$ , then by taking $\lambda\in\left(\frac{7}{\tau}\varepsilon^{*},\frac{C}{\tau}\varepsilon^{*}\right)$ the sign consistency $\operatorname{sign}(\widehat{\bm{H}})=\operatorname{sign}\left(\bm{H}^{*}\right)$ is achieved.

Theorem 2 guarantees the selection consistency of the functional covariates by developing sign consistency under some mild conditions. Furthermore, it demonstrates the relationship between $\left\|\widehat{\bm{H}}-\bm{H}^{*}\right\|$ and $\lambda$ whose value depends on $\varepsilon^{*}$ and $\tau$ , where $\varepsilon^{*}$ comes from the first-order partial derivative of the empirical loss function (without penalty) when the parameters are known and $\tau$ satisfying a generalized irrepresentable condition from (22) (Lee et al., 2015). When $\varepsilon^{*}$ is small, selecting an appropriate $\lambda$ leads to satisfactory performance in the model estimation. Furthermore, when $\bm{A}$ is not available, $\widehat{\bm{A}}$ will be estimated and employed. Accordingly, the sign consistency is further guaranteed by showing that the loss function (3.1) based on $\widehat{\bm{A}}$ will have an asymptotic property as follows.

Theorem 3

Let $L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)|_{\bm{A}}$ be the loss function in (3.1), based on which $\widehat{\bm{f}}_{i}$ and $\widehat{\bm{u}}_{i}$ are obtained using the true $\bm{A}$ , and $L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)|_{\widehat{\bm{A}}}$ using ${\widehat{\bm{A}}}$ . Then for any fixed $\bm{\theta}$ , by Assumption 3 and Lemma 4,

\displaystyle\left|L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)|_{\bm{A}}-L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)|_{\widehat{\bm{A}}}\right|\rightarrow 0\text{ in probability}.

Corollary 1

Suppose $L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)$ has a unique global minimum. Then $L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)|_{\bm{A}}$ and $L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)|_{\widehat{\bm{A}}}$ will have an identical global minimizer when $n$ is sufficiently large, so that the estimates from $L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)|_{\widehat{\bm{A}}}$ will also have sign consistency.

4 Simulation Studies

In this section, the performance of the proposed fFASM estimator is examined for correlated functional data using simulation studies. To start with, the functional covariates $\{X_{i}^{(g)}(t)\}_{g=1}^{G}$ are generated by $X_{i}^{(g)}(t)=\bm{a}_{0,i}^{(g)}{}^{\top}\cdot\bm{\phi}_{0,m_{g}}^{(g)}(t)+\varepsilon_{i}(t)$ with $\varepsilon_{i}(\cdot)$ independently and identically (i.i.d) drawn from $N(0,0.25)$ , and sampled over 51 uniformly distributed grid points for $t$ from $\mathcal{T}=[0,1]$ with interval lengths of 0.02, where $G$ is set as 20, 50, 100 and 150, respectively. The basis functions $\bm{\phi}_{0,m_{g}}^{(g)}(t)=\bm{\phi}_{m}(t)=(\phi_{1}(t),\ldots,\phi_{m}(t))^{\top}$ are Fourier basis with $m$ set as 10. To blend correlations into functional covariates $X_{i}^{(g)}(t)$ , two scenarios are considered when $\bm{a}_{i}$ is generated:

•

Scenario I: the factor model case, i.e., $\bm{a}_{0,i}=\left(\bm{a}_{0,i}^{(1)\top}{},\ldots,\bm{a}_{0,i}^{(G)\top}{}\right)^{\top}=\bm{B}\bm{f}_{i}+\bm{u}_{i},$ where the elements of $\bm{f}_{i}\in\mathbb{R}^{K}$ , $\bm{u}_{i}\in\mathbb{R}^{Gm}$ and $\bm{B}\in\mathbb{R}^{Gm\times K}$ are generated from $N(0,25)$ , $N(0,1)$ , and $N(0,1)$ , respectively. The true number of factors $K$ is set as $1,2,\ldots,6$ , respectively;
•

Scenario II: the equal correlation case, i.e., $\bm{a}_{0,i}$ are generated from a multivariate normal distribution $\bm{a}_{0,i}\sim N(\bm{0},\Sigma_{Equal})$ where $\Sigma_{Equal}$ has diagonal elements as 1 and off-diagonal elements as $\rho$ . The true value of $\rho$ is set as $0,0.1,0.2\ldots,0.9$ , respectively.

Furthermore, the true functional coefficients $\beta^{*(g)}$ in (15) are generated $\beta^{*(1)}(t)=\beta^{*(2)}(t)=\beta^{*(3)}(t)=\beta^{*(4)}(t)=(1,1/2^{2},\ldots,1/m^{2})\cdot\bm{\phi}_{m}(t)$ (the Harmonic attenuation type), $\beta^{*(5)}(t)=\beta^{*(6)}(t)=(0,1/2^{2},0,\ldots,0)\cdot\bm{\phi}_{m}(t)$ (the weak single signal type), and $\beta^{*(g)}(t)=0$ for $g=7,8,\ldots,G$ . Accordingly, the response $y_{i}$ is generated by (15) with $\epsilon_{i}\sim N(0,0.1)$ .

We compare the performance of the proposed fFASM method (MCP penalty) with other functional variable selection methods, including the MCP (MCP) and the group MCP (grMCP) methods on the synthetic data in different scenarios. To evaluate the performance, three different measurements are adopted, i.e, the model size defined as the cardinality of the set $\left\{g:\widehat{\beta}^{(g)}(t)\neq 0\right\}$ , the integrated mean squared error (IMSE)

IMSE=\sum_{g=1}^{G}\mathbb{E}\left(\left\|\widehat{\beta}^{(g)}(t)-\beta^{*(g)}(t)\right\|^{2}_{L^{2}}\right)=\sum_{g=1}^{G}\mathbb{E}\int_{\mathcal{T}}\left(\widehat{\beta}^{(g)}(t)-\beta^{*(g)}(t)\right)^{2}dt,

and the true positive rate (TPR)

\text{TPR}=\frac{\text{TP}}{\text{TP}+\text{FN}},

where TP represents the number of correctly predicted nonzero instances of $\widehat{\beta}^{(g)}(t)$ (events that are actually positive and predicted as positive), and FN for that of instances that are actually nonzero but are predicted as zero. Note that a model with the TPR of 1 and the model size of 6 indicates perfect recovery. The proposed fFASM method employs the MCP regularization, and the hyperparameters are tuned using cross-validation by minimizing the cross-validation error with a randomly selected subset of $[n/3]$ training samples, where the sample size $n$ is set as 100. The whole experiment is repeated for 500 times, and the averaged results are reported.

Refer to caption — Figure 1: Estimation and selection performance of three methods with different numbers of factors $K$ for Scenario I. In each row, four subfigures are displayed with $G=20,50,100,150$ , respectively, and in each column with a given $K$ , the IMSE (with 0.5 standard deviations), the TPR, and model size are presented, respectively.

Figure 1 shows the estimation and selection performance for Scenario I (the factor model structure case) using different methods. As is easily observed, the proposed fFASM method significantly and consistently outperforms MCP and grMCP in the sense that the proposed fFASM has the smallest IMSE for each $K$ , demonstrating a more accurate estimation of the functional coefficients. Furthermore, the TPR for the proposed fFASM method is always greater than those of MCP and grMCP as with the increase of $K$ , and MCP and grMCP even show dramatic drops when $K$ increases. In terms of the model size, the proposed fFASM tends to select more functional covariates into the model, slightly greater than the true model size 6, while the MCP and grMCP methods tend to select irrelevant functional covariates into the model when looking deeper at the selection results, and the grMCP method tends to be more conservative in selecting functional covariates than MCP.

Furthermore, Figure 2 shows the estimation and selection results in Scenario II (the equal correlation case). The IMSE values of the proposed fFASM method turn out to remain at a low level when the correlation $\rho$ increases in each subfigure in the first row, indicating a robust estimation performance from fFASM, while MCP and grMCP tend to increase, demonstrating less satisfactory estimation performance. Although the two methods show slightly better results when functional covariates are of small correlation with each other, this is expected since the equal correlation indicates a small number of factors, which may be captured easily. In terms of TPR, the fFASM method is again consistently close to 1 when the correlation varies, while TPRs for the competitor methods drop dramatically. Also, when the number of functional covariates $G$ increases, the TPRs from all three methods turn out to drop without surprise. When looking at the model size, the proposed fFASM remains stable, which may slightly overestimate the model size, while the competitor models show a much less stable model size of little consistency.

Additionally, the selection frequencies of the two types of functional coefficients $\beta^{*(g)}(t)$ are examined, namely the Harmonic attenuation type (Type1) and the weak single signal type (Type2), displayed in Figure 3 and 4, respectively. Note that the number of true Type1 functional covariates is 4 and that of Type2 is 2 in both scenarios. On one hand, Figure 3 shows that the selection frequency of the proposed fFASM method decreases first and then increases when the number of factors $K$ increases in both two types of functional covariates, while those for MCP and grMCP tend to be 0, especially in the Type2 setting where the two competitor methods show lower starting points. Note that the selection frequencies of the proposed method for both two types are not extraordinarily high, since the functional linear model is somewhat lack of fit due to the complexity of the true model. On the other hand, in Scenario II, Figure 4 shows a much improved selection consistency of the proposed method than grMCP and MCP for both two types, even though Type2 may still select fewer covariates because of slight violation of norm requirement in Theorem 2.

5 Real Data Application

5.1 Effects of macroeconomic covariates on lifetime expectancy

In this section, the effects of macroeconomic covariates on national expected lifetime are explored for European and American countries using the open-sourced EPS data platform ¹¹1https://www.epsnet.com.cn/. The data have been collected and documented annually over a span of 21 years, from 2000 to 2020, for 40 countries or regions. Our focus is the average lifetime expectancy as a scalar dependent variable against 33 macroeconomic functional covariates, such as gross domestic product, healthcare expenditure, and educational attainment. Functional covariates are correlated with each other. For example, higher levels of educational attainment often correlate with higher GDP, since a more educated workforce may usually contribute to greater economic productivity. To estimate lifetime expectancy and find out useful functional covariates, the proposed fFASM method, the MCP and grMCP methods are employed. Specifically, to recover the functional trajectories, year points from 2000 to 2022 are used as grid points. During the study, the collected data for 30 randomly selected countries (or regions) out of 40 are used as the train set to build the model, and the remaining 10 are used as the test set. The whole procedure is repeated 200 times, and the average out-of-sample $R^{2}$ and model size are reported in Table 1 and the top five selected functional covariates in Table 2.

The results from both the prediction performance in Table 1 and the variable selection in Table 2 may strongly support the benefits of the proposed fFASM. The average out-of-sample $R^{2}$ from the fFASM is the greatest among the three competitor models, indicating its competitive prediction accuracy. Further, the variables selected by fFASM are closely related to economic development and social well-being, providing insights into education and labor market conditions, which may show significant influence on the average life expectancy of a nation. Although all three methods successfully selected the key variable “Employment” with the highest selection frequencies, “Adult literacy rate” and “total fertility rate” are both selected with relatively large frequencies by fFASM and MCP, while grMCP fails to show high selection frequencies of them. Additionally, the fFASM method selected “Gross national income per capita” and “Unemployment rate”, which are of little correlation with other variables, while MCP selected “Labor force” highly correlated with “Employment” and “Unemployment rate”. The grMCP only focused on “Employment” and ignored almost all other variables with selection rates lower than 5% .

average out-of-sample $R^{2}$			average model size
fFASM	MCP	grMCP	fFASM	MCP	grMCP
0.288	0.223	0.227	1.95	1.53	1.04

Table 1: The average out-of-sample

R^{2}

and model size for EPS data.

Rank	fFASM
Rank	Name	Frequency
1	Employment (million people)	160
2	Adult literacy rate (age 15 and above) (%)	56
3	Total fertility rate	54
4	Unemployment rate (% of total labor force)	21
5	Gross national income per capita (current USD)	17
Rank	MCP
Rank	Name	Frequency
1	Employment (million people)	166
2	Total fertility rate	57
3	Adult literacy rate (age 15 and above, %)	37
4	Labor force (people)	16
5	Unemployment rate (% of total labor force)	7
Rank	grMCP
Rank	Name	Frequency
1	Employment (million people)	187
2	Prevalence of undernourishment (% of total population)	7
3	Adult literacy rate (age 15 and above) (%)	3
4	Gross national income per capita (current USD)	2
5	Per capita health expenditure (current USD)	1

Table 2: Top 5 variables with highest selection frequency.

5.2 Prediction of sales areas of commercial houses

In this section, we focus on the prediction of annual average sales area of commercial houses in each province in China. Data are collected from the National Bureau of Statistics of China²²2https://data.stats.gov.cn/ on a monthly basis for 31 administrative regions of mainland China in the year 2022 and 2023. The dependent variable is the annual average sales area of commercial houses against 60 functional covariates, including cement, medium tractors, engines, aluminum materials, and other variables related to industrial output. To recover the functional trajectories, 24 grid points are predetermined on a monthly basis from 2022 to 2023. To predict the response, the proposed fFASM, the MCP and grMCP methods are employed. During the study, the collected data for 20 randomly selected provinces are used as the training set to build the model, and the remaining 11 as the test set to evaluate the prediction performance, with 200 repetitions of the whole procedure. The average out-of-sample $R^{2}$ and the average model sizes are reported in Table 3, and the most frequently selected functional covariates in Table 4.

From Table 3, the proposed fFASM method is found to show a better performance than MCP and grMCP in terms of a larger average out-of-sample $R^{2}$ and a slightly lower average model size. In Table 4, all three methods select “Cement” as the most frequently selected functional covariates out of 200 repetitions, which is expected as the most important raw material required for house construction. Similarly, the functional covariate “Aluminum Materials” is selected among top five. Specifically, the fFASM and MCP methods select the functional covariates “Medium Tractors” and “Engine”, which are closely related to construction machinery used in the house building industry, while grMCP selects “Mechanized Paper”, “Hydropower Generation” and “Phosphate Rock”, which seem to have limited connection with house construction. This may explain why grMCP has a larger model size but a relatively smaller average out-of-sample performance.

average out-of-sample $R^{2}$			average model size
fFASM	MCP	grMCP	fFASM	MCP	grMCP
0.612	0.601	0.355	1.35	2.20	2.55

Table 3: The average out-of-sample

R^{2}

and model size for sales area of commercial houses data.

Rank	fFASM
Rank	Name	Frequency
1	Cement	198
2	Medium Tractors	13
3	Engines	12
3	Aluminum Materials	12
5	Chemical Pesticides (Active Ingredients)	4
Rank	MCP
Rank	Name	Frequency
1	Cement	199
2	Medium Tractors	48
3	Aluminum Materials	19
4	Computers	18
4	Engines	18
Rank	grMCP
Rank	Name	Frequency
1	Cement	137
2	Aluminum Materials	97
3	Mechanized Paper	68
4	Hydropower Generation	30
5	Phosphate Rock	20

Table 4: Top 5 variables with highest selection frequency.

6 Conclusion

In this article, a novel functional factor augmentation structure (fFAS) is proposed to capture associations for correlated functional processes, and further a functional factor augmentation selection model (fFASM) is developed to select useful functional covariates in high dimensions with correlated functional covariates. Note only is the correlation between functional covariates addressed without assuming an explicit covariance structure, theoretical properties of the estimated functional factors are established. We primarily discuss the rationale for constructing a fFAS, how to estimate the fFAS and its estimation error, and the impact of truncating functional data on the validity and estimation of the factor model. Due to the unique characteristics of functional data, the assumed factor model and the actually estimated factor model may differ.

Numerical investigations on both simulated and real datasets support the superior performance of the proposed fFASM method. It is found that our proposed method performs better than general functional data variable selection methods when dealing with the variable selection problem of correlated multivariate functional covariates. A practical issue may be how to determine the model size, as our method may slightly select more functional covariates in simulation studies, which may be a trade-off for modeling functional processes with correlations in high dimensions.

References

(1)
Ahn & Horenstein (2013) Ahn, S. C. & Horenstein, A. R. (2013), ‘Eigenvalue ratio test for the number of factors’, Econometrica 81(3), 1203–1227.
Aneiros, Horová, Hušková & Vieu (2022) Aneiros, G., Horová, I., Hušková, M. & Vieu, P. (2022), ‘On functional data analysis and related topics’, Journal of Multivariate Analysis 189, 104861.
Aneiros, Novo & Vieu (2022) Aneiros, G., Novo, S. & Vieu, P. (2022), ‘Variable selection in functional regression models: A review’, Journal of Multivariate Analysis 188, 104871. 50th Anniversary Jubilee Edition.
https://www.sciencedirect.com/science/article/pii/S0047259X21001494
Bai (2003) Bai, J. (2003), ‘Inferential theory for factor models of large dimensions’, Econometrica 71(1), 135–171.
Bai & Ng (2002) Bai, J. & Ng, S. (2002), ‘Determining the number of factors in approximate factor models’, Econometrica 70(1), 191–221.
Cardot et al. (2003) Cardot, H., Ferraty, F. & Sarda, P. (2003), ‘Spline estimators for the functional linear model’, Statistica Sinica pp. 571–591.
Castellanos et al. (2015) Castellanos, L., Vu, V. Q., Perel, S., Schwartz, A. B. & Kass, R. E. (2015), ‘A multivariate gaussian process factor model for hand shape during reach-to-grasp movements’, Statistica Sinica 25(1), 5.
Centofanti et al. (2021) Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B. & Vantini, S. (2021), ‘Functional regression control chart’, Technometrics 63(3), 281–294.
Chen et al. (2022) Chen, C., Guo, S. & Qiao, X. (2022), ‘Functional linear regression: dependence and error contamination’, Journal of Business & Economic Statistics 40(1), 444–457.
Chen et al. (2011) Chen, D., Hall, P. & Müller, H.-G. (2011), ‘Single and multiple index functional regression models with nonparametric link’.
Chen et al. (2021) Chen, L., Wang, W. & Wu, W. B. (2021), ‘Dynamic semiparametric factor model with structural breaks’, Journal of Business & Economic Statistics 39(3), 757–771.
Chiou & Müller (2016) Chiou, J.-M. & Müller, H.-G. (2016), ‘A pairwise interaction model for multivariate functional and longitudinal data’, Biometrika 103(2), 377–396.
Cuevas (2014) Cuevas, A. (2014), ‘A partial overview of the theory of statistics with functional data’, Journal of Statistical Planning and Inference 147, 1–23.
Fan et al. (2020) Fan, J., Ke, Y. & Wang, K. (2020), ‘Factor-adjusted regularized model selection’, Journal of Econometrics 216(1), 71–85.
Fan & Li (2001) Fan, J. & Li, R. (2001), ‘Variable selection via nonconcave penalized likelihood and its oracle properties’, Journal of the American statistical Association 96(456), 1348–1360.
Fan et al. (2013) Fan, J., Liao, Y. & Mincheva, M. (2013), ‘Large covariance estimation by thresholding principal orthogonal complements’, Journal of the Royal Statistical Society Series B: Statistical Methodology 75(4), 603–680.
Fang et al. (2015) Fang, L., Zhao, H., Wang, P., Yu, M., Yan, J., Cheng, W. & Chen, P. (2015), ‘Feature selection method based on mutual information and class separability for dimension reduction in multidimensional time series for clinical data’, Biomedical Signal Processing and Control 21, 82–89.
Ferraty (2006) Ferraty, F. (2006), Nonparametric functional data analysis, Springer.
Gonzalez-Vidal et al. (2019) Gonzalez-Vidal, A., Jimenez, F. & Gomez-Skarmeta, A. F. (2019), ‘A methodology for energy multivariate time series forecasting in smart buildings based on feature selection’, Energy and Buildings 196, 71–82.
Hall & Horowitz (2007) Hall, P. & Horowitz, J. L. (2007), ‘Methodology and convergence rates for functional linear regression’.
Hays et al. (2012) Hays, S., Shen, H. & Huang, J. Z. (2012), ‘Functional dynamic factor models with application to yield curve forecasting’, The Annals of Applied Statistics pp. 870–894.
Hörmann et al. (2015) Hörmann, S., Kidziński, Ł. & Hallin, M. (2015), ‘Dynamic functional principal components’, Journal of the Royal Statistical Society Series B: Statistical Methodology 77(2), 319–348.
Hörmann & Kokoszka (2010) Hörmann, S. & Kokoszka, P. (2010), ‘Weakly dependent functional data’.
Htun et al. (2023) Htun, H. H., Biehl, M. & Petkov, N. (2023), ‘Survey of feature selection and extraction techniques for stock market prediction’, Financial Innovation 9(1), 26.
Jiménez et al. (2020) Jiménez, F., Palma, J., Sánchez, G., Marín, D., Francisco Palacios, M. & Lucía López, M. (2020), ‘Feature selection based multivariate time series forecasting: An application to antibiotic resistance outbreaks prediction’, Artificial Intelligence in Medicine 104, 101818.
https://www.sciencedirect.com/science/article/pii/S0933365719306608
Kong et al. (2016) Kong, D., Xue, K., Yao, F. & Zhang, H. H. (2016), ‘Partially functional linear regression in high dimensions’, Biometrika 103(1), 147–159.
Lee et al. (2015) Lee, J. D., Sun, Y. & Taylor, J. E. (2015), ‘On model selection consistency of regularized m-estimators’.
Li & Hsing (2007) Li, Y. & Hsing, T. (2007), ‘On rates of convergence in functional linear regression’, Journal of Multivariate Analysis 98(9), 1782–1804.
Matsui & Konishi (2011) Matsui, H. & Konishi, S. (2011), ‘Variable selection for functional regression models via the l1 regularization’, Computational Statistics & Data Analysis 55(12), 3304–3310.
Merlevède et al. (2011) Merlevède, F., Peligrad, M. & Rio, E. (2011), ‘A bernstein type inequality and moderate deviations for weakly dependent sequences’, Probability Theory and Related Fields 151, 435–474.
Müller & Yao (2008) Müller, H.-G. & Yao, F. (2008), ‘Functional additive models’, Journal of the American Statistical Association 103(484), 1534–1544.
Nti et al. (2019) Nti, K. O., Adekoya, A. & Weyori, B. (2019), ‘Random forest based feature selection of macroeconomic variables for stock market prediction’, American Journal of Applied Sciences 16(7), 200–212.
Peng et al. (2005) Peng, R. D., Dominici, F., Pastor-Barriuso, R., Zeger, S. L. & Samet, J. M. (2005), ‘Seasonal analyses of air pollution and mortality in 100 us cities’, American journal of epidemiology 161(6), 585–594.
Petersen (2024) Petersen, A. (2024), ‘Mean and covariance estimation for discretely observed high-dimensional functional data: Rates of convergence and division of observational regimes’, Journal of Multivariate Analysis p. 105355.
Ramsay & Silverman (2005) Ramsay, J. & Silverman, B. (2005), Functional Data Analysis, Springer Series in Statistics, Springer.
Shi (1997) Shi, L. (1997), ‘Local influence in principal components analysis’, Biometrika 84(1), 175–186.
Tibshirani (1996) Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society Series B: Statistical Methodology 58(1), 267–288.
Yang & Ling (2024) Yang, S. & Ling, N. (2024), ‘Robust estimation of functional factor models with functional pairwise spatial signs’, Computational Statistics pp. 1–24.
Yao et al. (2005) Yao, F., Müller, H.-G. & Wang, J.-L. (2005), ‘Functional linear regression analysis for longitudinal data’, The Annals of Statistics 33(6), 2873 – 2903.
Yuan & Cai (2010) Yuan, M. & Cai, T. T. (2010), ‘A reproducing kernel hilbert space approach to functional linear regression’.
Zhang (2010) Zhang, C.-H. (2010), ‘Nearly unbiased variable selection under minimax concave penalty’, Annals of statistics 38(2), 894–942.

SUPPLEMENTARY MATERIAL

7 Regularity Condition

In this section, we introduce some regularity conditions (Kong et al. 2016). Without loss of generality, we assume that $\{X^{(g)}(\cdot),g=1,\ldots,G\}$ have been centred to have zero mean. With $W_{i}^{(g)}(t_{ij})=X^{(g)}_{i}(t_{ij})+\epsilon_{i}^{(g)}(t_{ij})$ , for definiteness, we consider the local linear smoother for each set of subjects using bandwidths $\{h_{i}^{(g)},g=1,\ldots,G\}$ , and denote the smoothed trajectories by $\widehat{X}_{i}^{(g)}(\cdot)$ .

Condition (B1) consists of regularity assumptions for functional data. Condition (B2) is standard for local linear smoother, (B3)–(B4) concern how the functional predictors are sampled and smoothed.

(B1)

For $g=1,\ldots,G$ , for any $C>0$ there exists an $\epsilon>0$ such that

\sup_{t\in\mathcal{T}}\left[\mathbb{E}\left\{\left|X^{(g)}(t)\right|^{C}\right\}\right]<\infty,\quad\sup_{s,t\in\mathcal{T}}\left(\mathbb{E}\left[\left\{|s-t|^{-\epsilon}\left|X^{(g)}(s)-X^{(g)}(t)\right|\right\}^{C}\right]\right)<\infty.

For each integer $r\geq 1,\left(\tau_{k}^{(g)}\right)^{-r}E\left({a_{k}^{(g)}}\right)^{2r}$ is bounded uniformly in $k$ .

(B2)

For $g=1,\ldots,G,X^{(g)}(\cdot)$ is twice continuously differentiable on $\mathcal{T}$ with probability 1 , and $\int\mathbb{E}\left\{X^{(g)}{}^{\prime\prime}(t)\right\}^{4}dt<\infty$ , where $X^{(g)}{}^{\prime\prime}(\cdot)$ denotes the second derivative of $X^{(g)}(\cdot)$ . The following condition concerns the design on which $X_{i}^{(g)}(\cdot)$ is observed and the local linear smoother $\widehat{X}_{i}^{(g)}(\cdot)$ . When a function is said to be smooth, we mean that it is continuously differentiable to an adequate order.
(B3)

For $g=1,\ldots,G,\left\{t_{ij}^{(g)},j=1,\ldots,n_{i}^{(g)}\right\}$ are considered deterministic and ordered increasingly for $i=1,\ldots,n$ . There exist densities $p_{i}^{(g)}$ uniformly smooth over $\mathcal{T}$ , satisfying $\int_{\mathcal{T}}p_{i}^{(g)}(t)dt=1$ and $0<c_{1}<\inf_{i}\left\{\inf_{t\in T}p_{i}^{(g)}(t)\right\}<\sup_{i}\left\{\sup_{t\in T}p_{i}^{(g)}(t)\right\}<c_{2}<\infty$ that generate $t_{ij}^{(g)}$ according to $t_{ij}^{(g)}=P_{i}^{(g)}{}^{-1}\{j/\left(n_{i}^{(g)}+1\right)\}$ , where $P_{i}^{(g)}{}^{-1}$ is the inverse of $P_{i}^{(g)}(t)=\int_{-\infty}^{t}p_{i}^{(g)}(s)ds$ . For each $g=1,\ldots,G$ , there exist a common sequence of bandwidths $h^{(g)}$ such that $0<c_{1}<$ $\inf_{i}h_{i}^{(g)}/h^{(g)}<\sup_{i}h_{i}^{(g)}/h^{(g)}<c_{2}<\infty$ , where $h_{i}^{(g)}$ is the bandwidth for $\widehat{X}_{i}^{(g)}$ . The kernel density function is smooth and compactly supported.

Let $\mathcal{T}=\left[a_{0},b_{0}\right],t_{i0}^{(g)}=a_{0},t_{i,n_{i}^{(g)}+1}^{(g)}=b_{0}$ , let $\Delta_{i}^{(g)}=\sup\left\{t_{i,j+1}^{(g)}-t_{i,j}^{(g)},j=0,\ldots,n_{i}^{(g)}\right\}$ and $n^{(g)}=n^{(g)}(n)=\inf_{i=1,\ldots,n}n_{i}^{(g)}$ . The condition below is to let the smooth estimate $\widehat{X}_{i}^{(g)}$ serve as well as the true $X_{i}^{(g)}$ in the asymptotic analysis, denoting $0<\lim a_{n}/b_{n}<\infty$ by $a_{n}\sim b_{n}$ .

(B4)

For $g=1,\ldots,G,~{}\sup_{i}\Delta_{i}^{(g)}=O\left(\left(n^{(g)}\right)^{-1}\right),~{}h^{(g)}\sim\left(n^{(g)}\right)^{-1/5},n^{(g)}n^{-5/4}\rightarrow\infty$ .

8 Proof of Theorem

Proof of Lemma 3:

Proof: $\Lambda_{C,K}=\operatorname{diag}\{\lambda_{1}^{(C)},\ldots,\lambda_{K}^{(C)}\}$ consists of first $K$ eigenvalues of $C_{\bm{a}}(\omega)$ . $\Lambda_{B}=\operatorname{diag}\{\lambda_{1}^{(B)},\ldots,\lambda_{K}^{(B)}\}$ consists of first $K$ eigenvalues $BB^{\top}$ and note $i$ -th eigenvalues of $\Lambda_{u}$ as $\lambda_{i}^{(u)}$ .

	$\displaystyle{\bm{f}}_{i}^{\top}(\omega)-\bm{f}_{i}^{\top}=$	$\displaystyle\bm{a}_{i}^{\top}(\bm{\xi}_{1}^{(C)}(\omega),\ldots,\bm{\xi}_{K}^{(C)}(\omega))\Lambda_{C,K}^{-1/2}(\omega)-(\bm{a}_{i}^{\top}-\bm{u}_{i}^{\top})\bm{B}\Lambda_{B}^{-1}$
	$\displaystyle=$	$\displaystyle\bm{a}_{i}^{\top}(\bm{\xi}_{1}^{(C)}(\omega),\ldots,\bm{\xi}_{K}^{(C)}(\omega))\Lambda_{C,K}^{-1/2}(\omega)-(\bm{a}_{i}^{\top}-\bm{u}_{i}^{\top})(\bm{\xi}_{1}^{(B)},\ldots,\bm{\xi}_{K}^{(B)})\Lambda_{B}^{-1/2}$
	$\displaystyle=$	$\displaystyle\bm{a}_{i}^{\top}(\bm{\xi}_{1}^{(B)},\ldots,\bm{\xi}_{K}^{(B)})(\Lambda_{C,K}^{-1/2}(\omega)-\Lambda_{B}^{-1/2})-\bm{a}_{i}^{\top}((\mathbf{B}\mathbf{B}^{\top}-\lambda_{k}^{(B)}I)^{+}\Lambda_{\bm{u}}\bm{\xi}_{k}^{(B)}\omega)_{k=1,\ldots,K}\Lambda_{C,K}^{-1/2}$
		$\displaystyle+\bm{u}_{i}^{\top}\mathbf{B}\Lambda_{B}^{-1}+\bm{a}_{i}^{\top}O(\omega^{2}).$

Then we give the bound separately: For the first part, as $\lambda_{i}^{(C)}=\lambda_{i}^{(B)}+(\bm{\xi}^{(B)}_{k}{}^{\top}\Lambda_{\bm{u}}\bm{\xi}^{(B)}_{k})\omega+O(\omega^{2})$ , we have

	$\displaystyle\frac{1}{\sqrt{\lambda_{i}^{(C)}}}-\frac{1}{\sqrt{\lambda_{i}^{(B)}}}$	$\displaystyle=\frac{\sqrt{\lambda_{i}^{(B)}}-\sqrt{\lambda_{i}^{(C)}}}{\sqrt{\lambda_{i}^{(B)}\lambda_{i}^{(C)}}}$
		$\displaystyle=\frac{{\lambda_{i}^{(B)}}-{\lambda_{i}^{(C)}}}{\sqrt{\lambda_{i}^{(B)}\lambda_{i}^{(C)}}\left(\sqrt{\lambda_{i}^{(B)}}+\sqrt{\lambda_{i}^{(C)}}\right)}$
		$\displaystyle=\frac{-\left((\bm{\xi}^{(B)}_{k}{}^{\top}\Lambda_{\bm{u}}\bm{\xi}^{(B)}_{k})\omega+O(\omega^{2})\right)}{\sqrt{\lambda_{i}^{(B)}\lambda_{i}^{(C)}}\left(\sqrt{\lambda_{i}^{(B)}}+\sqrt{\lambda_{i}^{(C)}}\right)}$

	$\displaystyle\\|\bm{a}_{i}^{\top}(\bm{\xi}_{1}^{(B)},\ldots,\bm{\xi}_{K}^{(B)})(\Lambda_{C,K}^{-1/2}(\omega)-\Lambda_{B}^{-1/2})\\|_{2}^{2}$	$\displaystyle\leq\\|\|\Lambda_{C,K}^{-1/2}(\omega)-\Lambda_{B}^{-1/2}\|\\|_{\mathop{\rm max}}^{2}\\|\bm{a}_{i}\\|_{2}^{2}$
		$\displaystyle\leq C_{s_{1}}^{2}\omega^{2}\\|\bm{a}_{i}\\|_{2}^{2}+O(\omega^{3})\\|\bm{a}_{i}\\|_{2}^{2}$

where $C_{s_{1}}=\frac{\lambda_{\mathop{\rm max}}(\Lambda_{\bm{u}})}{\min_{i\leq K}^{3/2}(\lambda_{i}^{(C)},\lambda_{i}^{(B)})}$ is a constant. For the second part,

	$\displaystyle\|\bm{a}_{i}^{\top}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{k}^{(B)}I)^{+}\Lambda_{\bm{u}}\bm{\xi}_{k}^{(B)}\omega\|^{2}$	$\displaystyle\leq(\frac{\lambda_{\mathop{\rm max}}^{\bm{u}}}{\min_{i<j\leq K}(\lambda_{i}^{(B)}-\lambda_{j}^{(B)})})^{2}\\|\bm{a}_{i}\\|_{2}^{2}\omega^{2}$
	$\displaystyle\\|\bm{a}_{i}^{\top}((\mathbf{B}\mathbf{B}^{\top}-\lambda_{k}^{(B)}I)^{+}\Lambda_{\bm{u}}\bm{\xi}_{k}^{(B)}\omega)_{k=1,\ldots,K}{\Lambda^{C}}^{-1/2}\\|_{2}^{2}$	$\displaystyle\leq\frac{1}{\lambda^{(C)}_{\min}(\omega)}(\frac{\lambda_{\mathop{\rm max}}^{\bm{u}}}{\min_{i<j\leq K}(\lambda_{i}^{(B)}-\lambda_{j}^{(B)})})^{2}\\|\bm{a}_{i}\\|_{2}^{2}\omega^{2}.$

And the third part, with $\operatorname{Cov}(\bm{u}_{i})=\omega\Lambda_{u}$

\displaystyle\|\bm{u}_{i}^{\top}\mathbf{B}\Lambda_{B}^{-1}\|_{2}^{2}\leq\frac{\|\bm{u}_{i}\|_{2}^{2}}{\lambda_{\min}^{(B)}}\leq O_{P}(\omega).

To sum up, we obtain that

\|{\bm{f}}_{i}^{\top}(\omega)-\bm{f}_{i}^{\top}\|_{2}^{2}=O_{P}(\omega^{2}\|\bm{a}_{i}\|_{2}^{2}+\|\bm{u}_{i}^{\top}\mathbf{B}\Lambda_{B}^{-1}\|_{2}^{2})

For $\bm{u}_{i}$ , with

\displaystyle\bm{F}

\displaystyle=(\bm{A}-\bm{U})\mathbf{B}\Lambda^{B}{}^{-1},

that we have

	$\displaystyle\bm{U}$	$\displaystyle=\bm{A}-(\bm{A}-\bm{U})B\Lambda^{B}{}^{-1}\mathbf{B}^{\top}$
		$\displaystyle=\bm{A}-(\bm{A}-\bm{U})\sum_{j=1}^{K}\bm{\xi}_{j}{}^{(B)}\bm{\xi}_{j}^{(B)\top}.$

It means

\displaystyle\bm{U}\sum_{j=K+1}^{p}\bm{\xi}_{j}{}^{(B)}\bm{\xi}_{j}^{(B)\top}=\bm{A}\sum_{j=K+1}^{p}\bm{\xi}_{j}{}^{(B)}\bm{\xi}_{j}^{(B)\top}.

\displaystyle{\bm{u}_{i}}^{\top}(\omega)-\bm{u}_{i}^{\top}=

\displaystyle\bm{a}_{i}^{\top}\sum_{j=K+1}^{p}\bm{\xi}_{j}^{(C)}(\omega)\bm{\xi}_{j}^{(C)}(\omega)^{\top}-\bm{a}_{i}^{\top}\sum_{j=K+1}^{p}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)}{}^{\top}-\bm{u}_{i}^{\top}\sum_{j=1}^{K}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)}{}^{\top}.

Where

		$\displaystyle\bm{a}_{i}^{\top}\sum_{j=K+1}^{p}\bm{\xi}_{j}^{(C)}(\omega)\bm{\xi}_{j}^{(C)}(\omega)^{\top}-\bm{a}_{i}^{\top}\sum_{j=K+1}^{p}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)}{}^{\top}$
	$\displaystyle=$	$\displaystyle\bm{a}_{i}^{\top}\left(\bm{I}_{p}-\sum_{j=1}^{K}\bm{\xi}_{j}^{(C)}(\omega)\bm{\xi}_{j}^{(C)}(\omega)^{\top}\right)-\bm{a}_{i}^{\top}\left(\bm{I}_{p}-\sum_{j=1}^{K}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)}{}^{\top}\right)$
	$\displaystyle=$	$\displaystyle-\bm{a}_{i}^{\top}\sum_{j=1}^{K}\bm{\xi}_{j}^{(C)}(\omega)[\bm{\xi}_{j}^{(C)}(\omega)-\bm{\xi}_{j}^{(B)}]^{\top}+\bm{a}_{i}^{\top}\sum_{j=1}^{K}[\bm{\xi}_{j}^{(B)}-\bm{\xi}_{j}^{(C)}(\omega)]\bm{\xi}_{j}^{(B)}{}^{\top}$
	$\displaystyle=$	$\displaystyle-\bm{a}_{i}^{\top}\sum_{j=1}^{K}[\bm{\xi}_{j}^{(C)}(\omega)-\bm{\xi}_{j}^{(B)}][\bm{\xi}_{j}^{(C)}(\omega)-\bm{\xi}_{j}^{(B)}]^{\top}-\bm{a}_{i}^{\top}\sum_{j=1}^{K}\bm{\xi}_{j}^{(B)}[\bm{\xi}_{j}^{(C)}(\omega)-\bm{\xi}_{j}^{(B)}]^{\top}$
		$\displaystyle+\bm{a}_{i}^{\top}\sum_{j=1}^{K}[\bm{\xi}_{j}^{(B)}-\bm{\xi}_{j}^{(C)}(\omega)]\bm{\xi}_{j}^{(B)}{}^{\top}$
	$\displaystyle=$	$\displaystyle-\omega^{2}\bm{a}_{i}^{\top}\sum_{j=1}^{K}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\Lambda_{u}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)\top}\Lambda_{u}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}+\bm{a}_{i}^{\top}O(\omega^{3})$
		$\displaystyle+\omega\bm{a}_{i}^{\top}\sum_{j=1}^{K}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)}{}^{\top}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\Lambda_{u}+\omega\bm{a}_{i}^{\top}\sum_{j=1}^{K}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\Lambda_{u}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)\top}$
		$\displaystyle+\bm{a}_{i}^{\top}O(\omega^{2}).$

Then we give the bound separately:

	$\displaystyle\\|\bm{a}_{i}^{\top}\sum_{j=1}^{K}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\Lambda_{u}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)\top}\Lambda_{u}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\\|_{2}^{2}$
	$\displaystyle\leq\\|\sum_{j=1}^{K}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\Lambda_{u}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)\top}\Lambda_{u}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\\|_{2}^{2}\\|\bm{a}_{i}\\|_{2}^{2}$
	$\displaystyle\leq(\lambda^{u}_{max})^{4}\\|\sum_{j=1}^{K}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)\top}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\\|_{2}^{2}\\|\bm{a}_{i}\\|_{2}^{2}$
	$\displaystyle\leq(p-K)(\frac{\lambda_{\mathop{\rm max}}^{\bm{u}}}{\min_{i<j\leq K}(\lambda_{i}^{(B)}-\lambda_{j}^{(B)})})^{4}\\|\bm{a}_{i}\\|_{2}^{2}.$

\displaystyle\|\bm{a}_{i}^{\top}\sum_{j=1}^{K}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)}{}^{\top}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\Lambda_{u}\|_{2}^{2}\leq(m_{1}+m_{2}-r)(\frac{\lambda_{\mathop{\rm max}}^{\bm{u}}}{\min_{i<j\leq K}(\lambda_{i}^{(B)}-\lambda_{j}^{(B)})})^{2}\|\bm{a}_{i}\|_{2}^{2}

To sum up, remove the constant then we can obtain that

\|{\bm{u}}_{i}^{\top}(\omega)-\bm{u}_{i}^{\top}\|_{2}^{2}=O_{P}(\omega^{2}\|\bm{a}_{i}\|_{2}^{2}+\|\bm{u}_{i}^{\top}\mathbf{B}\Lambda_{B}^{-1}\mathbf{B}^{\top}\|_{2}^{2})

Proof of Lemma 4:

Proof:

Under (A1,A2) (B1-B4), using the Lemma 2(a) in Kong et al. (2016). $i=1,\ldots,n;k_{1},k_{2}=1,\ldots,s_{n};j=1,\ldots,G$

	$\displaystyle\hat{\bm{A}}^{\top}\hat{\bm{A}}$	$\displaystyle=(\sum_{i=1}^{n}\bm{a}_{i}\bm{a}_{i}^{\top})_{Gs_{n}\times Gs_{n}}=\sum_{i=1}^{n}\left(\begin{array}[]{ccc}\bm{a}_{i}^{(1)}\bm{a}_{i}^{(1)\top}&\ldots&\bm{a}_{i}^{(1)}\bm{a}_{i}^{(G)\top}\\ \ldots&\ldots&\ldots\\ \bm{a}_{i}^{(G)}\bm{a}_{i}^{(1)\top}&\ldots&\bm{a}_{i}^{(G)}\bm{a}_{i}^{(G)\top}\end{array}\right)$
		$\displaystyle=\sum_{i=1}^{n}\left(\begin{array}[]{ccc}(a_{ik_{1}}^{(1)}a_{ik_{2}}^{(1)})_{s_{n}\times s_{n}}&\ldots&(a_{ik_{1}}^{(1)}a_{ik_{2}}^{(G)})_{s_{n}\times s_{n}}\\ \ldots&\ldots&\ldots\\ (a_{ik_{1}}^{(G)}a_{ik_{2}}^{(1)})_{s_{n}\times s_{n}}&\ldots&(a_{ik_{1}}^{(G)}a_{ik_{2}}^{(G)})_{s_{n}\times s_{n}}\end{array}\right)$
		$\displaystyle=\left(\begin{array}[]{ccc}(\sum_{i=1}^{n}a_{ik_{1}}^{(1)}a_{ik_{2}}^{(1)})_{s_{n}\times s_{n}}&\ldots&(\sum_{i=1}^{n}a_{ik_{1}}^{(1)}a_{ik_{2}}^{(G)})_{s_{n}\times s_{n}}\\ \ldots&\ldots&\ldots\\ (\sum_{i=1}^{n}a_{ik_{1}}^{(G)}a_{ik_{2}}^{(1)})_{s_{n}\times s_{n}}&\ldots&(\sum_{i=1}^{n}a_{ik_{1}}^{(G)}a_{ik_{2}}^{(G)})_{s_{n}\times s_{n}}\end{array}\right),$

then the analysis for $\{(\bm{A}^{\top}\bm{A})_{ij}-E(\bm{A}^{\top}\bm{A})_{ij}\}/n$ is is equivalent to analyzing the following .

	$\displaystyle n^{-1}\sum_{i=1}^{n}[\hat{a}_{ik_{1}}^{(j)}\hat{a}_{ik_{2}}^{(j)}-E(a_{ik_{1}}^{(j)}a_{ik_{2}}^{(j)})]$	$\displaystyle=n^{-1}\sum_{i=1}^{n}[(\hat{a}_{ik_{1}}^{(j)}\hat{a}_{ik_{2}}^{(j)}-a_{ik_{1}}^{(j)}a_{ik_{2}}^{(j)})+(a_{ik_{1}}^{(j)}a_{ik_{2}}^{(j)}-E(a_{ik_{1}}^{(j)}a_{ik_{2}}^{(j)}))]$
		$\displaystyle=O_{p}(k_{1}^{a/2+1}n^{-1/2}+k_{2}^{a/2+1}n^{-1/2}+n^{-1/2})$
		$\displaystyle\leq O_{p}(s_{n}^{a/2+1}n^{-1/2})$

Proof of Theorem 1

Proof: We view the estimation error as perturbation

\displaystyle\frac{1}{n}\bm{A}^{\top}\bm{A}=\mathbb{E}(\bm{a}_{i}\bm{a}_{i}^{\top})+\frac{1}{\sqrt{n}}\Lambda_{e}+O(\frac{1}{n})

where $\Lambda_{e}/\sqrt{n}+O(1/n)$ represents the estimation error, and the equation holds true with probability. $\Lambda_{e}$ is a bounded random matrix since $\bm{A}$ is bounded. Consider

	$\displaystyle\widehat{\bm{f}}_{i}^{\top}-\bm{f}_{i}^{\top}(\omega)=$	$\displaystyle\widehat{\bm{a}}_{i}^{\top}(\widehat{\bm{\xi}}_{1},\ldots,\widehat{\bm{\xi}}_{K}){\widehat{\Lambda}_{K}}^{-1/2}-\bm{a}_{i}^{\top}(\bm{\xi}_{1}^{(C)}(\omega),\ldots,\bm{\xi}_{K}^{(C)}(\omega))\Lambda_{C,K}^{-1/2}(\omega),$
	$\displaystyle\widehat{\bm{u}}_{i}^{\top}-\bm{u}_{i}^{\top}(\omega)=$	$\displaystyle\widehat{\bm{a}}_{i}^{\top}\sum_{j=K+1}^{p}\widehat{\bm{\xi}}_{j}\widehat{\bm{\xi}}_{j}^{\top}-\bm{a}_{i}^{\top}\sum_{j=K+1}^{p}\bm{\xi}_{j}^{(C)}(\omega)\bm{\xi}_{j}^{(C)}(\omega){}^{\top}.$

Follow the steps in proof of Lemma 3, we have

	$\displaystyle\\|\widehat{\bm{f}}_{i}^{\top}-\bm{f}_{i}^{\top}(\omega)\\|_{2}^{2}$	$\displaystyle\leq O_{P}(\frac{1}{n}\\|\bm{a}_{i}\\|_{2}^{2})$
	$\displaystyle\\|\widehat{\bm{u}}_{i}^{\top}-\bm{u}_{i}^{\top}(\omega)\\|_{2}^{2}$	$\displaystyle\leq O_{P}(\frac{1}{n}\\|\bm{a}_{i}\\|_{2}^{2})$

Proof of Theorem 2

Proof: Recall that

\widehat{\boldsymbol{\theta}}=\left(\begin{array}[]{c}\widehat{\mathbf{H}}_{1}\\ \widehat{\mathbf{\gamma}}\end{array}\right)=\operatorname{argmin}_{\boldsymbol{\theta}}\left\{L_{n}(\mathbf{y},\widehat{\mathbf{W}}\boldsymbol{\theta})+\lambda\left\|\boldsymbol{\theta}_{[p_{1}]}\right\|_{1}\right\}.

(23)

Also, Assumption 5 tells us $\mathbf{V}_{0}$ is nonsingular and so is $\mathbf{V}=\left(\begin{array}[]{cc}\mathbf{I}_{p_{1}}&\mathbf{0}_{p_{1}\times K}\\ \mathbf{0}_{K\times p_{1}}&\mathbf{V}_{0}\end{array}\right)$ .

Define $\overline{\mathbf{W}}=\left(1,\bar{\bm{u}}_{t}^{\top},\bar{\bm{f}}_{t}^{\top}\right)=\widehat{\mathbf{W}}\mathbf{V}=\left(1,\widehat{\bm{u}}_{t}^{\top},\widehat{\bm{f}}_{t}^{\top}\right)\cdot\mathbf{V}$ , $\bar{\boldsymbol{\theta}}=\mathbf{V}^{-1}\widehat{\boldsymbol{\theta}}$ , $\widehat{\mathbf{B}}_{0}=\left(\mathbf{0}_{K},\widehat{\mathbf{B}}^{\top}\right)^{\top},\widehat{\boldsymbol{\theta}}^{*}=\left(\begin{array}[]{c}\mathbf{H}_{1}^{*}\\ \widehat{\mathbf{B}}_{0}^{\top}\mathbf{H}_{1}^{*}\end{array}\right)$ and $\bar{\boldsymbol{\theta}}^{*}=\mathbf{V}^{-1}\widehat{\boldsymbol{\theta}}^{*}$ . We easily see that $\widehat{\mathbf{H}}=\widehat{\boldsymbol{\theta}}_{[p_{1}]}=\bar{\boldsymbol{\theta}}_{[p_{1}]}$ and

\bar{\boldsymbol{\theta}}=\operatorname{argmin}_{\boldsymbol{\theta}}\left\{L_{n}(\mathbf{y},\overline{\mathbf{W}}\boldsymbol{\theta})+\lambda\left\|\boldsymbol{\theta}_{[p_{1}]}\right\|_{1}\right\}.

Then it follows that $\operatorname{supp}(\widehat{\mathbf{H}}_{1})=\operatorname{supp}\left(\bar{\boldsymbol{\theta}}_{[p_{1}]}\right)$ and $\left\|\widehat{\mathbf{H}}_{1}-\mathbf{H}_{1}{}^{*}\right\|=\left\|\bar{\boldsymbol{\theta}}_{[p_{1}]}-\bar{\boldsymbol{\theta}}_{[p_{1}]}^{*}\right\|\leq\left\|\bar{\boldsymbol{\theta}}-\bar{\boldsymbol{\theta}}^{*}\right\|$ for any norm $\|\cdot\|$ .

Consequently, Theorem 2 is reduced to studying $\bar{\boldsymbol{\theta}}$ and the loss function $L_{n}(\mathbf{y},\overline{\mathbf{W}}\boldsymbol{\theta})$ . The Lemma 5 below implies that all the regularity conditions (with $A=\infty$ ) in Lemma 6 are satisfied.

Let $w_{tj}$ and $\bar{w}_{tj}$ be the $(t,j)$ -th element of $\mathbf{W}$ and $\overline{\mathbf{W}}$ , respectively. Observe that $L_{n}(\mathbf{y},\overline{\mathbf{W}}\boldsymbol{\theta})=$ $\frac{1}{n}\sum_{t=1}^{n}\left[-y_{t}\overline{\mathbf{w}}_{t}^{\top}\boldsymbol{\theta}+b\left(\overline{\mathbf{w}}_{t}^{\top}\boldsymbol{\theta}\right)\right],\nabla L_{n}(\mathbf{y},\overline{\mathbf{W}}\boldsymbol{\theta})=\frac{1}{n}\sum_{t=1}^{n}\left[-y_{t}+b^{\prime}\left(\overline{\mathbf{w}}_{t}^{T}\boldsymbol{\theta}\right)\right]\overline{\mathbf{w}}_{t}$ and $\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}=\mathbf{A}\boldsymbol{H}^{*}+(\beta_{0}^{*})_{n\times 1}$ . Hence $\left\|\nabla L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}\right)\right\|_{\infty}=\varepsilon^{*}$ and consequently, $\left\|\nabla_{S}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}\right)\right\|_{\infty}\leq\varepsilon^{*},\left\|\nabla_{S}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}\right)\right\|_{2}\leq\varepsilon^{*}\sqrt{|S|}$ and $\left\|\nabla_{S}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}\right)\right\|_{1}\leq\varepsilon^{*}|S|$ . In addition, $\lambda>7\varepsilon^{*}/\tau\geq\varepsilon^{*}$ .

Based on these estimates, all the results follow from Theorem 2 and some simple algebra.

Here we present the Lemma 5 used above and its proof.

Lemma 5

(Lemma C.1 of Fan et al. (2020)) Let Assumptions 3, 4 and 5 hold. Treat $L_{n}(\mathbf{y},\overline{\mathbf{W}}\boldsymbol{\theta})$ as a function of $\boldsymbol{\theta}$ , and the derivatives below are taken with respect to it. Firstly, by Assumption 5, we have $\|\bm{W}\|_{\mathop{\rm max}}\leq\frac{M_{0}}{2}$ . Then

(i)

$\left\|\nabla_{.S}^{2}L(\mathbf{y},\overline{\mathbf{W}}\boldsymbol{\theta})-\nabla_{.S}^{2}L\left(\mathbf{y},\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}\right)\right\|_{\infty}\leq M\left\|\boldsymbol{\theta}-\bar{\boldsymbol{\theta}}^{*}\right\|_{2},\forall\bm{\theta}$ with $\operatorname{supp}(\boldsymbol{\theta})\subseteq S$ ,
(ii)

$\left\|\left(\nabla_{SS}^{2}L\left(\mathbf{y},\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}\right)\right)^{-1}\right\|_{\infty}\leq\frac{1}{2\kappa_{\infty}}$ ,
(iii)

$\left\|\left(\nabla_{SS}^{2}L\left(\mathbf{y},\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}\right)\right)^{-1}\right\|_{2}\leq\frac{1}{2\kappa_{2}}$ ,
(iv)

$\left\|\nabla_{S_{2}S}^{2}L\left(\mathbf{y},\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}\right)\left(\nabla_{SS}^{2}L\left(\mathbf{y},\overline{\mathbf{W}}\bar{\boldsymbol{\theta}}^{*}\right)\right)^{-1}\right\|_{\infty}\leq 1-\tau$ .

Proof of Lemma 5 (i) Based on the fact that $\mathbf{W}\boldsymbol{\theta}^{*}=\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}={\bm{A}}\bm{H}^{*}+(\beta_{0}^{*})_{n\times 1}$ , we have $\nabla^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)=\frac{1}{n}\sum_{t=1}^{n}b^{\prime\prime}\left(\overline{\mathbf{w}}_{t}^{T}\overline{\boldsymbol{\theta}}^{*}\right)\mathbf{w}_{t}\mathbf{w}_{t}^{T}$ and $\nabla^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)=\frac{1}{n}\sum_{t=1}^{n}b^{\prime\prime}\left(\overline{\mathbf{w}}_{t}^{T}\overline{\boldsymbol{\theta}}^{*}\right)\overline{\mathbf{w}}_{t}\overline{\mathbf{w}}_{t}^{T}$ . For any $j,k\in[p_{1}+K]$ and $\operatorname{supp}(\boldsymbol{\theta})\subseteq S$ ,

		$\displaystyle\left\|\nabla_{jk}^{2}L_{n}(\mathbf{y},\overline{\mathbf{W}}\boldsymbol{\theta})-\nabla_{jk}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\right\|\leq\frac{1}{n}\sum_{t=1}^{n}\left\|b^{\prime\prime}\left(\overline{\mathbf{w}}_{t}^{T}\boldsymbol{\theta}\right)-b^{\prime\prime}\left(\overline{\mathbf{w}}_{t}^{T}\overline{\boldsymbol{\theta}}^{}\right)\right\|\cdot\left\|\overline{\mathbf{w}}_{tj}\overline{\mathbf{w}}_{tk}\right\|$
		$\displaystyle\leq\frac{1}{n}\sum_{t=1}^{n}M_{3}\left\|\overline{\mathbf{w}}_{t}^{T}\left(\boldsymbol{\theta}-\overline{\boldsymbol{\theta}}^{*}\right)\right\|\cdot\\|\overline{\mathbf{W}}\\|_{\mathop{\rm max}}^{2}$		(24)

By the Cauchy-Schwarz inequality and $\|\overline{\mathbf{W}}\|_{\mathop{\rm max}}\leq\|\mathbf{W}\|_{\mathop{\rm max}}+\|\overline{\mathbf{W}}-\mathbf{W}\|_{\mathop{\rm max}}\leq M_{0}$ , we obtain that for $i\in[n],\left|\overline{\mathbf{w}}_{t}^{T}\left(\boldsymbol{\theta}-\overline{\boldsymbol{\theta}}^{*}\right)\right|=\left|\overline{\mathbf{w}}_{tS}^{T}\left(\boldsymbol{\theta}-\overline{\boldsymbol{\theta}}^{*}\right)_{S}\right|\leq\left\|\overline{\mathbf{w}}_{tS}\right\|_{2}\left\|\boldsymbol{\theta}-\overline{\boldsymbol{\theta}}^{*}\right\|_{2}\leq\sqrt{|S|}M_{0}\left\|\boldsymbol{\theta}-\overline{\boldsymbol{\theta}}^{*}\right\|_{2}$ . Plugging this result back to (8), we get

		$\displaystyle\left\|\nabla_{jk}^{2}L_{n}(\mathbf{y},\overline{\mathbf{W}}\boldsymbol{\theta})-\nabla_{jk}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\right\|\leq\sqrt{\|S\|}M_{3}M_{0}^{3}\left\\|\boldsymbol{\theta}-\overline{\boldsymbol{\theta}}^{}\right\\|_{2},\quad\forall j,k\in[p+K]$
		$\displaystyle\left\\|\nabla_{\cdot S}^{2}L_{n}(\mathbf{y},\overline{\mathbf{W}}\boldsymbol{\theta})-\nabla_{\cdot S}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\right\\|_{\infty}\leq\|S\|^{3/2}M_{3}M_{0}^{3}\left\\|\boldsymbol{\theta}-\overline{\boldsymbol{\theta}}^{}\right\\|_{2}=M\left\\|\boldsymbol{\theta}-\overline{\boldsymbol{\theta}}^{*}\right\\|_{2}$

(ii) Now we come to the second claim. For any $k\in[p+K]$ ,

		$\displaystyle\left\\|\nabla_{kS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)-\nabla_{kS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\right\\|_{\infty}\leq\frac{1}{n}\sum_{t=1}^{n}b^{\prime\prime}\left(\mathbf{x}_{t}^{T}\boldsymbol{\beta}^{*}\right)\left\\|\bar{w}_{tk}\overline{\mathbf{w}}_{tS}^{T}-w_{tk}\mathbf{w}_{tS}^{T}\right\\|_{\infty}$
		$\displaystyle\leq\frac{M_{2}\sqrt{\|S\|}}{n}\sum_{t=1}^{n}\left\\|\bar{w}_{tk}\overline{\mathbf{w}}_{tS}^{T}-w_{tk}\mathbf{w}_{tS}^{T}\right\\|_{2}$

Also, by $\|\mathbf{W}\|_{\mathop{\rm max}}\leq M_{0}/2$ and $\|\overline{\mathbf{W}}\|_{\mathop{\rm max}}\leq M_{0}$ we have

		$\displaystyle\left\\|\bar{w}_{tk}\overline{\mathbf{w}}_{tS}^{T}-w_{tk}\mathbf{w}_{tS}^{T}\right\\|_{2}\leq\left\|w_{tk}\right\|\cdot\left\\|\left(\overline{\mathbf{w}}_{tS}-\mathbf{w}_{tS}\right)^{T}\right\\|_{2}+\left\|\bar{w}_{tk}-w_{tk}\right\|\cdot\left\\|\overline{\mathbf{w}}_{tS}^{T}\right\\|_{2}$
		$\displaystyle\leq\\|\mathbf{W}\\|_{\mathop{\rm max}}\left\\|\overline{\mathbf{w}}_{tS}-\mathbf{w}_{tS}\right\\|_{2}+\left\|\bar{w}_{tk}-w_{tk}\right\|\cdot\sqrt{\|S\|}\\|\overline{\mathbf{W}}\\|_{\mathop{\rm max}}$
		$\displaystyle\leq\frac{M_{0}}{2}\left\\|\overline{\mathbf{w}}_{tS}-\mathbf{w}_{tS}\right\\|_{2}+M_{0}\sqrt{\|S\|}\cdot\left\|\bar{w}_{tk}-w_{tk}\right\|$

Define $\delta=\mathop{\rm max}_{j\in[p+K]}\left(\frac{1}{n}\sum_{t=1}^{n}\left|\bar{w}_{tj}-w_{tj}\right|^{2}\right)^{1/2}$ . By the Jensen’s inequality, $\forall J\subseteq[p+K]$ ,

\frac{1}{n}\sum_{t=1}^{n}\left\|\overline{\mathbf{w}}_{tJ}-\mathbf{w}_{tJ}\right\|_{2}\leq\left(\frac{1}{n}\sum_{t=1}^{n}\left\|\overline{\mathbf{w}}_{tJ}-\mathbf{w}_{tJ}\right\|_{2}^{2}\right)^{1/2}\leq\left(\frac{|J|}{n}\mathop{\rm max}_{j\in[p+K]}\sum_{t=1}^{n}\left|\bar{w}_{tj}-w_{tj}\right|^{2}\right)^{1/2}\leq\sqrt{|J|}\delta.

As a result,

		$\displaystyle\left\\|\nabla_{\cdot S}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)-\nabla_{\cdot S}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\right\\|_{\infty}$
		$\displaystyle=\mathop{\rm max}_{k\in[p+K]}\left\\|\nabla_{kS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)-\nabla_{kS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\right\\|_{\infty}\leq\frac{3}{2}M_{0}M_{2}\|S\|\delta$		(25)

Let $\alpha=\left\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right)^{-1}\left[\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)-\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right]\right\|_{\infty}$ . Then

	$\displaystyle\alpha$	$\displaystyle\leq\left\\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\right)^{-1}\right\\|_{\infty}\left\\|\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}^{T}\overline{\boldsymbol{\theta}}^{}\right)-\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right\\|_{\infty}$
		$\displaystyle\leq\frac{3}{8\kappa_{\infty}}M_{0}M_{2}\|S\|\delta\leq\frac{1}{2}.$		(26)

Lemma 7 yields

		$\displaystyle\left\\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\right)^{-1}-\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\right)^{-1}\right\\|_{\infty}\leq\left\\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right)^{-1}\right\\|_{\infty}\frac{\alpha}{1-\alpha}$
		$\displaystyle\leq\frac{1}{4\kappa_{\infty}}\cdot\frac{\alpha}{1-\frac{1}{2}}\leq\frac{3}{16\kappa_{\infty}^{2}}M_{0}M_{2}\|S\|\delta.$

We also have a cruder bound $\left\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)\right)^{-1}-\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right)^{-1}\right\|_{\infty}\leq\frac{1}{4\kappa_{\infty}}$ , which leads to

\displaystyle\left\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)\right)^{-1}\right\|_{\infty}\leq\left\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right)^{-1}\right\|_{\infty}+\frac{1}{4\kappa_{\infty}}\leq\frac{1}{2\kappa_{\infty}}

(27)

(iii) The third argument follows (27) easily. Since $\|\mathbf{A}\|_{2}\leq\|\mathbf{A}\|_{\infty}$ holds for any symmetric matrix $\mathbf{A}$ , we have $\left\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)\right)^{-1}-\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right)^{-1}\right\|_{2}\leq\frac{1}{4\kappa_{\infty}}\leq\frac{1}{4\kappa_{2}}$ and thus $\left\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)\right)^{-1}\right\|_{2}\leq\frac{1}{2\kappa_{2}}$ . (iv) Finally we prove the last inequality. On the one hand,

		$\displaystyle\left\\|\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\right)^{-1}-\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\right)^{-1}\right\\|_{\infty}$
		$\displaystyle\leq\left\\|\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)-\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\right\\|_{\infty}\left\\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)\right)^{-1}\right\\|_{\infty}$
		$\displaystyle+\left\\|\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\left[\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\right)^{-1}-\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right)^{-1}\right]\right\\|_{\infty}$

From claim (ii) and (8) it is easy to see that

\left\|\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)-\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right\|_{\infty}\left\|\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)\right)^{-1}\right\|_{\infty}\leq\frac{1}{4\kappa_{\infty}}3M_{0}M_{2}|S|\delta.

On the other hand, we can take $\mathbf{A}=\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right),\mathbf{B}=\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)$ and $\mathbf{C}=\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)-$ $\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)$ . By Assumption 4, $\left\|\mathbf{AB}^{-1}\right\|_{\infty}\leq 1-2\tau\leq 1$ . Lemma 7 forces that

		$\displaystyle\left\\|\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\left[\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\right)^{-1}-\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right)^{-1}\right]\right\\|_{\infty}$
		$\displaystyle=\left\\|\mathbf{A}\left[(\mathbf{B}+\mathbf{C})^{-1}-\mathbf{B}^{-1}\right]\right\\|_{\infty}\leq\left\\|\mathbf{A}B^{-1}\right\\|_{\infty}\frac{\left\\|\mathbf{CB}^{-1}\right\\|_{\infty}}{1-\left\\|\mathbf{CB}^{-1}\right\\|_{\infty}}\leq\frac{\\|\mathbf{C}\\|_{\infty}\left\\|\mathbf{B}^{-1}\right\\|_{\infty}}{1-\\|\mathbf{C}\\|_{\infty}\left\\|\mathbf{B}^{-1}\right\\|_{\infty}}$

We have shown above in (8) that $\|\mathbf{C}\|_{\infty}\left\|\mathbf{B}^{-1}\right\|_{\infty}\leq\frac{3}{8\kappa_{\infty}}M_{0}M_{2}|S|\delta\leq 1/2$ . As a result,

\left\|\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\left[\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)\right)^{-1}-\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{*}\right)\right)^{-1}\right]\right\|_{\infty}\leq\frac{3}{4\kappa_{\infty}}M_{0}M_{2}|S|\delta.

By combining these estimates, we have

		$\displaystyle\left\\|\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{}\right)\right)^{-1}-\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\mathbf{W}\boldsymbol{\theta}^{}\right)\right)^{-1}\right\\|_{\infty}$
		$\displaystyle\leq\frac{3}{2\kappa_{\infty}}M_{0}M_{2}\|S\|\delta\leq\tau.$

Therefore $\left\|\nabla_{S_{2}S}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)\left(\nabla_{SS}^{2}L_{n}\left(\mathbf{y},\overline{\mathbf{W}}\overline{\boldsymbol{\theta}}^{*}\right)\right)^{-1}\right\|_{\infty}\leq(1-2\tau)+\tau=1-\tau$ .

Finally, we use the following lemma to prove the theorem.

Lemma 6

(Theorem B.1 of Fan et al. (2020)) For optimization 23, if the Lemma 5 is statisfied, then we have that

Error bounds: if

\frac{7}{\tau}\left\|\nabla L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\|_{\infty}<\lambda<\frac{\kappa_{2}}{4\sqrt{|S|}}\min\left\{A,\frac{\kappa_{\infty}\tau}{3M}\right\},

then $\operatorname{supp}(\widehat{\boldsymbol{\theta}})\subseteq S$ and

		$\displaystyle\left\\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{\infty}\leq\frac{3}{5\kappa_{\infty}}\left(\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}+\lambda\right),$
		$\displaystyle\left\\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{2}\leq\frac{2}{\kappa_{2}}\left(\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{2}+\lambda\sqrt{\left\|S_{1}\right\|}\right),$
		$\displaystyle\left\\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{1}\leq\min\left\{\frac{3}{5\kappa_{\infty}}\left(\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{1}+\lambda\left\|S_{1}\right\|\right),\frac{2\sqrt{\|S\|}}{\kappa_{2}}\left(\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\\|_{2}+\lambda\sqrt{\left\|S_{1}\right\|}\right)\right\}.$

Sign consistency : In addition, if the following two conditions

		$\displaystyle\min\left\{\left\|{\mathbf{H}_{1}^{}}_{j}\right\|:{\mathbf{H}_{1}^{}}_{j}\neq 0,j\in[p_{1}]\right\}>\frac{C}{\kappa_{\infty}\tau}\left\\|\nabla L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\\|_{\infty}$
		$\displaystyle\left\\|\nabla L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\\|_{\infty}<\frac{\kappa_{2}\tau}{7C\sqrt{\|S\|}}\min\left\{A,\frac{\kappa_{\infty}\tau}{3M}\right\}$

hold for some $C\geq 5$ , then by taking $\lambda\in\left(\frac{7}{\tau}\left\|\nabla L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\|_{\infty},\frac{1}{\tau}\left(\frac{5C}{3}-1\right)\left\|\nabla L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\|_{\infty}\right)$ , the estimator achieves the sign consistency $\operatorname{sign}(\widehat{\mathbf{H}})=\operatorname{sign}\left(\mathbf{H}^{*}\right)$ .

Proof of Lemma 6

We need the following two lemma to prove the Lemma 6.

Define $B_{S}\left(\boldsymbol{\theta}^{*},r\right)=\left\{\boldsymbol{\theta}:\left\|\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\right\|_{2}\leq r,\operatorname{supp}(\boldsymbol{\theta})\subseteq S\right\}$ for $r>0$ . We first introduce two useful lemmas.

Lemma 7

Suppose $\mathbf{A}\in\mathbb{R}^{q\times r}$ and $\mathbf{B},\mathbf{C}\in\mathbb{R}^{r\times r}$ and $\left\|\mathbf{CB}^{-1}\right\|<1$ , where $\|\cdot\|$ is an induced norm. Then $\left\|\mathbf{A}\left[(\mathbf{B}+\mathbf{C})^{-1}-\mathbf{B}^{-1}\right]\right\|\leq\frac{\left\|\mathbf{AB}^{-1}\right\|\cdot\left\|\mathbf{CB}^{-1}\right\|}{1-\left\|\mathbf{CB}^{-1}\right\|}$ .

Lemma 8

Under (i) (ii) (iii) in Lemma 5, we have $\left\|\left(\nabla_{SS}^{2}L_{n}(\boldsymbol{\theta})\right)^{-1}\right\|_{2}<\kappa_{2}^{-1}$ and $\left\|\left(\nabla_{SS}^{2}L_{n}(\boldsymbol{\theta})\right)^{-1}\right\|_{\infty}$ $<\kappa_{\infty}^{-1}$ over $B_{S}\left(\boldsymbol{\theta}^{*},\min\left\{A,\frac{\kappa_{\infty}}{M}\right\}\right)$ .

Lemma 9

Suppose $\lambda\geq 0,\mathcal{M}$ is a Euclidean space, $\boldsymbol{\theta}_{0}\in\mathcal{M},L(\boldsymbol{\theta})\in C^{2}(\mathcal{M})$ and is convex, and $R(\boldsymbol{\theta})$ is convex. In addition, there exist $\kappa,A>0$ such that $\nabla^{2}L(\boldsymbol{\theta})\succeq\kappa I$ as long as $\left\|\boldsymbol{\theta}-\boldsymbol{\theta}_{0}\right\|_{2}\leq$ A. If $\left\|\nabla L\left(\boldsymbol{\theta}_{0}\right)\right\|_{2}+\lambda\inf_{\mathbf{h}\in\partial R\left(\boldsymbol{\theta}_{0}\right)}\|\mathbf{h}\|_{2}<\frac{1}{2}\kappa A$ , then $L_{\lambda}(\boldsymbol{\theta})=L(\boldsymbol{\theta})+\lambda R(\boldsymbol{\theta})$ has unique minimizer $\hat{\boldsymbol{\theta}}$ and $\left\|\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\right\|_{2}\leq\frac{2}{\kappa}\left(\left\|\nabla L\left(\boldsymbol{\theta}_{0}\right)\right\|_{2}+\lambda\inf_{\mathbf{h}\in\partial R\left(\boldsymbol{\theta}_{0}\right)}\|\mathbf{h}\|_{2}\right)$ .

Lemma 10

Suppose $L(\boldsymbol{\theta})\in C^{2}\left(\mathbb{R}^{p}\right)$ and is convex. $R(\boldsymbol{\theta})$ is convex and $R(\boldsymbol{\alpha}+\boldsymbol{\beta})=R(\boldsymbol{\alpha})+$ $R(\boldsymbol{\beta})$ for $\boldsymbol{\alpha}\in\mathcal{M}$ and $\boldsymbol{\beta}\in\mathcal{M}^{\perp}$ , where $\mathcal{M}$ is a linear subspace of $\mathbb{R}^{p}$ and $\mathcal{M}^{\perp}$ is its orthonormal complement. In addition, there exists $R^{*}(\boldsymbol{\theta})\in C\left(\mathbb{R}^{p}\right)$ such that $|\langle\boldsymbol{\alpha},\boldsymbol{\beta}\rangle|\leq R(\boldsymbol{\alpha})R^{*}(\boldsymbol{\beta})$ for $\boldsymbol{\alpha}\in\mathcal{M}^{\perp}$ and $\boldsymbol{\beta}\in\mathbb{R}^{p}$ . Let $L_{\lambda}(\boldsymbol{\theta})=L(\boldsymbol{\theta})+\lambda R(\boldsymbol{\theta})$ where $\lambda\geq 0$ , and $\hat{\boldsymbol{\theta}}\in\operatorname{argmin}_{\boldsymbol{\theta}\in\mathcal{M}}L_{\lambda}(\boldsymbol{\theta})$ .

If $R^{*}(\nabla L(\hat{\boldsymbol{\theta}}))<\lambda$ and $\boldsymbol{\theta}^{T}\nabla^{2}L(\hat{\boldsymbol{\theta}})\boldsymbol{\theta}>0$ for all $\boldsymbol{\theta}\in\mathcal{M}$ , then $\hat{\boldsymbol{\theta}}$ is the unique global minimizer of $L_{\lambda}(\boldsymbol{\theta})$ .

Now we start the proof: First we study the restricted problem

\bar{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}\in\mathcal{M}}\left\{L_{n}(\boldsymbol{\theta})+\lambda R(\boldsymbol{\theta})\right\},

where $\mathcal{M}=\{\boldsymbol{\theta}:\boldsymbol{\theta}\in\mathbb{R}^{p+K},\boldsymbol{\theta}_{S_{2}}=\boldsymbol{0}\}$ is oracle parameter set. Take $R(\boldsymbol{\theta})=\left\|\boldsymbol{\theta}_{[p]}\right\|_{1}$ and $R^{*}(\boldsymbol{\theta})=\left\|\boldsymbol{\theta}_{S_{2}}\right\|_{\infty}$ . Let $A_{1}=\min\left\{A,\frac{\kappa_{\infty}\tau}{3M}\right\}$ and hence $A_{1}\leq\min\left\{A,\frac{\kappa_{\infty}}{M}\right\}$ . Lemma 8 shows that $\left\|\left(\nabla_{SS}^{2}L_{n}(\boldsymbol{\theta})\right)^{-1}\right\|_{2}<\kappa_{2}^{-1}$ and $\left\|\left(\nabla_{SS}^{2}L_{n}(\boldsymbol{\theta})\right)^{-1}\right\|_{\infty}<\kappa_{\infty}^{-1}$ over $B_{S}\left(\boldsymbol{\theta}^{*},A_{1}\right)$ .

Since $\operatorname{supp}\left(\boldsymbol{\theta}^{*}\right)\subseteq S$ , any $\mathbf{h}\in\partial R\left(\boldsymbol{\theta}^{*}\right)$ satisfies $\|\mathbf{h}\|_{2}\leq\sqrt{\left|S_{1}\right|}$ . Therefore

\left\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\|_{2}+\lambda\|\mathbf{h}\|_{2}\leq\frac{1}{2}\kappa_{2}A_{1}\leq\frac{1}{2}\kappa_{2}A.

Then Lemma 9 implies that $\left\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\|_{2}\leq\frac{2}{\kappa_{2}}\left(\left\|\nabla_{S}L\left(\boldsymbol{\theta}^{*}\right)\right\|_{2}+\lambda\sqrt{\left|S_{1}\right|}\right)\leq A_{1}$ .

Second, we study the $L_{\infty}$ bound. On the one hand, the optimality condition yields $\nabla_{S}L_{n}(\bar{\boldsymbol{\theta}})\in-\lambda\partial\left\|\bar{\boldsymbol{\theta}}_{S}\right\|_{\infty}$ and hence $\left\|\nabla_{S}L_{n}(\bar{\boldsymbol{\theta}})\right\|_{\infty}\leq\lambda$ . On the other hand, by letting $\boldsymbol{\theta}_{t}=(1-t)\boldsymbol{\theta}^{*}+t\bar{\boldsymbol{\theta}}(0\leq t\leq 1)$ we have

		$\displaystyle\nabla_{S}L_{n}(\bar{\boldsymbol{\theta}})-\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)=\int_{0}^{1}\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}_{t}\right)\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)_{S}\mathrm{d}t$
		$\displaystyle=\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)_{S}+\int_{0}^{1}\left[\nabla_{SS}^{2}L_{n}\left(\bar{\boldsymbol{\theta}}_{t}\right)_{S}-\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right]\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)_{S}\mathrm{d}t.$

Hence

		$\displaystyle\left\\|\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)_{S}-\left(\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right)^{-1}\left[\nabla_{S}L_{n}(\bar{\boldsymbol{\theta}})-\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{*}\right)\right]\right\\|_{\infty}$
		$\displaystyle\leq\int_{0}^{1}\left\\|\left(\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right)^{-1}\left[\nabla_{SS}^{2}L_{n}\left(\bar{\boldsymbol{\theta}}_{t}\right)-\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right]\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right)_{S}\right\\|_{\infty}\mathrm{d}t$
		$\displaystyle\leq\left\\|\left(\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right)^{-1}\right\\|_{\infty}\sup_{t\in[0,1]}\left\\|\nabla_{SS}^{2}L_{n}\left(\bar{\boldsymbol{\theta}}_{t}\right)-\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\\|_{\infty}$

By (i) (ii) (iii) in Lemma 5, we obtain that

\left\|\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right)_{S}-\left(\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{*}\right)\right)^{-1}\left[\nabla_{S}L_{n}(\bar{\boldsymbol{\theta}})-\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{*}\right)\right]\right\|_{\infty}\leq\frac{M}{2\kappa_{\infty}}\left\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\|_{2}\left\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\|_{\infty}.

By $\bar{\boldsymbol{\theta}}\in B_{S}\left(\boldsymbol{\theta}^{*},A_{1}\right)$ we have

	$\displaystyle\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\\|_{\infty}$	$\displaystyle\leq\left\\|\left(\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right)^{-1}\right\\|_{\infty}\left\\|\nabla_{S}L_{n}(\bar{\boldsymbol{\theta}})-\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}+\frac{M}{2\kappa_{\infty}}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{2}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{\infty}$
		$\displaystyle\leq\frac{1}{2\kappa_{\infty}}\left(\lambda+\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}\right)+\frac{1}{6}\left\\|\overline{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{\infty}.$

Therefore,

\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\|_{\infty}\leq\frac{3}{5\kappa_{\infty}}\left(\left\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\|_{\infty}+\lambda\right).

(28)

Third we study the $L_{1}$ bound. The bound on $\left\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\|_{1}$ can be obtained in a similar way. Using the fact that $\|\cdot\|_{1}=\|\cdot\|_{\infty}$ for symmetric matrices,

	$\displaystyle\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\\|_{1}$	$\displaystyle\leq\left\\|\left(\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right)^{-1}\right\\|_{1}\left\\|\nabla_{S}L_{n}(\bar{\boldsymbol{\theta}})-\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{1}+\frac{M}{2\kappa_{\infty}}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{2}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{1}$
		$\displaystyle\leq\frac{1}{2\kappa_{\infty}}\left(\lambda\left\|S_{1}\right\|+\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{1}\right)+\frac{1}{6}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{1}.$

Hence $\left\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\|_{1}\leq\frac{3}{5\kappa_{\infty}}\left(\left\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\|_{1}+\lambda\left|S_{1}\right|\right)$ . Since $\operatorname{supp}(\bar{\boldsymbol{\theta}})\subseteq S$ , we also have

\left\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\|_{1}\leq\sqrt{|S|}\left\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\|_{2}\leq\frac{2\sqrt{|S|}}{\kappa_{2}}\left(\left\|\nabla_{S}L\left(\boldsymbol{\theta}^{*}\right)\right\|_{2}+\lambda\sqrt{\left|S_{1}\right|}\right).

This gives another $L_{1}$ bound. By Lemma 10, to derive $\hat{\boldsymbol{\theta}}=\bar{\boldsymbol{\theta}}$ it remains to show that $\left\|\nabla_{S_{2}}L_{n}(\bar{\boldsymbol{\theta}})\right\|_{\infty}<\lambda$ . Using the Taylor expansion we have

		$\displaystyle\nabla_{S_{2}}L_{n}(\bar{\boldsymbol{\theta}})-\nabla_{S_{2}}L_{n}\left(\boldsymbol{\theta}^{}\right)=\int_{0}^{1}\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}_{t}\right)\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)_{S}\mathrm{d}t$		(29)
		$\displaystyle=\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)+\int_{0}^{1}\left[\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}_{t}\right)-\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right]\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)_{S}\mathrm{d}t.$		(29)

On the one hand, the first term in (29) follows,

	$\displaystyle\left\\|\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)\right\\|_{\infty}$	$\displaystyle=\left\\|\left[\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\left(\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right)^{-1}\right]\left[\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)\right]\right\\|_{\infty}$
		$\displaystyle\leq(1-\tau)\left\\|\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)\right\\|_{\infty}.$

By the Taylor expansion, triangle’s inequality, (i) in Lemma 5 and the fact that $\bar{\boldsymbol{\theta}}\in B_{S}\left(\boldsymbol{\theta}^{*},A_{1}\right)$ ,

	$\displaystyle\left\\|\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)\right\\|_{\infty}$	$\displaystyle\leq\left\\|\nabla_{S}L_{n}(\bar{\boldsymbol{\theta}})-\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\\|_{\infty}$
		$\displaystyle+\int_{0}^{1}\left\\|\left[\nabla_{SS}^{2}L_{n}\left(\bar{\boldsymbol{\theta}}_{t}\right)-\nabla_{SS}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right]\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)\right\\|_{\infty}\mathrm{d}t$
		$\displaystyle\leq\left\\|\nabla_{S}L_{n}(\bar{\boldsymbol{\theta}})\right\\|_{\infty}+\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}+M\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{2}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\\|_{\infty}$
		$\displaystyle\leq\lambda+\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}+\frac{\kappa_{\infty}\tau}{3}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{\infty}.$

On the other hand, we bound the second term in (B.4). Note that $\boldsymbol{\theta}_{t}\in B_{S}\left(\boldsymbol{\theta}^{*},A_{1}\right)$ for all $t\in[0,1]$ . (i) in Lemma 5 yields

		$\displaystyle\left\\|\int_{0}^{1}\left[\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}_{t}\right)-\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right]\left(\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right)\mathrm{d}t\right\\|_{\infty}$
		$\displaystyle\leq\sup_{t\in[0,1]}\left\\|\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}_{t}\right)-\nabla_{S_{2}S}^{2}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{\infty}\leq\frac{\kappa_{\infty}\tau}{3}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\\|_{\infty}.$

As a result,

	$\displaystyle\left\\|\nabla_{S_{2}}L_{n}(\overline{\boldsymbol{\theta}})\right\\|_{\infty}$	$\displaystyle\leq\left\\|\nabla_{S_{2}}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}+(1-\tau)\left(\lambda+\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}+\frac{\kappa_{\infty}\tau}{3}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\\|_{\infty}\right)$
		$\displaystyle+\frac{\kappa_{\infty}\tau}{3}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*}\right\\|_{\infty}$
		$\displaystyle\leq\lambda-\tau\left(\lambda-\frac{2\kappa_{\infty}}{3}\left\\|\bar{\boldsymbol{\theta}}-\boldsymbol{\theta}^{}\right\\|_{\infty}-\frac{2}{\tau}\left\\|\nabla L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}\right).$

Recall that the $L_{\infty}$ bound in (28). By plugging in this estimate, and using the assumptions $0<\tau<1$ and $\lambda>\frac{20}{3\tau}\left\|\nabla L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\|_{\infty}$ , we derive that

	$\displaystyle\left\\|\nabla_{S_{2}}L_{n}(\bar{\boldsymbol{\theta}})\right\\|_{\infty}$	$\displaystyle\leq\lambda-\tau\left(\lambda-\frac{2}{5}\left(\left\\|\nabla_{S}L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}+\lambda\right)-\frac{2}{\tau}\left\\|\nabla L_{n}\left(\boldsymbol{\theta}^{}\right)\right\\|_{\infty}\right)$
		$\displaystyle\leq\lambda-\tau\left(\frac{3}{5}\lambda-\frac{4}{\tau}\left\\|\nabla L_{n}\left(\boldsymbol{\theta}^{*}\right)\right\\|_{\infty}\right)<\lambda.$

This implies $\widehat{\boldsymbol{\theta}}=\bar{\boldsymbol{\theta}}$ and translates all the bounds for $\bar{\boldsymbol{\theta}}$ to the ones for $\widehat{\boldsymbol{\theta}}$ . The proposition on sign consistency follows from elementary computation, thus we omit its proof.

Proof of Theorem 3

Proof: Write the loss function 3.1 as $L_{n}\left(\bm{y},\widehat{\bm{W}}\bm{\theta}\right)=\frac{1}{n}\sum_{i=1}^{n}L\left(y_{i},\widehat{\bm{w}}_{i}^{\top}\bm{\theta}\right)$ . We note $L\left(y_{i},\widehat{\bm{w}}_{i}^{\top}\bm{\theta}\right)|_{\bm{A}}$ to represent that $\widehat{\bm{f}}_{i}$ and $\widehat{\bm{u}}_{i}$ of the loss function are calculated by $\bm{A}$ , so as $L\left(y_{i},\widehat{\bm{w}}_{i}^{\top}\bm{\theta}\right)|_{\widehat{\bm{A}}}$ . Note $\zeta_{i}=L\left(y_{i},\widehat{\bm{w}}_{i}^{\top}\bm{\theta}\right)|_{\widehat{\bm{A}}}-L\left(y_{i},\widehat{\bm{w}}_{i}^{\top}\bm{\theta}\right)|_{{\bm{A}}}$ . For a fixed $n$ , $\zeta_{i},i=1,\ldots,n$ have the same distribution.

\displaystyle\text{Var}(\frac{1}{n}\sum_{i=1}^{n}\zeta_{i})\leq\mathbb{E}(\zeta_{i}^{2}).

With $\zeta_{i}\rightarrow 0$ in probability, we have $\mathbb{E}(\zeta_{i}^{2})\rightarrow 0$ .

	$\displaystyle\\|{\bm{f}}_{i}^{\top}(\omega)-\bm{f}_{i}^{\top}\\|_{2}^{2}$	$\displaystyle=O_{P}(\omega^{2}\\|\bm{a}_{i}\\|_{2}^{2}+\\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\\|_{2}^{2}),~{}\quad{\mbox{and}}$
	$\displaystyle\\|{\bm{u}}_{i}^{\top}(\omega)-\bm{u}_{i}^{\top}\\|_{2}^{2}$	$\displaystyle=O_{P}(\omega^{2}\\|\bm{a}_{i}\\|_{2}^{2}+\\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\bm{B}^{\top}\\|_{2}^{2})$

$\displaystyle\bm{y}$	$\displaystyle=\left(\widehat{\bm{A}}-\mathbb{E}(\widehat{\bm{A}})\right)\bm{H}^{}+\left(\bm{A}-\widehat{\bm{A}}+\mathbb{E}(\widehat{\bm{A}})\right)\bm{H}^{}{}+\tilde{\bm{\epsilon}}$
	$\displaystyle=\widehat{\bm{F}}\widehat{\bm{B}}^{\top}\bm{H}^{}{}+\widehat{\bm{U}}\bm{H}^{}{}+\left(\bm{A}-\widehat{\bm{A}}+\mathbb{E}(\widehat{\bm{A}})\right)\bm{H}^{*}{}+\tilde{\bm{\epsilon}}$
	$\displaystyle=\widehat{\bm{F}}\widehat{\bm{B}}^{\top}\bm{H}^{}{}+\widehat{\bm{U}}\bm{H}^{}{}+\tilde{\bm{e}}+\tilde{\bm{\epsilon}},$	(18)

	$\displaystyle\\|\widehat{\bm{f}}_{i}^{\top}-\bm{f}_{i}^{\top}\\|_{2}^{2}$	$\displaystyle=O(\omega^{2}+\\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\\|_{2}^{2}+1/n),$
	$\displaystyle\\|\widehat{\bm{u}}_{i}^{\top}-\bm{u}_{i}^{\top}\\|_{2}^{2}$	$\displaystyle=O(\omega^{2}+\\|\bm{u}_{i}^{\top}\bm{B}\Lambda_{B}^{-1}\bm{B}^{\top}\\|_{2}^{2}+1/n),$

	$\displaystyle\|\bm{a}_{i}^{\top}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{k}^{(B)}I)^{+}\Lambda_{\bm{u}}\bm{\xi}_{k}^{(B)}\omega\|^{2}$	$\displaystyle\leq(\frac{\lambda_{\mathop{\rm max}}^{\bm{u}}}{\min_{i<j\leq K}(\lambda_{i}^{(B)}-\lambda_{j}^{(B)})})^{2}\\|\bm{a}_{i}\\|_{2}^{2}\omega^{2}$
	$\displaystyle\\|\bm{a}_{i}^{\top}((\mathbf{B}\mathbf{B}^{\top}-\lambda_{k}^{(B)}I)^{+}\Lambda_{\bm{u}}\bm{\xi}_{k}^{(B)}\omega)_{k=1,\ldots,K}{\Lambda^{C}}^{-1/2}\\|_{2}^{2}$	$\displaystyle\leq\frac{1}{\lambda^{(C)}_{\min}(\omega)}(\frac{\lambda_{\mathop{\rm max}}^{\bm{u}}}{\min_{i<j\leq K}(\lambda_{i}^{(B)}-\lambda_{j}^{(B)})})^{2}\\|\bm{a}_{i}\\|_{2}^{2}\omega^{2}.$

	$\displaystyle\\|\bm{a}_{i}^{\top}\sum_{j=1}^{K}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\Lambda_{u}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)\top}\Lambda_{u}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\\|_{2}^{2}$
	$\displaystyle\leq\\|\sum_{j=1}^{K}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\Lambda_{u}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)\top}\Lambda_{u}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\\|_{2}^{2}\\|\bm{a}_{i}\\|_{2}^{2}$
	$\displaystyle\leq(\lambda^{u}_{max})^{4}\\|\sum_{j=1}^{K}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\bm{\xi}_{j}^{(B)}\bm{\xi}_{j}^{(B)\top}(\mathbf{B}\mathbf{B}^{\top}-\lambda_{j}^{(B)}I)^{+}\\|_{2}^{2}\\|\bm{a}_{i}\\|_{2}^{2}$
	$\displaystyle\leq(p-K)(\frac{\lambda_{\mathop{\rm max}}^{\bm{u}}}{\min_{i<j\leq K}(\lambda_{i}^{(B)}-\lambda_{j}^{(B)})})^{4}\\|\bm{a}_{i}\\|_{2}^{2}.$

Sparsity learning via structured functional factor augmentation

Abstract

1 Introduction

2 A novel functional factor augmentation structure

2.1 Functional observations and its expansion

Lemma 1

Remark 1

2.2 A functional factor augmentation structure

Example 1

Lemma 2

Remark 2

2.3 Properties of functional factor augmentation structure

Lemma 3

Remark 3

Assumption 1

Assumption 2

Lemma 4

Theorem 1

2.4 Truncation analysis and the relationship between KK and mgm_{g}

3 The functional factor augmentation selection model

3.1 The functional linear regression model

3.2 Theoretical justifications for functional variable selection

Assumption 3

Assumption 4

Assumption 5

Theorem 2

Theorem 3

Corollary 1

4 Simulation Studies

5 Real Data Application

5.1 Effects of macroeconomic covariates on lifetime expectancy

5.2 Prediction of sales areas of commercial houses

6 Conclusion

References

7 Regularity Condition

8 Proof of Theorem

Lemma 5

Lemma 6

Lemma 7

Lemma 8

Lemma 9

Lemma 10

2.4 Truncation analysis and the relationship between $K$ and $m_{g}$