Kernel-based estimation for partially functional linear model: Minimax rates and randomized sketches

This paper makes three main contributions to this functional modeling literature.

Our first contribution is to establish Theorem 1 stating that with high probability, under mild regularity conditions, the prediction error of our procedure under the squared $L_{2}$-norm is bounded by $\big{(}\frac{p_{0}\log p}{n}+n^{-\frac{2r}{2r+1}}\big{)}$, where the quantity $r>1/2$ corresponds to the kernel complexity of one composition kernel $K^{1/2}CK^{1/2}$. The proof of this upper bound involves two different penalties for analyzing the obtained estimator in high dimensions, and we want to emphasize that it is very hard to prove constraint cone set that has often been used to define a critical condition (constraint eigenvalues constant) for high-dimensional problems (Bickel, Ritov, and Tsybakov, 2009; Verzelen, 2012). To handle this technical difficulty, we combine the methods used in Müller and Van de Geer (2015) for high dimensional partial linear models with various techniques in empirical process theory for analyzing kernel classes (Aronszajn, 1950; Cai and Yuan, 2012; Yuan and Cai, 2010; Zhu et al., 2014).

Our second contribution is to establish algorithm-independent minimax lower bounds under the squared $L_{2}$ norm. These minimax lower bounds, stated in Theorem 2, are determined in terms of the metric entropy of the composition kernel $K^{1/2}CK^{1/2}$ and the sparsity structure of high dimensional scalar coefficients. For the commonly-used kernels, including the Sobolev classes, these lower bounds match our achievable results, showing optimality of our estimator for PFLM. It is worthy noting that, the lower bound of parametric part does not depend on nonparametric smoothness indices, coinciding with the classical sparse estimation rate in the high dimensional linear models (Verzelen, 2012). By contrast, the lower bound for estimating $f^{*}$ turns out to be affected by the regression coefficient $\boldsymbol{\gamma}^{*}$. The proof of Theorem 2 is based on characterizing the packing entropies of the class of nonparametric kernel models, interaction between the composition kernel and high dimensional scalar vector, combined with classical information theoretic techniques involving Fano’s inequality and variants (Yang and Barron, 1999; Van. de. Geer, 2000; Tsybakov, 2009).

Our third contribution is to consider randomized sketches for our original estimation with statistical dimension. Despite these attractive statistical properties stated as above, the computational complexity of computing our original estimate prevents it from being routinely used in large-scale problems. In fact, a standard implementation for any kernel estimator leads to the time complexity $O(n^{3})$ and space complexity $O(n^{2})$ respectively. To this end, we employ the random projection and sketching techniques developed in Yang et al. (2017); Mahoney (2011), where it is proposed to approximate $n$-dimensional kernel matrix by projecting its row and column subspaces to a randomly chosen m-dimensional subspace with $m\ll n$. We give the sketch dimension $m$ proportional to the statistical dimension, under which the resulting estimator has a comparable numerical performance.

A class of conventional estimation procedures for functional linear regressions in the statistical literature are based on functional principal components regression (FPCA) or spline functions; see (Ramsay and Silverman, 2005; Ferraty and Vieu, 2006; Kong, Xue, Yao, and Zhang, 2016) and (Cardot, Ferraty, and Sarda, 2003) for details. These truncation approaches to handle an infinity-dimensional function only depend on the information of the feature $X$. In particular, commonly-used FPCA methods that form an available basis for the slope function $f^{*}$ are determined solely by empirical covariance of the observed feature $X$, and these basis may not act as an efficient representation to approximate $f^{*}$, since the slope function $f^{*}$ and the leading functional components are essentially unrelated. Similar problems also arise when spline-based finite representation are used.

To avoid inappropriate representation for the slope function, reproducing kernel methods have been known to be a family of powerful tools for directly estimating infinity-dimensional functions. When the slope function is assumed to reside in a reproducing kernel Hilbert Space (RKHS), denoted by $\mathcal{H}_{K}$, several existing work (Yuan and Cai, 2010; Cai and Yuan, 2012; Zhu, Yao, and Zhang, 2014) for functional linear or additive regression have proved that the minimax rate of convergence depends on both the kernel $K$ and the covariance function $C$ of the functional feature $X$. In particular, the alignment of $K$ and $C$ can significantly affect the optimal rate of convergence. However, it is well known that kernel-based methods suffer a lot from storage cost and computational burden. Specially, kernel-based methods need to store a $n\times n$ matrix before running algorithms and are limited to small-scale problems.

The rest of this paper is organized as follows. Section 2 introduces some notations and the basic knowledge on kernel method, and formulates the proposed kernel-based regularized estimation method. Section 3 is devoted to establish the minimax rate of the prediction problem for PFLM and provide detailed discussion on the obtained results, including the desired convergence rate of the upper bounds and a matching set of minimax lower bounds. In Section 4, a general sketching-based strategy is provided, and an approximate algorithm for solving (2.2) is employed. Several numerical experiments are implemented in Section 5 to support the proposed approach and the employed optimization strategy. A brief summary of this paper is provided in Section 6. Appendix A contains several core proof procedures of the main results, including the technical proofs of Theorems 1–3. Some useful lemmas and more technical details are provided in Appendix B.

Let $u,v$ be two general random variables, and denote the joint distribution of $(u,v)$ by $Q$ and the marginal distribution of $u(z)$ by $Q_{u}(Q_{v})$. For a measurable function $f:\,u\times v\rightarrow\mathbb{R}$, we define the squared $L_{2}$-norm by $\|f\|^{2}:=\mathbb{E}_{Q}f^{2}(u,v)$, and the squared empirical norm is given by $\|f\|_{n}^{2}:=\frac{1}{n}\sum_{i=1}^{n}f^{2}(u_{i},v_{i})$, where $\{(u_{i},v_{i})\}_{i=1}^{n}$ are i.i.d. copies of $(u,v)$. Note that $Q$ may differ from line to line. For a vector $\boldsymbol{\gamma}\in\mathbb{R}^{p}$, the $\ell_{1}$-norm and $\ell_{2}$-norm are given by $\|\boldsymbol{\gamma}\|_{1}:=\sum_{j=1}^{p}|\gamma_{j}|$ and $\|\boldsymbol{\gamma}\|_{2}:=\big{(}\sum_{j=1}^{p}\gamma_{j}^{2}\big{)}^{1/2}$, respectively. With a slight abuse of notation, we write $\|f\|_{\mathcal{L}_{2}}:=\langle f,f\rangle_{\mathcal{L}_{2}}$ with $\langle f,g\rangle_{\mathcal{L}_{2}}=\int_{\mathcal{T}}f(t)g(t)dt$. For two sequences $\{a_{k}:k\geq 1\}$ and $\{b_{k}:k\geq 1\}$, $a_{k}\lesssim b_{k}$ (or $a_{k}=O(b_{k})$) means that there exists some constant $c$ such that $a_{k}\leq cb_{k}$ for all $k\geq 1$. Also, we write $a_{k}\gtrsim b_{k}$ if there is some positive constant $c$ such that $a_{k}\geq cb_{k}$ for all $k\geq 1$. Accordingly, we write $a_{k}\asymp b_{k}$ if both $a_{k}\lesssim b_{k}$ and $a_{k}\gtrsim b_{k}$ are satisfied.

Kernel methods are one of the most powerful learning schemes in machine learning, which often take the form of regularization schemes in a reproducing kernel Hilbert space (RKHS) associated with a Mercer kernel (Aronszajn, 1950). A major advantage of employing the kernel methods is that the corresponding optimization task over an infinite dimensional RKHS are equivalent to a $n$-dimensional optimization problems, benefiting from the so-called reproducing property.

Recall that a kernel $K(\cdot,\cdot):\mathcal{T}\times\mathcal{T}\rightarrow{\cal R}$ is a continuous, symmetric, and positive semi-definite function. Let $\mathcal{H}_{{K}}$ be the closure of the linear span of functions $\{K_{t}(\cdot):={K}(t,\cdot),t\in\mathcal{T}\}$ endowed with the inner product $\langle\sum_{i=1}^{n}\alpha_{i}{K}_{t_{i}},\,\sum_{j=1}^{n}\beta_{j}{K}_{t_{j}}\rangle_{{K}}:=\sum_{i,j=1}^{n}\alpha_{i}\beta_{j}{K}(t_{i},t_{j})$, for any $\{t_{i}\}_{i=1}^{n},\{t_{i}\}_{i=1}^{n}\in\mathcal{T}^{n}$ and $n\in\mathcal{N^{+}}$. An important property on $\mathcal{H}_{{K}}$ is the reproducing property stating that

$$f(t)=\langle f,{K}_{t}\rangle_{{K}},\,\,\,\hbox{for any}\,f\in\mathcal{H}_{{K}}.$$

This property ensures that an RKHS inherits many nice properties from the standard finite dimensional Euclidean spaces. Throughout this paper, we assume that the slope function $f^{*}$ resides in a specified RKHS, still denoted by $\mathcal{H}_{K}$. In addition, another RKHS can be naturally induced by the stochastic process of $X(\cdot)$. Without loss of generality, we assume that $X(\cdot)$ is square integrable over $\mathcal{T}$ with zero-mean, ant thus the covariance function of $X$, defining as

$$C(s,t)=\mathbb{E}[X(s)X(t)],\quad\forall\,t,\,s\in\mathcal{T},$$

is also a real, semi-definite kernel.

Note that the kernel complexity is characterized explicitly by a kernel-induced integral operator. Precisely, for any kernel ${K(\cdot,\cdot)}:\mathcal{T}\times\mathcal{T}\rightarrow{\cal R}$, we define the integral operator $L_{{K}}:\mathcal{L}_{2}(\mathcal{T})\rightarrow\mathcal{L}_{2}(\mathcal{T})$ by

$$L_{{K}}(f)(\cdot)=\int_{\mathcal{T}}{K}(s,\cdot)f(s)ds.$$

By the reproducing property, $L_{{K}}$ can be equivalently defined as

$$\langle f,L_{{K}}(g)\rangle_{K}=\langle f,g\rangle_{\mathcal{L}_{2}},\quad\forall\,f\in\mathcal{H}_{{K}},\,g\in\mathcal{L}_{2}(\mathcal{T}).$$

Since the operator $L_{{K}}$ is linear, bounded and self-adjoint in $\mathcal{L}_{2}(\mathcal{T})$, the spectral theorem implies that there exists a family of orthonormalized eigenfunctions $\{\phi^{{K}}_{\ell}:\,\ell\geq 1\}$ and a sequence of eigenvalues $\theta_{1}^{{K}}\geq\theta_{2}^{{K}}\geq...>0$ such that

$${K}(s,t)=\sum_{\ell\geq 1}\theta_{\ell}^{{K}}\phi^{{K}}_{\ell}(s)\phi^{{K}}_{\ell}(t),\quad s,\,t\in\mathcal{T},$$

and thus by definition, it holds

$$L_{{K}}(\phi^{{K}}_{\ell})=\theta_{\ell}^{{K}}\phi^{{K}}_{\ell},\quad\ell=1,2,...$$

Based on the semi-definiteness of $L_{{K}}$, we can always decompose it into the following form

$$L_{{K}}=L_{{K}}^{1/2}\circ L_{{K}}^{1/2},$$

where $L_{{K}^{1/2}}$ is also a kernel-induced integral operator associated with a fractional kernel ${K}^{1/2}$ that

$${K}^{1/2}(s,t):=\sum_{\ell\geq 1}\sqrt{\theta_{\ell}^{{K}}}\phi^{{K}}_{\ell}(s)\phi^{{K}}_{\ell}(t),\quad s,\,t\in\mathcal{T}.$$

Also, it holds

$$L_{{K}^{1/2}}(\phi^{K}_{\ell}):=\sqrt{\theta_{\ell}^{{K}}}\phi^{{K}}_{\ell}.$$

Given two kernels $K_{1},K_{2}$, we define

$$(K_{1}K_{2})(s,t):=\int_{\mathcal{T}}K_{1}(s,u)K_{2}(t,u)du,$$

and then it holds $L_{K_{1}K_{2}}=L_{K_{1}}\circ L_{K_{2}}$. Note that $K_{1}K_{2}$ is not necessarily a symmetric kernel.

In the rest of this paper, we focus on the RKHS $\mathcal{H}_{K}$ in which the slope function $f^{*}$ in (1.2) resides. Given the kernel $K$, the covariance function $C$ and by using the above notation, we define the linear operator $L_{K^{1/2}C_{k}K^{1/2}}$ by

$$L_{K^{1/2}CK^{1/2}}:=L_{K^{1/2}}\circ L_{C}\circ L_{K^{1/2}}.$$

If the both operators $L_{K^{1/2}}$ and $L_{C}$ are linear, bounded and self-adjoint, so is $L_{K^{1/2}CK^{1/2}}$. By the spectral theorem, there exist a sequence of positive eigenvalues $s_{1}\geq s_{2}\geq...>0$ and a set of orthonormalied eigenfunctions $\{\varphi_{\ell}:\ell\geq 1\}$ such that

$$K^{1/2}CK^{1/2}(s,t)=\sum_{\ell\geq 1}s_{\ell}\varphi_{\ell}(s)\varphi_{\ell}(t),\quad\forall\,s,t\in\mathcal{T},$$

and particularly

$$L_{K^{1/2}CK^{1/2}}(\varphi_{\ell})=s_{\ell}\varphi_{\ell},\quad\ell=1,2,...$$

It is worthwhile to note that the eigenvalues $\{s_{\ell}:\ell\geq 1\}$ of the linear operator $L_{K^{1/2}CK^{1/2}}$ depend on the eigenvalues of both the reproducing kernel $K$ and the covariance function $C$. We shall show in Section 3 that the minimax rate of convergence of the excess prediction risk is determined by the decay rate of the eigenvalues $\{s_{\ell}:\ell\geq 1\}$.

Given the sample $\{Y_{i},(X_{i},{\bf Z}_{i})\}_{i=1}^{n}$ which are independently drawn from (1.2), the proposed estimation procedure for PFLM (1.2) is formulated in a least square regularization scheme by solving

$$(\widehat{f},\widehat{\boldsymbol{\gamma}})=\mathop{\rm argmin}_{f\in\mathcal{H}_{K},\boldsymbol{\gamma}\in{\cal R}^{p}}\Big{\{}\frac{1}{n}\sum_{i=1}^{n}\big{(}Y_{i}-\langle X_{i},f\rangle_{\mathcal{L}_{2}}-{\bf Z}^{T}_{i}\boldsymbol{\gamma}\big{)}^{2}+\mu^{2}\|f\|^{2}_{K}+\lambda\|\boldsymbol{\gamma}\|_{1}\Big{\}},$$

(2.1)

where the parameter $\mu^{2}>0$ is used to control the smoothness of the nonparametric component and $\lambda>0$ associated with the $\ell_{1}$-type penalty is used to generate sparsity with respect to the scalar covariates.

Note that although the proposed estimation procedure (2.1) is formulated within an infinite-dimensional Hilbert space, the following lemma shows that this optimization task is equivalent to a finite-dimensional minimization problem.

To rewrite the minimization problem (2.1) into a matrix form, we define a $n\times n$ semi-definite matrix $\mathbb{K}^{c}=(K^{c}_{ik})_{i,k=1}^{n}$ with $K^{c}_{ik}:=\langle X_{i},B_{k}\rangle_{\mathcal{L}_{2}}=\iint X_{k}(u)X_{i}(t)K(t,u)dudt$, and by the reproducing property on $K$, we also get $\langle B_{i},B_{k}\rangle_{K}=K^{c}_{ik}$, $i,k=1,...,n$. Thus, by Lemma 1, the matrix form of (2.1) is given as

$\displaystyle\min_{\boldsymbol{\alpha}\in{\cal R}^{n},\boldsymbol{\gamma}\in{\cal R}^{p}}\frac{1}{n}\big{\|}\mathbf{y}-\mathbb{K}^{c}\boldsymbol{\alpha}-\mathbb{Z}\boldsymbol{\gamma}\big{\|}_{2}^{2}+\mu^{2}\boldsymbol{\alpha}^{T}\mathbb{K}^{c}\boldsymbol{\alpha}+\lambda\|\boldsymbol{\gamma}\|_{1},$

(2.2)

where $\mathbb{Z}\in{\cal R}^{n\times p}$ denotes the design matrix of $\bf Z$. Since the unconstrained problem (2.2) is convex for both $\boldsymbol{\alpha}$ and $\boldsymbol{\gamma}$, the standard alternative optimization (Boyd and Vandenberghe, 2004) can be applied directly to approximate a global minimizer of (2.2). Yet, due to the fact that $\mathbb{K}^{c}$ is a $n\times n$ matrix, both computation cost and storage burden are very heavy in standard implementation, with the orders $O(n^{3})$ and $O(n^{2})$ respectively. To alleviate the computational issue, we propose an approximate numerical optimization instead of (2.2) in Section 4. Precisely, a class of general random projections are adopted to compress the original kernel matrix $\mathbb{K}^{c}$ and improve the computational efficiency.

We denote the short-hand notation

$$\mathcal{G}:=\big{\{}g=\langle X,f\rangle_{\mathcal{L}_{2}}+{\bf Z}^{T}\boldsymbol{\gamma},\,\,f\in\mathcal{H}_{K},\,\boldsymbol{\gamma}\in{\cal R}^{p}\big{\}},$$

and the functional $g^{*}({\bf U}):=\langle X,f^{*}\rangle_{\mathcal{L}_{2}}+{\bf Z}^{T}\boldsymbol{\gamma}^{*}$ for ${\bf U}=(X,{\bf Z})$. With a slight confusion of notation, we sometimes also write $\mathcal{G}:=\big{\{}g=(f,\boldsymbol{\gamma}),\,\,f\in\mathcal{H}_{K},\,\boldsymbol{\gamma}\in{\cal R}^{p}\big{\}}$. To split the scalar components and the functional component involved in our analysis, we define the projection of $Z_{j}$ concerning $\mathcal{H}_{K}$ as $\Pi(Z_{j}|\mathcal{H}_{K})=\mathop{\rm argmin}_{f\in\mathcal{H}_{K}}\|Z_{j}-\langle X,f\rangle_{\mathcal{L}_{2}}\|^{2}$. Let $\Pi(Z_{j}|X)=\langle X,\Pi(Z_{j}|\mathcal{H}_{K})\rangle_{\mathcal{L}_{2}}$ and $\Pi_{{\bf Z}|X}=(\Pi(Z_{1}|X),...,\Pi(Z_{p}|X))^{T}$, and then we denote $\widetilde{\bf Z}:={\bf Z}-\Pi_{{\bf Z}|X}$ as a random vector of ${\cal R}^{p}$. For any $g_{1}({\bf U}):=\langle X,f_{1}\rangle_{\mathcal{L}_{2}}+{\bf Z}^{T}\boldsymbol{\gamma}_{1}\in\mathcal{G}$ and $g_{2}({\bf U}):=\langle X,f_{2}\rangle_{\mathcal{L}_{2}}+{\bf Z}^{T}\boldsymbol{\gamma}_{2}\in\mathcal{G}$, we have the following orthogonal decomposition that

	$\displaystyle g_{1}({\bf U})-g_{2}({\bf U})$	$\displaystyle={\bf Z}^{T}(\boldsymbol{\gamma}_{1}-\boldsymbol{\gamma}_{2})+\langle X,f_{2}-f_{1}\rangle_{\mathcal{L}_{2}}$
		$\displaystyle=\widetilde{\bf Z}^{T}(\boldsymbol{\gamma}_{1}-\boldsymbol{\gamma}_{2})+\Pi_{{\bf Z}^{T}|X}^{T}(\boldsymbol{\gamma}_{1}-\boldsymbol{\gamma}_{2})+\langle X,f_{2}-f_{1}\rangle_{\mathcal{L}_{2}},$

and by the definition of projection, it holds

$\displaystyle\|g_{1}-g_{2}\|^{2}=\|\widetilde{\bf Z}^{T}(\boldsymbol{\gamma}_{1}-\boldsymbol{\gamma}_{2})\|^{2}+\|\Pi_{{\bf Z}|X}^{T}(\boldsymbol{\gamma}_{1}-\boldsymbol{\gamma}_{2})+\langle X,f_{2}-f_{1}\rangle_{\mathcal{L}_{2}}\|^{2}.$

(3.1)

To establish the refined upper bounds of the prediction and estimation errors, we summarize and discuss the main conditions needed in the theoretical analysis below.

Condition A is commonly used in literature of semiparametric modelling ; see (Müller and Van de Geer, 2015) for reference. This condition ensures that there is enough information in the data to identify the parameters in the scalar part. Condition B imposes some boundedness assumptions, which are not essential and are used only for simplifying the technical proofs. Condition C implies that the random process $L_{K^{1/2}}X$ has an exponential decay rate and the same condition is also considered in Cai and Yuan (2012). Particularly, it is naturally satisfied if $X$ is a Gaussian process. In Condition D, the parameters $s_{\ell}$ are related to the alignment between $K$ and $C$, which plays an important role in determining the minimax optimal rates. Moreover, the decay of $s_{\ell}$ characterizes the kernel complexity and has close relation with various covering numbers and Radmeacher complexity. Specially, the polynomial decay assumed in Condition D is satisfied for the classical Sobolev class and Besov class.

The following theorem states that with an appropriately chosen $(\mu,\lambda)$, the predictor $\widehat{g}:=\langle X,\widehat{f}\rangle_{\mathcal{L}_{2}}+{\bf Z}^{T}\widehat{\boldsymbol{\gamma}}$ attains a sharp convergence rate under $L_{2}$-norm.

Theorem 1 shows that the proposed estimation (2.1) achieves a fast convergence rate in the term of prediction error. Note that the derived rate depends on the kernel complexity of $K^{1/2}CK^{1/2}$ and the sparsity of scalar components. It is interesting to note that even there exists some underlying correlation structure between the functional feature and the scalar covariates, the choice of hyper-parameter $\mu$ depends on the structural information of all the features, while the sparsity hyper-parameter $\lambda$ only depends on the scalar component.

It is worthy pointing out that the estimation error of the parametric estimator $\widehat{\boldsymbol{\gamma}}$ can achieve the optimal convergence rate in the high dimensional linear models (Verzelen, 2012), even in the presence of nonparametric components. This result in the functional literature is similar in spirit to the classical high dimensional partial linear models (Müller and Van de Geer, 2015; Yu, Levine, and Cheng, 2019).

In this part, we establish the lower bounds on the minimax risk of estimating $\boldsymbol{\gamma}^{*}$ and $\langle X,f^{*}\rangle_{\mathcal{L}_{2}}$ separately. Let $B[p_{0},p]$ be a set of $p$-dimensional vectors with at most $p_{0}$ non-zero coordinates, and ${\cal B}_{K}$ be the unit ball of $\mathcal{H}_{K}$. Moreover, we define the risk of estimating $\boldsymbol{\gamma}^{*}$ as

$$R_{\boldsymbol{\gamma}^{*}}(p_{0},p,{\cal B}_{K}):=\inf_{\widehat{\boldsymbol{\gamma}}}\sup_{\boldsymbol{\gamma}^{*}\in B[p_{0},p],f^{*}\in{\cal B}_{K}}\mathbb{E}[\|\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}^{*}\|_{2}^{2}],$$

where ${\inf}$ is taken over all possible estimators for $\boldsymbol{\gamma}^{*}$ in model (1.2). Similarly, we define the risk of estimating $\langle X,f^{*}\rangle_{\mathcal{L}_{2}}$ as

$$R_{f^{*}}(s_{0},p,{\cal B}_{K}):=\inf_{\widehat{f}}\sup_{\boldsymbol{\gamma}^{*}\in B[p_{0},p],f^{*}\in{\cal B}_{K}}\mathbb{E}[\langle X,\widehat{f}-f^{*}\rangle_{\mathcal{L}_{2}}^{2}]=\inf_{\widehat{f}}\sup_{\boldsymbol{\gamma}^{*}\in B[p_{0},p],f^{*}\in{\cal B}_{K}}\|L_{C^{1/2}}(\widehat{f}-f^{*})\|_{\mathcal{L}_{2}}^{2}.$$

The following theorem provides the lower bounds of the minimax optimal estimation error for $\boldsymbol{\gamma}^{*}$ and the predictor error for $f^{*}$, respectively.

The proof of Theorem 3 is provided in Appendix A. As mentioned previously, these results indicate that the best possible estimation of $\boldsymbol{\gamma}^{*}$ is not affected by the existence of nonparametric components, while the minimax risk for estimating the (nonparametric) slope function not only depends on the smoothness itself, but also on the dimensionality and sparsity of the scalar covariates. From the lower bound of $R_{f^{*}}(p_{0},p,{\cal B}_{K})$, we observe that a rate-switching phenomenon occurring between a sparse regime and a smooth regime. Particularly when $\frac{p_{0}\log(p/p_{0})}{n}$ dominates $n^{-\frac{2r}{2r+1}}$ corresponding to the sparse regime, the lower bound becomes the classical high dimensional parametric rate $\frac{p_{0}\log(p/p_{0})}{n}$. Otherwise, this corresponds to the smooth regime and thus has similar behaviors as classical nonparametric models. We also notice that the minimax lower bound obtained for the predictor error generalizes the previous results for the pure functional linar model (Cai and Yuan, 2012).

This section provides the detailed computational issues of the proposed approach. Precisely, we aim to solve the following optimization task that

		$\displaystyle(\widehat{\boldsymbol{\alpha}}_{s},\widehat{\boldsymbol{\gamma}}_{s}):=\mathop{\rm argmin}_{\boldsymbol{\alpha}\in{\cal R}^{m},\boldsymbol{\gamma}\in{\cal R}^{p}}\frac{1}{n}\boldsymbol{\alpha}^{T}(\mathbb{S}\mathbb{K}^{c})(\mathbb{S}\mathbb{K}^{c})^{T}\boldsymbol{\alpha}{-\frac{2}{n}\boldsymbol{\alpha}^{T}\mathbb{S}\mathbb{K}^{c}(\mathbf{y}-\mathbb{Z}\boldsymbol{\gamma})}+$
		$\displaystyle\hskip 113.81102pt{\frac{1}{n}(\mathbf{y}-\mathbb{Z}\boldsymbol{\gamma})^{T}(\mathbf{y}-\mathbb{Z}\boldsymbol{\gamma})}+\mu^{2}\boldsymbol{\alpha}^{T}\mathbb{S}\mathbb{K}^{c}\mathbb{S}^{T}\boldsymbol{\alpha}+\lambda\|\boldsymbol{\gamma}\|_{1}.$		(4.2)

To solve (4.1), a splitting algorithm with proximal operator is applied, which updates the representer coefficients ${\boldsymbol{\alpha}}$ and the linear coefficients ${\boldsymbol{\gamma}}$ sequentially. Specifically, at the $t$-th iteration with current solution $(\boldsymbol{\alpha}^{t},\mbox{\boldmath$\gamma$}^{t})$, the following two optimization tasks are solved sequentially to obtain the solution of the $(t+1)$-th iteration

	$\displaystyle\boldsymbol{\alpha}^{t+1}=\mathop{\rm argmin}_{\boldsymbol{\alpha}\in{\cal R}^{m}}\Big{\{}\frac{1}{n}\boldsymbol{\alpha}^{T}(\mathbb{S}\mathbb{K}^{c})(\mathbb{S}\mathbb{K}^{c})^{T}\boldsymbol{\alpha}-\frac{2}{n}\boldsymbol{\alpha}^{T}\mathbb{S}\mathbb{K}^{c}(\mathbf{y}-\mathbb{Z}\boldsymbol{\gamma}^{t})+\mu^{2}\boldsymbol{\alpha}^{T}\mathbb{S}\mathbb{K}^{c}\mathbb{S}^{T}\boldsymbol{\alpha}\Big{\}},$		(4.3)
	$\displaystyle\mbox{\boldmath$\gamma$}^{t+1}=\mathop{\rm argmin}_{\mbox{\boldmath$\gamma$}\in{\cal R}^{p}}\Big{\{}R_{n}(\boldsymbol{\alpha}^{t+1},\mbox{\boldmath$\gamma$})+\lambda\|\boldsymbol{\gamma}\|_{1}\Big{\}},$		(4.4)

where $R_{n}(\boldsymbol{\alpha}^{t+1},\mbox{\boldmath$\gamma$}):=\frac{2}{n}({\boldsymbol{\alpha}}^{t+1})^{T}\mathbb{S}\mathbb{K}^{c}\mathbb{Z}\boldsymbol{\gamma}+{\frac{1}{n}(\mathbf{y}-\mathbb{Z}\boldsymbol{\gamma})^{T}(\mathbf{y}-\mathbb{Z}\boldsymbol{\gamma})}$.

To update $\boldsymbol{\alpha}$, it is clear that the optimization task (4.3) has an analytic solution that

$$\boldsymbol{\alpha}^{t+1}=\big{(}(\mathbb{S}\mathbb{K}^{c})(\mathbb{S}\mathbb{K}^{c})^{T}+n\mu^{2}\mathbb{S}\mathbb{K}^{c}\mathbb{S}^{\prime}\big{)}^{-1}\mathbb{S}\mathbb{K}^{c}(\mathbf{y}-\mathbb{Z}\boldsymbol{\gamma}^{t}).$$

To update $\boldsymbol{\gamma}$, we first introduce the proximal operator (Moreau, 1962), which is defined as

$\displaystyle\mbox{Prox}_{{\lambda}\|\cdot\|_{1}}({\bf u}):=\mathop{\rm argmin}_{{\bf u}}\Big{\{}\frac{1}{2}\|{\bf u}-{\bf v}\|_{2}^{2}+\lambda\|{\bf u}\|_{1}\Big{\}}.$

(4.5)

Note that the solution of optimization task (4.5) is the well-known soft-thresholding operator with solution that

$${\big{(}\mbox{Prox}_{{\lambda}\|\cdot\|_{1}}(\mathop{\bf u})\big{)}_{i}=\mathop{\rm sign}(u_{i})(|u_{i}|-{\lambda})_{+}}.$$

Then, for the optimization task (4.4), we have

$$\mbox{\boldmath$\gamma$}^{t+1}=\mbox{Prox}_{\frac{\lambda}{D}\|\cdot\|_{1}}\Big{(}\mbox{\boldmath$\gamma$}^{t}-\frac{1}{D}\nabla_{\mbox{\boldmath$\gamma$}}R_{n}(\boldsymbol{\alpha}^{t+1},\mbox{\boldmath$\gamma$}^{t})\Big{)},$$

where $D$ denotes an upper bound of the Lipschitz constant of $R_{n}(\boldsymbol{\alpha}^{t+1},\mbox{\boldmath$\gamma$}^{t})$, and compute $\nabla_{\boldsymbol{\gamma}}R_{n}(\boldsymbol{\alpha}^{t+1},\mbox{\boldmath$\gamma$}^{t})=\frac{2}{n}\mathbb{Z}^{T}(\mathbb{S}\mathbb{K}^{c})^{T}{\boldsymbol{\alpha}}^{t+1}+\frac{2}{n}\mathbb{Z}^{T}\mathbb{Z}\boldsymbol{\gamma}^{t}-\frac{2}{n}\mathbb{Z}^{T}\mathbf{y}$. We repeat the above iteration steps until $(\boldsymbol{\alpha}^{t+1},\boldsymbol{\gamma}^{t+1})$ converges.

It should be pointed out that the exact value of $D$ is often difficult to determine in large-scale problems. A common way to handle this problem is to use a backtracking scheme (Boyd and Vandenberghe, 2004) as a more efficient alternative to approximately compute an upper bound of it.

In this paper, we consider three random sketch methods, including the sub-Gaussian random sketch (GRS), randomized orthogonal system sketch (ROS) and sub-sampling random sketch (SUB). Precisely, we denote the $i$-th row of the random matrix $\mathbb{S}$ as ${\mathop{\bf s}}_{i}$ and consider three different types of ${\mathop{\bf s}}_{i}$ as follows.

In practice, we are interested in the $m\times n$ sketch matrices with $m\ll n$ to enhance computational efficiency. Note that the existence of a $n\times n$ kernel matrix in Lemma 1 is only a sufficient condition for equivalent optimization. It has been shown theoretically in the kernel regression (Yang et al., 2017) that the kernel matrix can be compressed to be the one with small size, based on some intrinsic low-dimensional notations. Despite the model difference from Yang et al. (2017), our kernel matrix $\mathbb{K}^{c}$ does not depend on the scalar covariates $\bf Z$, and thus those derived results for the kernel regression are still applicable to our case.

Consider the eigen-decomposition $\mathbb{K}^{c}=\mathbb{U}\mathbb{D}\mathbb{U}^{T}$ of the kernel matrix, where $\mathbb{U}\in{\cal R}^{n\times n}$ is an orthonormal matrix of eigenvectors, and $\mathbb{D}=\hbox{diag}\{\hat{\mu}_{1},...,\hat{\mu}_{n}\}$ is a diagonal matrix of eigenvalues, where $\hat{\mu}_{1}\geq\hat{\mu}_{2}\geq...\geq\hat{\mu}_{n}\geq 0$. We define the kernel complexity function as

$$\widehat{\mathcal{R}}(\delta)=\sqrt{\frac{1}{n}\sum_{j=1}^{n}\min\{\delta,\hat{\mu}_{j}\}}.$$

The critical radius is defined to be the smallest positive solution of $\delta_{n}>0$ to the inequality

$$\widehat{\mathcal{R}}(\delta)\leq\delta^{2}/\sigma.$$

Note that the existence and uniqueness of this critical radius is guaranteed for any kernel class. Based on this, we define the statistical dimension of the kernel is

$$d_{n}:=\min\{j\in[n]:\hat{\mu}_{j}\leq\delta^{2}_{n}\}.$$

Recall that, Theorem 2 in Yang et al. (2017) shows that various forms of randomized sketches can achieve the minimax rate using a sketch dimension proportional to the statistical dimension $d_{n}$. In particular, for Gaussian sketches and ROS sketches, the sketch dimension $m$ is required satisfy a lower bound of the form

$$m\geq\begin{cases}cd_{n}&\hbox{for Gaussian sketches},\\ cd_{n}\log^{4}(n)&\hbox{for ROS sketches}.\end{cases}$$

Here $c$ is some constant. In this paper, we adopt this specified sketch dimension $m$ to implement our experiments.