Ensemble Projection Pursuit  for General Nonparametric Regression

Ensemble Projection Pursuit

Haoran<EMAIL_ADDRESS>Mingke<EMAIL_ADDRESS>Yingcun<EMAIL_ADDRESS>
Department of Statistics and Data Science, National University of Singapore

The projection pursuit regression (PPR) has played an important
role in the development of statistics and machine learning.
According to the two cultures of [8], PPR is an
algorithmic model that can be used to approximate any
general regression. Although PPR can achieve the almost optimal
consistency rate asymptotically as shown in this paper, its effectiveness in prediction is
rarely seen in practice. To improve the prediction, we propose an ensemble procedure, hereafter referred to as ePPR, by adopting the “feature bagging” of the
Random Forest (RF). In comparison, ePPR has several advantages over RF, and its theoretical consistency can be proved under more general settings than RF. Extensive comparisons based on real data sets show that ePPR is significantly more efficient in regression and classification than RF and other competitors.

[Primary ]	62G08
[; secondary ]62G20

algorithmic model

consistency of estimation
greedy algorithm
nonparametric smoothing
projection pursuit regression
random forest

§ INTRODUCTION

Projection pursuit techniques were first introduced by [41] and [27]. The idea is to use lower-dimensional projections of high-dimensional data to reveal consecutively the most significant structures in the data. After a projection is found, the structures along that projection are removed, and the process is repeated to find new projections.
[26] implemented this idea for regression and called it the projection pursuit regression (PPR), and popularized the techniques in the 1980s; see also [36].  PPR is an algorithm which is able to approximate any general functions [38, 45, 46], and thus can be used to estimate general nonparametric regressions.

More specifically, consider the general regression problem with response $ Y $ and predictor $ X \in \mathbb{R}^p $,

and the conditional mean  $ m(X) = \E(Y|X) $ is our interest.

Following the categorization of two cultures of [8], the projection pursuit regression can be used in two ways. The first way of using PPR is to assume
$  m(X) =  g_1(\theta_1^\top X) + ... + g_K(\theta_K^\top X),
where $ K $ is assumed to be fixed and finite but unknown. See, for example, [15],
[31] and  [53].  Obviously, these statistical models may not be correct for real data and are not recommended by [8]. The second way of using PPR is to approximate or estimate a general nonparametric regression based on the following fact: for any smooth function $ m $, there exist a sequence of projection directions $ \{\theta_k \in \mathbb{R}^p: k=1, 2,...\} $ and a sequence of ridge functions $\{g_k(.), k = 1, 2, ...\} $ such that

\begin{equation}
\mbox{PPR expansion:} \ \ \  m(X) =  \sum_{k=1}^\infty g_k(\theta_k^\top X).  \label{PPRidentity}
\end{equation}

Under mild conditions, [45] and [46] proved that the PPR identity holds  for general function $ m(x)$ in the $L^2$ sense.  Actually, the PPR algorithm  of [26] was proposed  for the approximation and can be stated as follows. Let

\begin{equation}
\theta_1 = arg \min_{ \theta: ||\theta||_2=1} \E(Y-\E(Y|X^\top \theta))^2 \label{ippr.1}
\end{equation}

and $ g_1(v)  = \E(Y|X^\top \theta_1 = v) $. For $ k \ge 2 $, let $ \delta_k = Y- \sum_{\iota=1}^{k-1} g_\iota(X^\top \theta_\iota) $, and  calculate

\begin{equation}
\theta_k = arg \min_{ \theta: ||\theta||_2=1} \E( \delta_k -\E(\delta_k |X^\top \theta))^2  \label{ippr.2}
\end{equation}

and $ g_k(v)  = \E(Y|X^\top \theta_k = v) $. [37] proved that the PPR algorithm converges in the sense that

\begin{equation}\label{ippr.3}
\E\Big[ \Big|m(X) - \sum_{k=1}^{K} g_k(X^\top \theta_k)\Big|^2\Big] \to 0, \ \  as \ K \to \infty.
\end{equation}

[38] further accelerated the convergence rate of the PPR approximation  by using a relaxed greedy algorithm.

Although PPR has good algorithmic and mathematical guarantees, it has not been favoured by many nor has it been extensively studied. Existing computer packages of PPR do not provide good prediction accuracy in comparing with the  other popular methods such as the Random Forest  <cit.> and the support vector machine  <cit.>. [32] attributed the low popularity of PPR to its computational complexity in the 1980s when computational power was much lower than it is today, while later on, more attention was paid to its replacement, the artificial neural network (ANN). In fact, PPR should be studied in more depth,  especially after the tremendous development of new algorithms and computational power, which is the motivation for this paper.

In this paper, we first show that PPR can indeed achieve the optimal consistency rate under the general nonparametric setting of [50]. Then, we  propose a boosting method to improve the prediction accuracy of PPR, called ensemble PPR (ePPR) hereafter. The basic idea is to search for projections in (<ref>) and (<ref>) in a randomly selected subset of $ X $, in a similar way as that RF selects split variables [9], which is also called “feature bagging” [34]. Thus, each run of the algorithm produces a random PPR approximation of the unknown regression function. The ensemble technique is then used to generate a final estimate of the regression function.  Theoretically, we will give a proof to the consistency of the proposed ePPR under some regularity conditions.

§.§ Related work
[4] studied the consistency of two algorithms for general function approximation: the orthogonal greedy algorithm (OGA) and the relaxed greedy algorithm (RGA). Their results can be applied to the PPR models
m(x)= \sum_{\iota=1}^L g_\iota(\theta_\iota^\top x), x\in \mathbb{B}^p
under the assumption that $\sum_{\iota=1}^L \|g_\iota\|_{H^r(\B^p)}<\infty$, where $ L $ can be either finite or infinite, i.e., $ L \le \infty$. When $L$ is finite, [6] further obtained the optimal consistency rate using a one-layer neural network. Note that the assumption is strong and cannot be satisfied generally by $ m(x) $ in the Sobolev space. In comparison, we will show in this paper that the PPR algorithm can achieve the optimal consistency rate up to a multiplier of $ \ln^4 n $ for a more general $m$, where $ n $ is the sample size.

As the main contribution of this paper, our proposed ensemble PPR (ePPR) is partially motivated by RF of [8]; see also [55] and [2] for the recent development of RF. ePPR shares one of the core elements of RF: the former uses recursive projection and the latter recursive partitioning. ePPR further borrows two of the other core elements of RF: random selection of predictors and averaging estimates from each member of the ensemble. The details will be illustrated below in Section <ref>. As noticed by [23], RF has weak ability to fit smooth signals because of the step-functions inherited from the classification and regression tree (CART). Obviously, PPR does not have this problem, and neither does our ePPR. [55] also considered linear combinations of the predictors as the splitting variables of RF, but its way in selecting the combinations is different.  Another problem with RF is that theoretically it only works for some special classes of regression functions; see, for example, [48] and [17]. Although [7] suggested using linear combinations of predictors as splitting variables in CART, but its extension to Random Forest does not see much improvement. The method became known as the oblique decision tree (ODT) and has received lots of attention.
In comparison, PPR or ePPR works for more general regression functions, and their statistical consistency is guaranteed.  Our extensive numerical study in Section <ref> also lends support to the superior performance of ePPR over RF and other competitors for both classification and regression problems.

The rest of this paper is organized as follows. Section <ref> introduces some notations and studies asymptotic efficiency of a full estimation algorithm for a general regression function based on the PPR identity.  Section <ref> discusses two existing greedy algorithms and a new one based on randomly selected variables, for the preparation of our main algorithm, ePPR. Details of ePPR are presented in  Section <ref>, statistical consistency of our estimation is also studied in the section. Extensive numerical studies of  real data sets, including comparison with the popular machine learning methods such as RF and SVM,  are presented in Section <ref>.  The technical proofs are provided in Section <ref>, while some of the auxiliary results  used in the proofs are given in a supplementary material (SM) in a separate file.
We summarize the advantages as well as the problems with the proposed ePPR in Section <ref>.

§ ESTIMATING $\E(Y|X)$ BY PPR

In this section, we consider the estimation of a general regression
function $m(X)=\E(Y|X)$ by a linear combination of ridge functions, i.e. PPR model, using
i.i.d. sample $\mathcal{D}_n = \{(X_i, Y_i), i = 1, ..., n\} $. Denote by $ \mu(x) $ the cumulative distribution function of $ X $. Let
$\B^p $ be the unit ball in $\mathbb{R}^p$ centered at the origin
and $L^2(\B^p)$ be the functional space on $\B^p$ with finite
$L^2$-norm. The function space which $m$ lies in is assumed to be a
Sobolev space $H^{r}(\B^p)$,
$$ H^{r}(\B^p) =\{u(\cdot)\in
L^2(\B^p): D^\kappa u\in L^2(\B^p), \ \text{for}\ 0\leq
|\kappa|\leq r \}, $$
where $r>0$ is an integer and
$|\kappa|=\sum_{i=1}^p\kappa_i$ for any tuple
$\kappa=(\kappa_1,\ldots,\kappa_p)\in\mathbb{Z}_+^{p}$, where $ \mathbb{Z}_+ $ is the collection of non-negative integers. The
derivative $D^\kappa m$ is in the weak sense and when $m\in
H^{r}(\B^p)$, the Sobolev norm of $m$ is denoted by
$\|m\|_{H^{r}(\B^p)}$; see [1] or Section <ref> of SM
for more details.

[45] showed that any smooth multivariate function $m \in H^{r}(\B^p)$ can be approximated by a linear combination of functions in  $L^2([-1,1])$. This result can be briefed as follows. Let $V_n$ be a linear space of  univariate functions in $L^2([-1,1])$ with approximation order $ s $, satisfying

\begin{equation}\label{sapproximation}
\inf_{v\in V_n}\|v- g\|_{L^2([-1,1])}\leq c(s, \|g\|_{H^{s}([-1,1])})\cdot J_n^{-s}
\end{equation}

for any $g\in H^{s}([-1,1])$, where $ H^{s}([-1,1]) $ is defined in the same way as $ H^r(\B^p) $ above, and $c(s, ||g||_{H^{s}([-1,1])})$ only depends on integer $ s $ and the Sobolev norm of $g$,  and $J_n$ denotes the dimension of $V_n$. Condition (<ref>) holds for many commonly used approximation spaces, such as splines defined in (<ref>) and the algebraic polynomials $V_n=\text{span}\{1,x,\ldots,x^n\}$; see [21]. [45] constructed a series of  finite subsets of the unit sphere $\Theta^{p} = \{ \theta: \theta \in \mathbb{R}^{p}, ||\theta||_2=1\}$, namely $\Theta_n\subseteq \Theta^{p}$, $n=1,2,\ldots$ with cardinality $K_n = card(\Theta_n) \asymp J_n^{p-1}$,

and proved that for any $m\in H^\lambda(\B^p)$, there is a function $\tilde m_n\in  \text{span}\{ v(x^\top \theta): v\in V_n,\ \theta \in\Theta_n\}$, i.e. $ \tilde m_n = \sum_{j=1}^{K_n}  \beta_j v_j( x^\top \theta_j) $ with $ v_j \in V_n $,  such that

||m-\tilde m_n ||_{L^2(\mathbb{B}^p)} \leq c(\lambda, ||m||_{H^{\lambda}(\B^p)})\cdot J_n^{-\lambda}, % \ \forall G\in H^{\lambda}(B^p)

where $\lambda= s+(p-1)/2$, and $c(\lambda, ||m||_{H^{\lambda}(\B^p)})$ only depends on $ \lambda $ and the Sobolev norm of $m$.

To apply the above result of [45] to PPR, we first approximate the ridge function $ g_k $ in (<ref>) by
$$g_k^*(\cdot)= \sum_{j=1}^{J_n} \alpha_{k,j} v_j(\cdot)$$

where $\{v_{j}\}_{j=1}^{J_n}$ is a basis of $V_n$. Minimize the following empirical loss function

\begin{equation}\label{Theoreticalestimator}
(\hat{\alpha}_{k,j},\hat{\theta}_k)  = \argmin_{\alpha_{k,j}\in\mathbb{R},\theta_i\in \Theta^{p}} \frac{1}{n}\sum_{i=1}^n \Big( Y_i-  \sum_{k=1}^{K_n} g_k^*(\theta_k^\top X_i) \Big)^2,
\end{equation}

and  estimate $m$ by

$$m_n(x)=\sum_{k=1}^{K_n}\hat{g}_k^*(\hat{\theta}_k^\top x),\ x\in\mathbb{B}^p,$$

$\hat{g}_k^*(\cdot)=\sum_{j=1}^{J_n} \hat{\alpha}_{k,j} v_j(\cdot)$.
Because the estimation (<ref>) is based on minimization with respect to  all the functions and parameters, we call it full algorithm of PPR.
To study the asymptotic property, we define a truncated version of $m_n$ by

\begin{equation}\label{trunc}
\hat{m}_n(x)=  \max\{\min\{m_{n}(x),t_n\},-t_n\}
%	\begin{array}{ll}
%	t_n, &  \ if \ \hat{m}_n>t_n,\\
%	\hat{m}_n, &  \ if \ -t_n\leq \hat{m}_n\leq t_n,\\
%	-t_n, &  \ if  \ \hat{m}_n<-t_n\\
%	\end{array}
%	\right
\end{equation}

for some $t_n\in\mathbb{R}^+$ which will be specified later.  The consistency of $m_n$ is as follows.

Assume $\E(e^{c\cdot Y^2})<\infty$ for a constant $c>0$ and that the probability density function  of $ X $ is bounded. Suppose $ V_n $ satisfies (<ref>) and each
$v_j$ consists of a finite number of pieces of polynomials with degrees not larger than $s$. If $J_n\asymp \left(\frac{n}{\ln^4 n}\right)^{\frac{1}{2r+p}}$ and $t_n\asymp \ln n$, then the truncated estimate $\hat{m}_n$ in (<ref>) satisfies

\begin{equation}
\E \int|\hat{m}_n(x)-m(x)|^2 d\mu(x)= O(\left(\frac{\ln^4 n}{n}\right)^{\frac{2r}{2r+p}})
\end{equation}

for any bounded $m\in H^r(\B^p)$ with  $ 0<r\leq s+\frac{p-1}{2}$, where $ s $ is the approximation order of $ V_n $, and $ \mu(x) $ is the cumulative distribution function of $ X $.

showed that the optimal minimax consistency rate for estimators of functions in $ H^r(\mathbb{B}^p)$ is $O(n^{-\frac{2r}{2r+p}})$; see also [50]. Therefore,  the consistency rate of our  estimator (<ref>) is  optimal up to a multiplier factor of $\ln^4 n$.  When $r\in (s,s+\frac{p-1}{2})$, the consistency rate of $\hat{m}_{n}$ is still almost optimal, while the rate of the estimator based on tensor product B-splines of order $s$ is at most $O(n^{\frac{2\cdot(s+1)}{2\cdot(s+1)+p}})$ (see, for example, [30]). Thus, Theorem <ref> shows an advantage of using decomposition of ridge functions; and such unexpected gain, which was also observed by [45], is because of the  approximation theory; see also   [45].

Next, we illustrate the implementation of  (<ref>) by using splines which include the popular rectified linear activation function (ReLU).

For each $n\geq 1$, define a mesh $\mathcal{T}_n: u_{-s}\leq u_{-s+1}\leq \cdots\leq u_{-1}=u_0=-1\leq u_1\leq\cdots \leq u_{J_n-1}=1\leq u_{J_n}\leq \cdots \leq u_{J_n+s} $ of the interval $[-1,1]$.   A commonly used assumption on the meshes is

\begin{equation}
\frac{\max_{0\leq i\leq J_n-2}\{u_{i+1}-u_i\}}{\min_{0\leq i\leq J_n-2}\{u_{i+1}-u_i\}}\leq c,  \label{mesh assumption}
\end{equation}

where $ c $ is a universal constant. We will assume (<ref>) throughout the rest of this paper.

Let $ \{B_{i, s}, i = 1, \ldots, J_n \}$ be the set of B-splines basis functions of order $ s $ based on the meshes above. With these basis functions, $V_n$ is realized by the following spline space

\begin{equation}\label{splinespace}
{\cal S}_{u,s}[-1,1]=\text{span}\{ B_{1, s},\ldots,  B_{J_n, s}\}.
\end{equation}

Note that  the order of spline $ s$ is also  the approximation order of $ V_n $ defined in (<ref>). Replacing each $ v_j $  by a B-spline, i.e. $ v_j = B_{j, s} $, $j=1,\ldots, J_n$,  minimizing (<ref>) will give the estimator of  $ m_{n}(x) $ and its truncated version, $ \hat m_n(x) $.

Let $ V_n =  {\cal S}_{u,s}[-1,1] $.  If the mesh is designed according to (<ref>), then $ V_n $ satisfies (<ref>), and thus the result of Theorem <ref> follows.

Under the condition that $\E(Y|X)$ follows a projection pursuit regression model with finite terms, i.e. $ \E(Y|X) = \sum_{k=1}^K g_k(\theta_k^\top X) $ for some finite $ K$, [6] gave a similar result as Theorem <ref> by using one hidden layer shallow neural network. Indeed, PPR can be regarded as a generalization of the shallow neural network. [39] showed that deep neural networks with finite neurons in each layer is dense in the space of continuous functions defined on a compact set in $\mathbb{R}^p$.  Following this idea, extending PPR to multi-layer is possible, which, however, is not our interest of this paper.

Note that the minimization of (<ref>)
involves a big number of parameters, its calculation is very unstable. [51] implemented the calculation by imposing penalties on the parameters. Although some improvement in prediction over the original PPR has been made by  [51], it still cannot compete with RF. More efficient algorithms need to be developed.

§ PPR APPROXIMATION USING RANDOMLY SELECTED VARIABLES

Although the full algorithm of PPR has a desirable consistency rate for the estimators as stated in Theorem <ref>, its  numerical performance is not satisfactory compared with the other estimation methods, such as RF and SVM; see [51] and  Section <ref>. Thus, new algorithms are needed. Next,  we will introduce the orthogonal greedy algorithm (OGA) and the relaxed greedy algorithm (RGA) to PPR.
[37, 38] showed that a RGA algorithm of  PPR has a faster convergence rate than the original PPR.  [44] and [20] found that OGA can  prevent the algorithm from selecting atoms that have already been added in earlier iterations.

To boost the estimation efficiency, we further use the “feature bagging” of [34] as RF does. Note that at each knot, RF selects the split variable from a random subset of the predictors rather than the full set of predictors $ X $ as CART does [10]. In the same manner, for each projection of (<ref>) and (<ref>) of PPR, we will randomly select a subset of the predictors for the projection. Because of the random selection, each run of the resultant algorithm  will generate different estimation and prediction. We can then take average of those  predictions to make a final estimation or prediction. We call this procedure ensemble PPR, denoted by  ePPR.

In this section, we start with the implementation of the PPR algorithm with the projection in a randomly selected subset of $X $, called rPPR hereafter. After the projections are estimated, the PPR model has an additive structure of the projected variables. Motivated by RGA and OGA, we first present two algorithms, called  rPPR.oga and  rPPR.rga respectively, and then we use the additive structure to propose another rPPR algorithm, called additive greedy algorithm (rPPR.aga).

For ease of exposition, we introduce several notations that are often used in the rest of this paper. Denote by $\mathcal{A}$ a (random) subset of $ \{ 1, ..., p\} $,  and $\mathcal{A}_\tau,\ \tau=1, 2, ...$ is a sequence of such subsets that may differ from one another. Let $ X_\A $ be a sub-vector of $ X $ consisting of coordinates indexed by elements of $\mathcal{A}$. Similarly, $ x_\A $ can be defined for vector $ x $, and $ X_{i, \A} $ for a sample $ X_i $ of $ X $. To simplify the notation, by writing $ \theta^\top X_\A $ we imply that $ \theta $ has the same dimension as $ X_\A $.  Let $\mathbf{\Omega}_q$ be the collection of all $\A$s with $card(\A) = q$. We can also think of $\mathbf{\Omega}_q $ as a probability space with probability equally assigned to every element. Thus, $ \A_\tau, \tau\geq 1 $, can be regarded as a (random) sample of $ \mathbf{\Omega}_q $. Let $\Xi_k=(\A_1,\ldots, \A_k)$  with  $k\geq 1$, which can be regarded as a random element in the product of probability spaces, $ \mathbf{\Omega}_q ^{\otimes k} $.

We first describe the rPPR algorithms for the approximation of any function $ m(x) $. We assume tentatively that $ m(x) $ is known to study the theory of the algorithms on the population level, while sample versions of the algorithms will be presented in Section <ref>.

rPPR.rga.  Starting with iteration $ k = 0 $, let $ \delta_0 = m $ and $ m_0 = 0 $. For $k \ge 1$, select a random subset $ \A_k \in \mathbf{\Omega_q} $ and

\begin{equation}
\theta_k = arg \min_{ \theta\in \Theta^q} \E( \delta_{k-1} -\E(\delta_{k-1} |X_{\A_k}^\top \theta))^2  \label{ppr.2}
\end{equation}

and $ m_k  = \alpha_k \cdot m_{k-1} + \E(\delta_{k-1}| X_{\A_k}^\top \theta_k) $, where $\delta_{k-1}=m-\alpha_k\cdot m_{k-1}$ and $ 1> \alpha_k > 0 $ is an increasing sequence, called relaxing coefficients.

rPPR.oga.  Let $ \delta_0 = m $ and $ m_0 = 0 $. For $k \ge 1$, randomly select $ \A_k $ and search for the projection direction $ \theta_k $ by  (<ref>)
with  $\delta_{k-1}=m-m_{k-1}$. Let $ g_k(v) = \E(\delta_{k-1}|X_{\A_k}^\top \theta_k = v) $ and update
$ m_k  = \displaystyle\sum_{\tau=1}^k c_\tau {g_\tau(X_{\A_\tau}^\top \theta_\tau)} $ where coefficients $ c_\tau $s are calculated by

(c_1, ..., c_k) = \argmin_{c_1,\ldots,c_k\in\mathbb{R}} \E\{( m(X) - c_1 \cdot g_{1}(X_{\A_1}^\top \theta_1) - ... - c_k\cdot g_k(X_{\A_k}^\top \theta_k))^2\}.

rPPR.aga.  Let $ \delta_0 = m $ and $ m_0 = 0 $. For $k \ge 1$, select $ \A_k $ and search for the projection direction $ \theta_k $ in the same way as (<ref>) with $\delta_{k-1}=m-m_{k-1}$. Let $ Z_k = X_{\A_k}^\top \theta_k $. Refit functions $ g_1, ..., g_k $ by

(g_1, \ldots , g_k) = \argmin_{g_\tau : \E g_\tau^2(Z_\tau)<\infty, \tau=1,\ldots, k }{\E|m-g_1(Z_1)-\ldots-g_k(Z_k)|^2}

and let $m_k(X) = g_1(Z_1) + ... + g_k(Z_k) $.

Note that in the above algorithms, $ \A_k $ is chosen at random, thus $ X_{\A_k} $ may be completely independent of $m(X)$. As a consequence, the projection with such $ X_{\A_k} $  makes no contribution to the approximation. To improve the efficiency of the algorithms, we can consider a few random candidates $ \A_{k,1}, ..., \A_{k,k'} $, where $ k' \ge 1 $ is pre-specified, and select $ \A_k $ that corresponds to the smallest of $ \min_{ \theta\in \Theta^q} \E( \delta_{k-1} -\E(\delta_{k-1} |X_{\A_{k,\iota}}^\top \theta))^2, \iota = 1, ..., k' $; see Algorithm <ref> for details.

Minimization (<ref>) can be implemented by the estimation method for the single-index model, such as [26] and  [31]. We shall discuss this later in more details for the sample versions of the algorithms. Note also that we do not mention the ensemble in the algorithms about the rPPRs so far, because this is an population level and the ensemble is equivalent to taking expectation over  $ \Xi_k = (\A_1, ..., \A_k) $, i.e. $ m_k^{ens}(x) = \E_{\Xi_k} m_k(x) $. We shall discuss the ensemble in detail for sample versions of the algorithms.

In the algorithms, the number of variables used in the projections, $ q$, plays an essential role as that used in RF where $ q = p/3 $ (or $\sqrt p$) is used [9]. Choice of $ q $ determines the diversity of individual estimates from different runs of the algorithms. Empirically, ensembles tend to yield better result when there is a significant diversity among the individual estimates; see [42] and [12]. Many ensemble methods, therefore, seek to promote diversity among the individual models they combine, which requires a relatively small $ q $ in our approach. However, a too small $ q $ may also lead to inefficient approximation of the function $ m(x) $  in each run.  Similarly, if $ \A_k $ is selected from $ k'$ candidates as discussed above, then larger $ k' $ will produce less diversity. As two extreme cases, $k'=1$ corresponds to the algorithms described above, and $ k' = {p \choose q} $ gives the same output for all runs of teh algorithm and has no diversity. More discussion can be found later. The relaxing coefficients  $ \alpha_k $ in rPPR.rga is an increasing sequence, a simple choice of it is $ \alpha_k = 1 - 2/(k+2) $ for any integer $k\geq 1$; see  [4] for more discussions about the relaxing coefficients.

For the three algorithms, our simulations suggest that ePPR based on rPPR.aga has slightly better predictions than those based on rPPR.rga and rPPR.oga, although their theoretical convergence rates are similar as shown below.

§.§ Approximation rates

In this section, we are interested in how the difference between   $ m_k(x) $  and  $ m(x) $ changes with  $ k $, i.e. convergence rates of the algorithms.

Because the number of variables $q $ selected for the projection and  “feature boosting”, the approximation space is thus

\begin{equation}\label{Hspace}
\mathcal{H}_{q}=\left\{ \sum_{\tau=1}^\infty \beta_{\tau} v_\tau(\theta_{\tau}^\top  x_{\C_\tau}) : \beta_{\tau}\in\mathbb{R},\ \theta_\tau\in \Theta^{p}, \ \C_\tau \in \mathbf{\Omega}_{q}   \right\}
\end{equation}

where $v_\tau\in V_{\theta_\tau, \C_\tau}=\{ v_\tau: \E (v_\tau^2(\theta_{\tau}^\top  X_{\C_\tau}))=1\}$, $\theta_\tau\in\mathbb{R}^{q}$, and $ x_{\C_\tau} $ is defined in the same way as  $ x_{\A_\tau} $, but $ {\C_\tau}$ is nonrandom and unknown. A function in $\mathcal{H}_{q}$ may have finite terms $ t $, i.e. $\beta_\tau=0$ for all $\tau>t$.
Note that outputs of the above three algorithms belong to $ \mathcal{H}_{q}$.  For any function $h(x)=\sum_{\tau=1}^\infty \beta_{\tau} v_\tau(\theta_{\tau}^\top  x_{\C_\tau}) \in \mathcal{H}_{ q} $, define its norm by
\|h\|_{EL^1}= \sum_{\tau=1}^\infty |\beta_{\tau}|.

For any function $m(x)$, not necessarily in $\mathcal{H}_{q} $,  the outputs of algorithms  rPPR.aga and rPPR.oga satisfy

\E_{\Xi_k, X}(m_k(X)-m(X))^2\leq \inf_{h\in \mathcal{H}_{q}}\{\E(m(X)-h(X))^2 + 4\left(c(p,q) \|h\|_{EL^1}\right)^2\cdot(k+1)^{-1}\},

for all $k= 1, 2,\ldots$ where $c(p,q)={\binom{p}{q}}$. For algorithm rPPR.rga, if the  relaxing coefficients $\alpha_k=1-2/(k+2),\ k \ge 1$, then the above equation also holds.

discussed the selection of relaxing coefficients, which takes the following form

\begin{equation}\label{relaxingcoefficient}
\alpha_k=\max\bigg\{ \left(1-\frac{\lambda}{k}\right),0 \bigg\}, \ k=1,2,\ldots,
\end{equation}

where $\lambda>1$. [4] argued that when $\lambda=2$, the upper bound of the approximation error given in their paper achieves the minimum. This is the reason why we choose $\lambda=2$ in Theorem <ref>. We also modify the relaxing coefficients $\alpha_k$ in (<ref>) in order to make them  positive.

Theorem <ref> indicates  that the approximation error depends on two parts. The first part is the distance between the target function $ m $ and the estimation space $ \mathcal{H}_{q} $. This is commonly the case for serial estimation methods including the splines: the bigger the space the better approximation to the target function. If $m(x) $ belongs to some specific function space,  for example $ m(x)  \in H^r(\B^q) $ or $ m $ is a sum of functions with each having $q$ variables as in (<ref>) and (<ref>) below, then the distance is $ O(card({\cal S}_{u,s}[-1,1])^{-r})$ under some mild conditions; see [45]. In those cases, Theorem <ref> implies

\lim_{k\to\infty}{\E_{\Xi_k, X}(m_k(X)-m(X))^2} =0,

which indicates that rPPR.aga, rPPR.oga and rPPR.rga produce estimators that converge with iterations to the true function $ m $.  However, for the sample versions of the algorithms, on the other hand, the larger $ \mathcal{H}_{q} $ is, the more unstable the estimator is.  The second part of the above approximation error depends on the choice of $ q $ and  the number of iterations. The choice of $ q $ will determine the diversity of the individual runs of the algorithm, and the term $ c(p, q) $ achieves the maximum when $ q = p/2 $. However, as we will see later, many other factors need to be taken into consideration when we choose $ q $ in the sample versions of the algorithms; see also [40],[48] and  [28].
For the number of iterations, the term $ (k+1)^{-1} $ means that more iterations will theoretically produce more accurate approximation. Note that this is also not true for sample version, which will be discussed in next section.

As mentioned before, the ensemble of outputs of each of the algorithms is actually the average the outputs over all possible choices of $ \Xi_k $, i.e. $ m_k^{ens}(X) = \E_{\Xi_k} m_k(X) $. By Jensen's inequality, we have

\E_X\{ (m_k^{ens}(X) - m(X))^2\}  \le \E_{\Xi_k, X}(m_k(X)-m(X))^2.

Thus,  it follows from Theorem <ref> that the ensemble of the outputs of each algorithm also converges and is possible to be more accurate than the individual approximation, mainly through the reduction of variance; see, for example, <cit.>.

§ ENSEMBLE PPR ESTIMATORS

Now,  consider the implementation of rPPR.aga, rPPR.rga and rPPR.oga based on observed data $ \mathcal{D}_n=\{(X_i, Y_i), i = 1, ..., n\} $, and their ensembles, called ePPR.aga, ePPR.oga and ePPR.rga respectively.  Here, we only give the details for ePPR.aga following rPPR.aga; the other two, ePPR.oga and ePPR.rga, can be implemented similarly. We use the spline space $ {\cal S}_{u,s}[-1,1]=\text{span}\{ B_{1, s},\ldots,  B_{J_n, s}\}$ defined in (<ref>) to approximate each unknown ridge functions $g_\iota, 1\leq \iota \leq k,$ in the calculation.

ePPR.aga (ensemble PPR via additive greedy algorithm)

Data $(X_i,Y_i), i=1, ..., n $; the size of each candidate subset, $q$ ($q\leq p$), for the projection; and the number of candidate subsets, $ \ell$; and maximum number of ridge functions $k$; and number of boosts $B$.

$b \in  \{1, 2, ..., B\}$

$\ e_i \leftarrow Y_i  \ \text{for}\ i=1,\ldots,n$

$\tau \in  \{1, 2, ..., k\}$

$ l \in  \{1, 2, ..., \ell \}$

randomly choose a subset $ \A_\tau^l \subset \{1, ..., p\} $ with $card(\A_\tau^l)=q$, calculate
\begin{aligned}
R^l(\theta^l_\tau, \beta^l_1, ..., \beta^l_{_{J_n}})=\min_{\theta\in \Theta^{p}, \beta_1, ..., \beta_{_{J_n}}\in\mathbb{R}}{\sum_{i=1}^n{\Big( e_i-\sum_{j=1}^{J_n}{\beta_{j}\cdot B_{j,s}(\theta^\top  X_{i, \A^l_\tau} )}\Big)^2}};
\end{aligned}

\theta_\tau \leftarrow \theta^{l_0}_\tau $ where $ l_0 = \argmin_{1 \le l \le \ell}R^l(\theta^l_\tau, \beta^l_1, ..., \beta^l_{_{J_n}})  $

$Z_{i,\tau} \leftarrow \theta_\tau^\top X_{i, \A_\tau}, \ \ i=1,\ldots,n$

\begin{aligned}
(\beta_{\iota,j}, \iota = 1, ..., \tau, j = 1, ..., J_n)  =
\argmin_{\beta_{\iota,j}\in\mathbb{R}} \  \sum_{i=1}^n \{Y_i-\sum_{\iota=1}^\tau \sum_{j=1}^{J_n} \beta_{\iota,j}\cdot B_{j,s}( Z_{i,\iota} )\}^2;
\end{aligned}

\begin{aligned}
g_\iota(v) \leftarrow  \sum_{j=1}^{J_n} \beta_{\iota,j}\cdot B_{j,s}(v), \ \ \iota=1,\ldots,\tau;
\end{aligned}

$m_{b,\tau}(x)\leftarrow  g_1(\theta_1^\top x_{\A_1})+\ldots+ g_\tau(\theta_p^\top x_{\A_\tau})$

$e_i \leftarrow  Y_i- m_{b,\tau}(X_i)$

estimator of $ m $: $m_{k}^{ens} = B^{-1} \sum_{b=1}^B m_{b,k}$.

Details of ePPR.aga algorithm are given in Algorithm <ref>.
In line 5 of the algorithm, estimating $ \theta_\tau $ is equivalent to the estimation of a single-index model, the algorithm of [26] can be used. See also  [31] and [35]. Unlike  (<ref>) where the projection directions $ \theta_1, ..., \theta_k $ are estimated simultaneously, Algorithm <ref> estimates the directions   in a sequential manner. Line 9 of the algorithm is to estimate an additive model, for which many other existing methods can be used; see, for example, [33] and [25].

There are four tuning parameters in the algorithm: the number of variables, $ q $, in each random subset for feature boosting, number of candidate subsets of variables, $ \ell $,   and number of iterations (or ridges),  $k$, of each individual run of the algorithm, and the number of runs, $ B $, for the ensemble. We have discussed the selection of $ q $ for rPPR.aga, but for sample version, we also need to consider the estimation efficiency of the projection $ \theta_\tau $. Because of the nonparametric estimation of the projections, $ q $ should be relatively much smaller than the sample size $ n $ to get a good estimator of the projection; see [31]. Thus, we select $ q = \min([2p/3], [n^{0.4}]) $ in our calculation.  For a similar reason, the choice of $ \ell $ also needs to consider the diversity of the estimators $  m_{1,k}, ...,  m_{B,k}$ and the efficiency of each individual $ m_{b,k} $. As an extreme case, if $ \ell = {q \choose p}$ then all the estimators are the same with no diversity. Our limited numerical study suggests that $ \ell \le p/q $ can be used. The selection of $ k $ depends on the model structure, and will be discussed in the next section. For $ B$, there are studies about its selection for RF; see, for example, [54]. However, we follow the common practice and take $ B= 500$ in our calculation.

§.§ Consistency of ePPR estimators
Due to technical complexity, we only study the asymptotic behaviors of  $\tilde{m}^{ens}$ for two models. The first model is the extended additive model:

\begin{equation}\label{Extendedmodel1}
\end{equation}

where each $m_\tau$ is a function of only $q_0$ variables $x_{\C_\tau}$ with $ card(\C_\tau) = q_0 $; and the second model is a PPR model:

\begin{equation}\label{Pprmode1}
m(x)=\sum_{\tau=1}^\infty {g_\tau(\theta_\tau^\top  x_{\C_\tau})},
\end{equation}

where $\theta_\tau\in\mathbb{R}^{q_0}$ and $ card(\C_\tau )= q_0 $. Note that the maximum terms $ L $ in model (<ref>) is $ q_0 \choose p $, and that (<ref>) is more general than (<ref>). When $ q_0 = 1$,  the two models above degenerate to an additive model. For convenience of derivation, the theoretical results are only given for the truncated version of ${m}^{ens}_{k_n}$,

\hat{m}^{ens}_{k_n}=\frac{1}{B}\sum_{b=1}^{B} {\hat{m}_{k_n,b}},

where $\hat{m}_{k_n, b}= \max\{\min\{{m}_{k_n,b},t_n\},-t_n\}, j=1,\ldots, B$, and $t_n$ is the threshold, and $t_n\to\infty$, while in the calculation this truncation is not necessary.

For model (<ref>), assume that the distribution of $(X,Y)$ satisfies
\E(e^{c\cdot Y^2})<\infty,
where the constant $c>0$ and that the joint density of $X$ is bounded. Suppose $m$ is bounded and $m_\tau\in H^r(\B^{q_0}), 1\leq \tau \leq L,$ where $r\leq s+(p-1)/2$, and $ s $ is the approximation order of $ S_{u, s}[-1, 1] $. If we select the meshes in an increasing order, $\mathcal{T}_1\subset \mathcal{T}_2\subset\cdots\subset \mathcal{T}_n \subset\cdots$, and $t_n\asymp \ln n$ and $k_n,J_n\to\infty$, then

\E_{\mathcal{D}_{n},\Xi_{k_n}} \int|\hat{m}^{ens}_{k_n}(x)-m(x)|^2d\mu(x)\to 0

as the sample size $n\to\infty$.

For  model (<ref>), assume that the distribution of $(X,Y)$ satisfies
\E(e^{c\cdot Y^2})<\infty,
where the constant $c>0$. Suppose $m$ is bounded and that $g_\tau\in H^r([-1,1]), 1\leq \tau <\infty $, satisfy $\displaystyle\sum_{\tau=1}^\infty\|g_\tau\|_{H^r([-1,1])}<\infty$ where $r\leq s+(p-1)/2$.
When $t_n\asymp \ln n$ , $k_n\asymp n^{\frac{2r}{4r+1}}$ and $J_n\asymp n^{\frac{1}{4r+1}}$, we have

\E_{\mathcal{D}_{n},\Xi_{k_n}} \int|\hat{m}^{ens}_{k_n}(x)-m(x)|^2d\mu(x) = O(
n^{-\frac{2r}{4r+1}}\cdot \ln^4 n).

Note that the consistency rate in Theorem <ref> is not optimal [50]. This is mainly because  ePPR is a greedy algorithm, which does not have very fast consistency rate [38]. Nevertheless, these results still reveal some attractive properties of $\hat{m}^{ens}_{k_n}$ that have not been proved for RF. Actually, it is known that RF is consistent under the model (<ref>) when $q_0=1$, i.e. the additive model [48]. However, to our best knowledge,  only [40]  showed that the consistency rate of RF is  $O_p\left((\ln n)^{-1}\right)$ which is much slower than Theorem <ref>. Therefore, ePPR is more efficient than RF in the minimax sense and its statistical consistency is guaranteed for more general regressions.

§.§ Stopping rule for the algorithms

In practice, we still need a criterion to stop the algorithm at the right time  of the iterations  to balance the estimation bias and over-fitting.
Next, we propose  a BIC criterion for this purpose. Consider the $b$-th run of the algorithm, denote by $m_{\tau,b}$ the output function of the inner-loop of ePPR.aga (or  ePPR.oga or ePPR.rga). Define a sequence of losses for each  individual estimator $m_{\tau,b},\tau=1,2,\ldots, $ by

\begin{equation}\label{BICeppr}
\text{BIC}_{b}(\tau)= n^{-1} \sum_{i=1}^n (Y_i -m_{\tau,b}(X_i))^2+\tau\cdot\frac{\ln{n}\cdot(q+J_n^{1+\nu})}{n},
\end{equation}

where $\nu>0$. In the BIC above, $\tau (q+J_n^{1+\nu}) $ can be regarded as the number of parameters in the model. We choose the stopping time $k^*_{b}$ for the $b$-th run by

\begin{equation}\label{BICchooseK}
k^*_{b}= \argmin_{1\leq \tau<\infty}\text{BIC}_{b}(\tau).
\end{equation}

In calculation, we select $k_b^*$ as the first $\tau\geq 1$ such that $BIC_b(\tau)<BIC_b(\tau+1)$. Finally, the ePPR estimator based on BIC criterion is given by

\begin{equation}\label{EppRBICestimator}
m^{ens}_{*}= \frac{1}{B}\sum_{b=1}^{B} {m_{k^*_{b},b}}.
\end{equation}

Again, to analyze asymptotic properties of estimator (<ref>), the truncated version  $m^{ens}_{*}$:
\hat{m}^{ens}_{*}=\frac{1}{B}\sum_{b=1}^{B} {\hat{m}_{k^*_{b},b}}
is considered where $\hat{m}_{k^*_{b},b}= \max\{\min\{m_{k^*_{b},b},t_n\},-t_n\}$ and $t_n>0$ is the threshold;  see (<ref>).

For  model (<ref>), assume that the distribution of $(X,Y)$ satisfies that
\E(e^{c\cdot Y^2})<\infty,
where the constant $c>0$ and that  the joint density of $X$  is bounded. Suppose  $m\in H^r(\B^p)$ is bounded where $r\leq s+(p-1)/2$ and $m_\tau\in H^r(\B^{q_0}), 1\leq \tau \leq L$. If we design the meshes in an increasing order: $\mathcal{T}_1\subset \mathcal{T}_2\subset\cdots\subset \mathcal{T}_n \subset\cdots$, and $t_n\asymp \ln n$ and $J_n \asymp n^\frac{1}{1+\nu}$, then

\int|\hat{m}^{ens}_{*}(x)-m(x)|^2d\mu(x) \xrightarrow{p} 0

as the sample size $n\to\infty$.

Note that Theorem <ref>  only requires that $k_n,J_n\to\infty$ as $n\to\infty$, while $J_n=o(n^\frac{1}{1+\nu})$ is required in Theorem <ref>. This is because $\hat{m}^{ens}_{*}$ is obtained by imposing a penalty in (<ref>) that  depends on $\nu$.

For  model (<ref>), assume that the distribution of $(X,Y)$ satisfies that
\E(e^{c\cdot Y^2})<\infty,
where the constant $c>0$ and that the joint density of $X$ is bounded. Suppose $m$ is bounded and that $g_\tau\in H^r([-1,1]), 1\leq \tau < \infty $, satisfy $\displaystyle\sum_{\tau=1}^\infty\|g_\tau\|_{H^r([-1,1])}<\infty$ where $r\leq s+(p-1)/2$.
When $t_n\asymp \ln n$, $J_n\asymp n^{\frac{1}{4r'+1+\nu}}$ and $r'=\min\{r,\nu/2\}$, we have

\int|\hat{m}^{ens}_{*}(x)-m(x)|^2d\mu(x)=O_p\left( n^{-\frac{2r'}{4r'+1+\nu}}\cdot \ln^4 n\right).

This theorem reveals that $\hat{m}^{ens}_{*}$ selected by BIC criterion has a similar consistency rate as $\hat{m}^{ens}_{k_n}$ in Theorem <ref>.
The difference is that the rate of $\hat{m}^{ens}_{k_n}$ depends on the parameter $\nu$ due to the penalty defined in (<ref>). It suggests that for a given smoothness $ r $, $ \nu = 2r $ is the optimal choice leading to the highest consistency rate. In other words,  the smoother $m(x)$ is, the fewer inner iterations are required for the estimator to achieve a predefined error, implying that the larger $\nu$ needs to be used. The value of  $ \nu$ also has impact on the  consistency rate as shown in Theorem <ref> because the larger $ \nu $ allows the larger $r'$. Intuitively, using fewer ridge functions also helps to reduce the variance of $\hat{m}^{ens}_{*}$. Therefore, the result in Theorem <ref> is in line with common understanding of statistical modelling. In comparing with RF, this is equivalent to requiring each leaf to have more observations. However, as noticed by [48], for RF the variation due to the small leaf can be averaged out in the procedure of ensemble. Thus,  $ \nu $ can be smaller than the “optimal”.   In our calculation, $ \nu = 0.2 $ is used.

§ NUMERICAL PERFORMANCE IN REAL DATA

In this section, we use 36 real data sets with continuous responses  and 36 data sets with binary categorical response (0 and 1) to demonstrate the performance of ePPR for regression prediction and classification respectively. These data sets can be found on one of the following websites:
(A) <https://archive.ics.uci.edu/ml/datasets>;
(B) <https://www.kaggle.com>;
(C) <http://lib.stat.cmu.edu/datasets>;
(D) <https://data.world>;
(E) <https://www.openml.org>.

Regression: average RPE based on 100 random partitions of each data set into training and test sets

Data sets (and sources)	   	N	   	p	    1cLM	   	1cPPR	   	SVM	   	1cRF	   	GRF	   	RLT	   	XGB	   	 ePPR
Concrete compressive strength (A)	   	1030	   	8	   	0.270	   	0.142	   	0.167	   	0.120	   	0.169	   	0.204	   	0.083	   	0.079

Wine quality (B)	   	1143	   	11	   	0.649	   	0.675	   	0.612	   	0.559	   	0.626	   	0.580	   	0.644	   	0.607

Boston house price (C)	   	506	   	13	   	0.252	   	0.227	   	0.184	   	0.143	   	0.206	   	0.175	   	0.154	   	0.125

Wild blueberry yield (B)	   	777	   	13	   	5.424	   	18.64	   	0.043	   	0.074	   	0.121	   	0.120	   	0.052	   	0.028

Body fat (B)	   	252	   	14	   	0.028	   	0.061	   	0.126	   	0.101	   	0.058	   	0.132	   	0.052	   	0.038

Baseball player salary (C)	   	263	   	16	   	0.483	   	0.605	   	0.302	   	0.249	   	0.304	   	0.270	   	0.315	   	0.247

Paris housing price (B)	   	10000	   	16	   	0.001	   	0.000	   	0.026	   	0.019	   	0.000	   	0.054	   	0.000	   	0.000

Insurance churn participants (B)	   	33908	   	16	   	0.794	   	0.921	   	0.886	   	0.715	   	0.778	   	0.743	   	0.857	   	0.758

House sales in King County, USA (B)	   	21613	   	18	   	0.228	   	0.217	   	0.226	   	0.161	   	0.199	   	0.195	   	0.157	   	0.134

Parkinsons (A)	   	5875	   	19	   	0.823	   	0.759	   	0.671	   	0.352	   	0.387	   	0.524	   	0.125	   	0.224

Life Expentency (WHO) (B)	   	1649	   	20	   	0.136	   	0.086	   	0.095	   	0.049	   	0.074	   	0.064	   	0.057	   	0.056

Bar crawl (B)	   	7590	   	21	   	0.231	   	0.225	   	0.193	   	0.170	   	0.235	   	0.192	   	0.172	   	0.144

Cancer death rate (D)	   	2332	   	32	   	0.522	   	0.614	   	0.536	   	0.502	   	0.556	   	0.524	   	0.524	   	0.187

Wave energy converters-Adelaide (A)	   	71998	   	32	   	0.841	   	0.614	   	0.203	   	0.280	   	0.453	   	0.302	   	0.326	   	0.237

Wave energy converters-Perth (A)	   	72000	   	32	   	0.885	   	0.708	   	0.258	   	0.345	   	0.561	   	0.375	   	0.395	   	0.290

Sidney house price (B)	   	30000	   	37	   	0.638	   	1.000	   	0.631	   	0.555	   	0.618	   	0.559	   	0.649	   	0.564

Russian house price (B)	   	21329	   	50	   	0.647	   	1.238	   	0.678	   	0.660	   	0.671	   	0.653	   	0.789	   	0.653

Facebook comment volume (A)	   	18370	   	52	   	0.727	   	2.225	   	0.808	   	0.567	   	0.630	   	0.625	   	0.648	   	0.536

Gold price prediction (B)	   	1718	   	74	   	0.018	   	0.011	   	0.009	   	0.005	   	0.010	   	0.005	   	0.008	   	0.006

Baseball player statistics (B)	   	4535	   	74	   	0.001	   	0.000	   	0.059	   	0.019	   	0.027	   	0.033	   	0.002	   	0.001

CNNpred (A)	   	1441	   	76	   	0.001	   	0.000	   	0.056	   	0.001	   	0.003	   	0.002	   	0.002	   	0.000

Superconductivity (A)	   	21263	   	81	   	0.285	   	0.257	   	0.226	   	0.161	   	0.221	   	0.169	   	0.180	   	0.184

Warsaw flat rent price (B)	   	3472	   	83	   	0.333	   	0.461	   	0.488	   	0.308	   	0.370	   	0.343	   	0.331	   	0.277

Year prediction MSD (A)	   	50000	   	90	   	0.785	   	1.635	   	0.796	   	0.853	   	0.906	   	0.850	   	1.001	   	0.754

Buzz in social media (A)	   	28179	   	96	   	0.200	   	0.460	   	0.655	   	0.103	   	0.204	   	0.132	   	0.132	   	0.103

Communities and crime (A)	   	1994	   	101	   	0.348	   	0.746	   	0.372	   	0.347	   	0.374	   	0.344	   	0.397	   	0.326

Residential building-Sales (A)	   	372	   	103	   	0.027	   	0.043	   	0.124	   	0.046	   	0.081	   	0.069	   	0.019	   	0.009

Residential building-Cost (A)	   	372	   	103	   	0.029	   	0.232	   	0.067	   	0.064	   	0.101	   	0.081	   	0.050	   	0.026

Geographical original-Latitude (A)	   	1059	   	116	   	0.827	   	1.471	   	0.765	   	0.737	   	0.831	   	0.759	   	0.823	   	0.748

Geographical original-Longitude (A)	   	1059	   	116	   	0.739	   	1.200	   	0.654	   	0.673	   	0.791	   	0.708	   	0.733	   	0.659

Blog feedback (A)	   	13063	   	249	   	0.645	   	1.136	   	0.742	   	0.556	   	0.613	   	0.588	   	0.650	   	0.559

Credit score (B)	   	80000	   	259	   	0.208	   	0.497	   	0.179	   	0.128	   	0.178	   	0.126	   	0.131	   	0.111

CT slices (A)	   	53500	   	380	   	0.206	   	0.214	   	0.087	   	0.080	   	0.173	   	0.075	   	0.094	   	0.088

UJIndoor-Longitude (A)	   	19937	   	465	   	0.099	   	0.143	   	0.087	   	0.020	   	0.043	   	0.030	   	0.037	   	0.030

UJIndoor-Latitude (A)	   	19937	   	465	   	0.131	   	0.341	   	0.121	   	0.033	   	0.083	   	0.066	   	0.052	   	0.054

Regression with categorical data (B)	   	169486	   	1032	   	0.621	   	6.957	   	0.701	   	0.576	   	0.659	   	0.583	   	0.635	   	0.556
3c|Average of RPE across all data sets   0.530	   	1.244	   	0.357	   	0.287	   	0.342	   	0.313	   	0.313	   	0.261

3c|no. of bests$^*$      2    1.75    3    9.5    0.25    1.5    1.25     16.75

11l * the number of data sets for which the method performs the best among all the methods;  if $Q$ methods have the

11lsame  best performance for one data set, count the no. as $1/Q$ for each of the methods.

Classification: average MR(%) based on 100 random partitions of each data set into training and test sets

Data sets (and sources)	   	N	   	p	   	GLM	   	PPR	   	SVM	   	RF	   	GRF	   	RLT	   	XGB	   	Ada	   	ANN	   	ePPR
Indian liver patient (A)	   	579	   	10	   	29.52	   	31.56	   	29.16	   	29.67	   	28.82	   	29.25	   	32.02	   	30.27	   	29.22	   	29.31

MAGIC Gamma telescope (A)	   	19020	   	10	   	20.58	   	15.85	   	15.62	   	14.68	   	16.42	   	14.73	   	17.90	   	14.89	   	36.97	   	13.73

Heart disease (A)	   	270	   	13	   	16.41	   	21.53	   	17.29	   	17.47	   	17.99	   	17.06	   	22.81	   	17.92	   	44.11	   	18.20

Australian credit approval (A)	   	690	   	14	   	13.75	   	14.76	   	14.57	   	13.12	   	14.23	   	13.50	   	15.68	   	13.27	   	14.96	   	13.44

EEG eye state (A)	   	14980	   	14	   	37.40	   	24.18	   	25.39	   	18.04	   	24.74	   	19.30	   	22.57	   	20.87	   	44.84	   	11.05

seismic-bumps (A)	   	2584	   	15	   	6.62	   	7.47	   	6.59	   	6.78	   	6.56	   	6.64	   	9.08	   	6.85	   	6.50	   	7.34

Autism screening adult (A)	   	702	   	16	   	0.00	   	0.04	   	2.36	   	0.00	   	0.00	   	0.28	   	0.00	   	0.00	   	15.90	   	0.01

Bank marketing (A)	   	45211	   	16	   	11.18	   	11.92	   	11.09	   	10.61	   	11.50	   	10.92	   	12.27	   	10.91	   	11.94	   	10.91

Diabetic retinopathy debrecen (A)	   	1151	   	19	   	26.05	   	26.80	   	30.78	   	31.74	   	33.51	   	32.38	   	34.87	   	31.41	   	48.39	   	25.12

Performance prediction (B)	   	1340	   	19	   	29.33	   	31.43	   	29.50	   	31.27	   	29.94	   	30.45	   	36.26	   	31.14	   	40.07	   	29.78

Climate model simulation crash (E)	   	540	   	20	   	8.99	   	12.91	   	8.58	   	8.71	   	8.58	   	8.59	   	11.21	   	9.12	   	8.72	   	9.27

Default of credit card clients (A)	   	30000	   	24	   	19.68	   	20.43	   	19.15	   	19.15	   	18.30	   	18.97	   	24.13	   	19.22	   	22.12	   	19.24

Parkinson multiple sound (B)	   	1208	   	26	   	33.06	   	33.85	   	29.60	   	27.61	   	30.39	   	28.37	   	30.83	   	28.49	   	43.25	   	26.02

Pistachio (B)	   	2148	   	28	   	7.50	   	8.02	   	8.19	   	10.80	   	12.63	   	10.91	   	11.69	   	9.34	   	44.36	   	6.57

Breast cancer (B)	   	569	   	30	   	2.95	   	5.46	   	2.67	   	4.24	   	5.97	   	4.02	   	5.52	   	3.85	   	37.95	   	2.41

Phishing websites (A)	   	11055	   	30	   	7.80	   	8.73	   	6.86	   	5.67	   	8.01	   	6.61	   	6.81	   	6.25	   	7.74	   	6.29

Ionosphere (A)	   	351	   	33	   	12.65	   	15.75	   	6.14	   	7.21	   	7.65	   	7.25	   	10.79	   	7.52	   	12.65	   	5.61

QSAR biodegradation (A)	   	1055	   	41	   	14.13	   	16.16	   	13.02	   	13.38	   	16.08	   	13.18	   	15.85	   	13.88	   	30.23	   	12.55

Spambase (A)	   	4601	   	57	   	7.04	   	9.19	   	8.47	   	6.16	   	8.30	   	6.29	   	8.35	   	6.32	   	41.62	   	5.87

Mice protein expression (A)	   	1047	   	70	   	2.94	   	3.89	   	1.00	   	1.35	   	7.47	   	1.28	   	6.89	   	3.10	   	6.70	   	0.54

Ozone level detection (A)	   	1847	   	72	   	6.66	   	11.72	   	7.00	   	6.44	   	7.02	   	6.70	   	7.45	   	6.37	   	7.06	   	6.35

Insurance company benchmark (A)	   	5822	   	85	   	5.99	   	15.22	   	5.99	   	6.54	   	5.99	   	6.20	   	9.03	   	6.33	   	6.16	   	6.77

Company bankruptcy (B)	   	6819	   	94	   	3.21	   	5.34	   	3.22	   	3.13	   	3.22	   	3.13	   	3.81	   	3.15	   	3.25	   	3.39

Hill-valley (A)	   	1212	   	100	   	50.48	   	4.32	   	50.00	   	40.76	   	46.00	   	41.93	   	38.13	   	43.15	   	51.93	   	1.58

Musk (A)	   	6598	   	166	   	7.22	   	10.24	   	5.68	   	5.63	   	9.63	   	6.21	   	7.06	   	6.15	   	15.45	   	4.34

ECG heartbeat categorization (B)	   	14550	   	186	   	19.18	   	25.45	   	18.06	   	8.88	   	14.86	   	8.88	   	13.36	   	10.33	   	13.26	   	7.30

Arrhythmia (A)	   	420	   	192	   	24.79	   	37.20	   	23.64	   	19.67	   	23.05	   	21.86	   	24.89	   	18.76	   	45.29	   	20.19

Financial indicators of US stocks (B)	   	986	   	216	   	0.62	   	12.93	   	13.94	   	0.94	   	1.33	   	4.74	   	0.05	   	0.05	   	23.53	   	0.31

Madelon (A)	   	2000	   	500	   	38.91	   	48.10	   	42.46	   	34.52	   	37.54	   	32.98	   	32.78	   	32.46	   	49.71	   	33.44

Human activity recognition (A)	   	2633	   	561	   	0.00	   	0.00	   	0.10	   	0.08	   	0.25	   	0.16	   	0.01	   	0.01	   	0.01	   	0.00

Gina agnostic (E)	   	3468	   	970	   	14.38	   	38.87	   	11.47	   	9.25	   	13.36	   	9.20	   	13.80	   	9.18	   	45.33	   	9.93

QSAR androgen receptor (A)	   	1687	   	1024	   	10.36	   	35.13	   	10.09	   	9.71	   	11.38	   	9.55	   	11.69	   	9.71	   	11.38	   	9.59

QSAR oral toxicity (A)	   	8992	   	1024	   	8.10	   	36.98	   	7.71	   	7.26	   	8.18	   	7.31	   	9.08	   	7.42	   	8.33	   	7.33

Predicting a biological response (B)	   	3751	   	1776	   	25.20	   	34.10	   	25.98	   	23.15	   	25.17	   	22.45	   	27.55	   	23.60	   	33.06	   	23.62

GISETTE (A)	   	6000	   	4955	   	4.84	   	16.36	   	5.29	   	5.16	   	8.34	   	4.77	   	7.72	   	4.69	   	39.35	   	4.30

Arcene (A)	   	200	   	9961	   	28.69	   	44.27	   	25.93	   	19.76	   	27.93	   	22.03	   	28.10	   	21.84	   	38.06	   	13.69
3c|average of MR(%) across all data sets   15.45	   	19.34	   	15.07	   	13.29	   	15.29	   	13.56	   	15.83	   	13.55	   	26.09	   	11.26

3c|no. of bests   2.87    0.33    0.83    4.7    3.03    2.5    0.7    3.7    1    16.33

Average RPE (or MR %) for each method across the 36 variable-augmented regression data sets (or 36 variable-augmented classification data sets) with p = 10000 and the number of data sets for which the method has the best performance (no. of bests) among all the methods

LM   SVM   RF   GRF   RLT   XGB   ePPR
the 36  regression       average RPE   0.469   0.860   0.394   0.748   0.446   0.407   0.348
data sets      no. of bests    1    0    6    0    1    2    26

the 36 classification         average MR(%)   17.95   28.12   22.97   25.45   17.23   20.17   15.90
2-9 data sets     no. of bests    5.06   0.76   0.76   0.76   6.56   4.5   17.56

For regression, we compare ePPR with the original PPR, the support vector machine (SVM), the random forest (RF), the generalized random forest (GRF), the reinforcement learning trees (RLT) and the extreme gradient boosting (XGB). For classification, we compare ePPR with PPR, SVM, RF, GRF, RLT, XGB, the adaptive boosting (Ada) and artificial neural network (ANN). To apply PPR, GRF and ePPR to classification, as they were developed for regression, we simply treat classes $0$ and $1$ as values and apply the methods directly to predict the response,  $ \hat y $, then classify the data according to $1(\hat y>0.5)$, where $ 1(.) $ is the indicator function. Extension to multiclass classification can be easily made.

The following functions and packages in R are used for the calculation of the methods:  [24] is used for Lasso, including LM for linear regression and GLM (logistic model) for classification,  function in package  [22] for SVM,   [43] for RF,   in package  [2] for GRF,   [55] for RLT,  package  [16] for XGB,  [19] for ada, and  function in package  [5] for ANN. For all the R functions  and  packages, their default values of tuning parameters are used. For ANN, however, there is no default size, i.e. numbers of hidden layers or nodes of each layer. We use 5-fold CV to select the size from candidate sizes of equal number of neurons $w$ in all hidden layers $L$ with $ L = 1, ..., 10 $ and $  w= 5, 6, ..., 50 $. For ePPR, $B=500$ is used.

For each data set, we randomly partition it into training set and test set. The training set consists of $n=\min(\lfloor 2N/3\rfloor,1000)$ randomly selected observations, where $N$ is the number of observations in the original data sets, and the remaining observations form the test set. For regression, the relative prediction error, defined as

$$RPE=\sum_{i\in \text{test set}}(\hat{y}_i-y_i)^2/\sum_{i\in \text{test set}}(\bar{y}_{\text{train}}-y_i)^2,$$

where $\bar{y}_{\text{train}}$ is naive predictions based on the average of $y$ in the training sets, is used to measure the performance of a method. For classification, the misclassification rate, defined as

$$MR=\sum_{i\in \text{test set}} 1(1(\hat{y}_i>0.5) \neq y_i) /(N-n),$$

is used to measure the performance. For each data set, the random partition is repeated 100 times,  and averages of the RPEs or MRs are calculated to compare different methods. The calculation results are listed in Table <ref> and Table <ref>. The smallest RPE or MR for each data set is highlighted in bold font.

By comparing the prediction errors, RPEs or MRs, of all the methods, we have the following observations. The classical methods SVM, RF and XGB have reasonably good performance on nearly all the data sets in comparison with LM and PPR. In particular, RF has been commonly recognized as the most robust and efficient learning method, especially for tabular data; our extensive calculation also supports the recognition. The recent variants of RF, such as GRF and RLT, also show their effectiveness in the prediction.
Although PPR has good performance on some data sets, it is unstable and has much bigger prediction errors than the other methods for most the data sets.  This could be the reason why PPR is not popularly received in practice.  The proposed ePPR is quite stable and attains the smallest RPE and MR in most data sets as indicated in Table <ref> and Table <ref>. The competency of ePPR is also verified by the fact that it has the smallest average of RPEs (or MSs) across all the 36 data sets amongst all the methods. The numbers of data sets for which a method performs the best among the all the competitors, denoted by no. of bests, also suggest that ePPR outperforms the other methods for most of the data sets. The reason for better performance of ePPR over RF might be the fact that ePPR is smooth, while RF is a step function [52].

It is also noticed that RF has good performance for high dimensional data [14, 13]. To see the performance of ePPR in high dimension, we add randomly generated i.i.d. normal variables, called augmented variables, to the original data  so that all together there are $p=10000$ predictors,  and then compare all the methods. Since PPR, Ada and ANN require huge memory with big $ p $, they are not included for comparison in Table <ref>. For the other methods, the overall performance across all the data sets are listed in Table <ref>. We can see that ePPR still has the best performance in most of the data sets for both regression and classification. Although all methods suffer from the high dimension, ePPR is less influenced: compared to RF, ePPR has about 9% improvement in prediction and 15% improvement in classification in Table <ref> and Table <ref> respectively, while ePPR is about 12% or 30% respectively better than RF in Table <ref>, indicating that the relative improvement increases  with the dimension of the data.  The number of data sets for which ePPR outperforms the competitors,  no. of bests, also increases with the dimension of the data.

§ PROOFS

The following notations are frequently used in the proofs. Let $c$ be a constant which may change from line to line. Function $c(a)$  only depends on $a $  and $ c(a, b) $ only depends on $a $ and $b$.
We use $\|f\|_{L^q(A)}$ to denote the $L^q$ norm of $f\in L^q(A), q> 1$, where $A$ is a subset of Euclidean space.

Define the population norm $\|f\|=\left(\E(f^2(X))\right)^\frac{1}{2}$ and  the inner product $\langle f,g\rangle=\E(f(X)\cdot g(X))$ if $\E\left( f^2(X)\right)<\infty$ and $\E\left( g^2(X)\right)<\infty$, and define the empirical norm $\|f\|_n=\left(\frac{1}{n}\sum_{i=1}^n{f^2(X_i)}\right)^\frac{1}{2}$ and empirical inner product $\langle f,g\rangle_n=\frac{1}{n}\sum_{i=1}^n{\left(f(X_i)\cdot g(X_i)\right)}$ for any real functions $f$ and $g$. Let $\|Y-f\|_n=\left(\frac{1}{n}\sum_{i=1}^n{(Y_i-f(X_i))^2}\right)^\frac{1}{2}$ for any real function $f$. Define $\|h\|_{n,EL^1}=\sum_{\iota=1}^{\infty} |\beta_{\iota}| \cdot \|S_\iota(\theta_{\iota}^\top  x_{\C_\iota})\|_n$ for any function $h=\sum_{\iota=1}^{\infty} \beta_{\iota} \cdot S_\iota(\theta_{\iota}^\top  x_{\C_\iota})\in \mathcal{H}_{\infty,q_0,n}$ with

\begin{equation*} %\label{sampleHspace}
\mathcal{H}_{\infty, q_0,n}=\Big\{\sum_{i=1}^{\infty}{S_i(\theta_i^\top  x_{\C_i})}:S_i\in {\cal S}_{u,s}[-1,1], \theta_i\in \Theta^{q_0}, \C_i\subset\{1,\ldots,p\}, card(\C_i)=q_0\Big\}.
\end{equation*}

Define $ \mathcal{H}_{k_n,p,n},\mathcal{H}_{N_n,p,n}$ in the same way as $\mathcal{H}_{\infty, q_0,n}$ except that their subscripts are changed respectively. Meanwhile, $\|h\|_{n,EL^1}$ for $h\in \mathcal{H}_{k_n,p,n}\cup\mathcal{H}_{N_n,p,n}$ is also defined in the same manner as  $h\in \mathcal{H}_{\infty, q_0,n}$.
For any estimator  $\hat m(x) $ of $ m(x) $, denote the integrated squared error (ISE) by

$$  {{\cal I}}(\hat m(x)) =  \int|\hat m(x) -m(x)|^2d\mu(x),  $$

and mean integrated squared error (MISE) by  $ \MSE(\hat m(x)) = \E_{\mathcal{D}_{n}} {{\cal I}}(\hat{m}(x)) $. For ensemble estimator, we also need to take average over the random subsets space $ \Xi_{k_n} $, the corresponding MISE is also denoted by

\MSE (\hat{m}(x)) = \E_{\mathcal{D}_{n},\Xi_{k_n}} {{\cal I}}(\hat m(x)).

Let $x_1,\ldots, x_n \in\mathbb{R}^p$ and set $x_1^n=\{x_1,\ldots, x_n\}$. Let $\mathcal{G}$ be a class of functions $g:\mathbb{R}^p\to\mathbb{R}$. An $L_\ell$-$\epsilon$-cover of $\mathcal{G}$ on $x_1^n$ is defined in [30]. Denote by $\mathcal{N}_\ell(\epsilon,\mathcal{G},x_1^n)$ the covering number. In this paper, $ \ell = 1 $ is used. We further simplify the covering number as  $\mathcal{N}(\epsilon,\mathcal{G},x_1^n)$. See [30] or Definition <ref> for more details in the supplementary material, abbreviated as  hereafter.

§.§ Proof of Theorem <ref>
Recall $ V_n = span\{ v_j, j=1, ..., J_n\} $.
Define a function space

\begin{equation*}\label{definitionMn}
\mathcal{M}_n :=\left\{ \sum_{\iota=1}^{K_n}{g_\iota(x^\top \theta_\iota)}: g_\iota\in V_n,\ \theta_\iota\in \Theta^{p}\right\}
\end{equation*}

and its truncated version ${\mathcal{M}^{t_n}_n}=\{\max\{\min\{g,t_n\},-t_n\}: g\in\mathcal{M}_n\}$, and  a subset of $ \mathcal{M}_n$,
\mathcal{F}_n:=\text{span}\{ g(x^\top \theta): g\in V_n,\ \theta \in\Theta_n\}.
It is easy to see that  $\mathcal{F}_n$ is also a closed sub-manifold of $L^2(\B^p)$ equipped with $L^2$ norm.

The proof of this theorem is established on Theorem 1 of [3] or Lemma <ref> in  which proved that

\begin{equation}\label{error}
\begin{aligned}
\MSE(\hat{m}_n(x)) \leq & c\cdot\inf_{v\in \mathcal{M}_n}\int|m(x)-v(x)|^2dx \\
& +\frac{c\ln^2 n}{n}\cdot \sup_{x_1^n}\ln\left(\mathcal{N}(1/(80n t_n),\mathcal{M}^{t_n}_n,x_1^n)\right),
\end{aligned}
\end{equation}

where  $\mathcal{N}$ is the covering number defined at the beginning of this section. The first part of the RHS of (<ref>) is about the approximation of $m$ by $\mathcal{M}_n$; and the second part of the RHS of (<ref>) is due to the size of approximation function space $\mathcal{M}_n$. Next, we consider these two parts separately.

For the first part,  according to [45], there is a function $R_n\in \mathcal{F}_n\subset \mathcal{M}_n$ such that

\begin{equation}\label{Tec1}
||m-R_n||_{L^2(\B^p)} \leq c(\lambda)J_n^{-\lambda}||m||_{H^\lambda(\B^p)}, \ \forall m\in H^{\lambda}(\B^p),
\end{equation}

where $\lambda= s +(p-1)/2$ if inequality (<ref>) is satisfied.  Next we extend (<ref>) to any function $m\in H^r(\B^p), r\leq \lambda$. For any $m\in L^2(\B^p)$, let $\mathcal{Q}_n(m)$ be the projector of $m$ onto the linear space $\mathcal{F}_n$ in terms of $L^2$ norm. It is easy to see that projection operator $\mathcal{Q}_n$ is linear from $L^2(\B^p)$ to $L^2(\B^p)$ and

\begin{equation}\label{Tec2}
\|m-\mathcal{Q}_n(m)\|_{L^2(\B^p)}\leq {\|m-v\|_{L^2(\B^p)}}
\end{equation}

for any $v\in \mathcal{F}_n$. Combination of (<ref>) and (<ref>) leads to

\begin{equation}\label{Tec3}
\|m-\mathcal{Q}_n(m)\|_{L^2(\B^p)}\leq c(\lambda)J_n^{-\lambda}||m||_{H^\lambda(\B^p)}
\end{equation}

provided that $m\in H^{\lambda}(\B^p)\subset L^2(\B^p)$. Now, define operator $\mathcal{L}:m\to m-\mathcal{Q}_n(m)$. Equation (<ref>) indicates that $\mathcal{L}$ is a linear map from $ H^{\lambda}(\B^p)$ to $L^2(\B^p)$ with norm smaller than $c(\lambda)J_n^{-\lambda}$, while  (<ref>) indicates that $\mathcal{L}$ is also a linear map from $L^2(\B^p)$ to $L^2(\B^p)$ with norm not larger than $1$, which equals 1 when $v=0$. Therefore, by the property of Banach-space interpolation  (see [11]), we know that $\mathcal{L}$ is a bounded linear operator from $H^{r}(\B^p)=[H^{\lambda}(\B^p), H^{0}(\B^p)]_{1-\frac{r}{\lambda},2}$ to
$L^2(\B^p)=[L^2(\B^p),L^2(\B^p)]_{1-\frac{r}{\lambda},2}$ with the norm smaller than $c(\lambda)\cdot J_n^{-r}$, where $[\cdot, \cdot]_{1-\frac{r}{\lambda},2}$ denotes the interpolation space (see [1] or [11] for more details). In other words, for any $m\in H^r(\B^p)$ we have

\|m-\mathcal{Q}_n(m)\|_{L^2(\B^p)}\leq c(\lambda)J_n^{-r}\|m\|_{H^r(\B^p)},

which implies

\begin{equation}\label{ConclusionPart1}
\inf_{v\in \mathcal{F}_n}\|m-v\|_{L^2(\B^p)}\leq c(\lambda)J_n^{-r}||m||_{H^r(\B^p)}.
\end{equation}

Now, we consider  the second part of the RHS of (<ref>). Let $VC(\cdot)$ denote the VC dimension of a function class; see [49] for the definition. With the standard calculations, as shown in Lemma <ref> in ,  we have

\begin{equation}\label{myVC}
VC({\mathcal{M}^{t_n}_n}) \leq VC({\mathcal{M}_n}) \leq c\cdot pK_nJ_n\ln(pK_nJ_n),
\end{equation}

where $c>0$. Taking $K_n\asymp J_n^{p-1}$, by Lemma 9.2 in [30] we have

\begin{align}\label{ConclusionPart2}
\mathcal{N}(1/(80n t_n),\mathcal{M}^{t_n}_n,x_1^n) &\leq 3\left( \frac{4e t_n}{1/(80n t_n)}\ln\left(\frac{6e t_n}{1/(80n t_n)} \right) \right)^{VC({\mathcal{M}^{t_n}_n})} \\
& \leq 3\left( \frac{4e t_n}{1/(80n t_n)}\ln\left(\frac{6e t_n}{1/(80n t_n)} \right) \right)^{VC({\mathcal{M}_n})}\nonumber \\
& \leq 3\left( \frac{4e t_n}{1/(80n t_n)}\ln\left(\frac{6e t_n}{1/(80n t_n)} \right) \right)^{c\cdot pK_nJ_n\ln(pK_nJ_n)}\nonumber \\
& \leq 3\left( 480e n t_n^2 \right)^{c\cdot J_n^p \ln(J_n^p)}.\nonumber
\end{align}

Substituting  (<ref>) and  (<ref>) into (<ref>) yields

\begin{aligned}
\MSE(\hat{m}_n(x)) &\leq c(\lambda)\cdot J_n^{-2r}+c\cdot\frac{\ln^2n \cdot J_n^p\cdot\ln(J_n)\cdot\ln(nt_n^2)}{n}\\
&\leq c(\lambda)\cdot J_n^{-2r}+c\cdot \frac{J_n^p\cdot\ln(J_n)\cdot\ln^3n}{n}.
\end{aligned}

If we set $J_n\asymp \left(\frac{n}{\ln^4 n}\right)^{\frac{1}{2r+p}}$, the RHS of above inequality is bounded by $c(\lambda,p,r)\cdot\left(\frac{\ln^4 n}{n}\right)^{\frac{2r}{2r+p}}$, where the constant $c(\lambda,p,r)$  only depends on $\lambda,p,r$. This completes the proof.   \(\Box\)

§.§ Proof of Corollary <ref>

The verification of (<ref>) for $S_{u,s}[-1,1]$ is a direct conclusion of Theorem 3 of [47]. Then, Corollary <ref> follows from Theorem <ref> and the fact that each $B_{s,j}$ consists of $s+2$ pieces of polynomials with degree not larger than $s$.  \(\Box\)

§.§ Proof of Theorem <ref>

Proofs for rPPR.oga, rPPR.aga and rPPR.rga are considered separately below.

§.§.§ rPPR.oga

Define a collection of functions, also called dictionary ([4]),

\mathcal{D}ic_{\A}=\{g(\theta^\top  X_{\A}): \theta\in \Theta^q, \|g(\theta^\top  X_{\A})\|=1\}

for any subset $\A\subset\{1,\ldots,p\}$ with $card(\A)=q$. Denote by $\delta_k=m-m_k$ the approximation error in the $k$-th step  for each $k\geq 0$. Let

g_{k}^o:=\argmax_{v\in \mathcal{D}ic_{\A_k}}{|\langle \delta_{k-1},v\rangle|}

be the best projection  at the $k$-th step, where $k\ge 1$, satisfying $g_k(\theta_{k}^\top X_{\A_k})=g_{k}^o\cdot \|g_k(\theta_{k}^\top X_{\A_k})\|$ and

g_k^o=\argmin_{v\in \mathcal{D}ic_{\A_k}}{\|\delta_{k-1}-\langle \delta_{k-1},v\rangle\cdot v\|}.

First, it follows from the property of projection that

\|\delta_k\|^2\leq \|\delta_{k-1}-\langle \delta_{k-1},g_{k}^o\rangle\cdot g_{k}^o\|^2=\|\delta_{k-1}\|^2-|\langle \delta_{k-1},g_{k}^o\rangle|^2.

Because conditional on $\Xi_{k-1}$, $\delta_{k-1}$ is independent with $\Xi_{k}$,   taking conditional expectation over $\Xi_{k}$ on both sides of the above inequality leads to

\begin{equation}\label{Math111}
\E_{\Xi_k|\Xi_{k-1}}\|\delta_k\|^2\leq \|\delta_{k-1}\|^2- \E_{\Xi_k|\Xi_{k-1}}|\langle \delta_{k-1},g_{k}^o\rangle|^2.
\end{equation}

Let ${\C}_0= \argmax_{{\C}\subset\{1,\ldots,p\},Card({\C})=q}\sup_{g\in \mathcal{D}ic_{\C}}|\langle \delta_{k-1},g\rangle|$.  Because the number of subsets $\mathcal{C}$ which contains ${\C}_0$ is only 1, we have

\begin{equation}\label{Math131}
\begin{aligned}
\E_{\Xi_k|\Xi_{k-1}} |\langle \delta_{k-1},g_{k}^o\rangle|^2 = \sum_{\mathcal{C}\subset\{1,\ldots, p\}}{\frac{1}{\binom{p}{q}}|\langle \delta_{k-1},g_{k,\mathcal{C}}^o\rangle|^2} \geq \frac{1}{\binom{p}{q}}\sup_{v\in \mathcal{D}ic_{{\C}_0}}{|\langle \delta_{k-1},v\rangle|^2},
\end{aligned}
\end{equation}

where $g_{k,\mathcal{C}}^o:=g_{k,\mathcal{C}}(\theta_{k,\C}^\top X_{\C})=\argmax_{v\in \mathcal{D}ic_{\mathcal{C}}}{|\langle \delta_{k-1},v \rangle|}$ for any fixed $\C\subset\{1,\ldots, p\}$ with $ card(\mathcal{C})=q$. In conclusion, (<ref>) and (<ref>) imply that

\begin{equation}\label{Math23}
\|\delta_{k-1}\|^2-\E_{\Xi_k|\Xi_{k-1}}\|\delta_k\|^2\geq \frac{\sup_{g\in \mathcal{D}ic_{{\C}_0}}{|\langle \delta_{k-1},g\rangle|^2}}{c(p,q)}.
\end{equation}

Second, we also have

\begin{align}
\|\delta_{k-1}\|^2&=\langle \delta_{k-1},h+m-h\rangle=\langle \delta_{k-1},h\rangle+\langle \delta_{k-1},m-h\rangle\label{Math21}\\
&\leq \|h\|_{EL^1}\cdot\sup_{{\C}\subset\{1,\ldots,p\},Card({\C})=q}\sup_{v\in \mathcal{D}ic_{\C}}|\langle \delta_{k-1},v\rangle|+\|\delta_{k-1}\|\cdot\|m-h\|,\nonumber
\end{align}

where in the first equation, the fact that $\delta_{k-1}$ is orthogonal to $m_{k-1}$ is used; and in the last inequality, Cauchy-Schwarz inequality is used. Using the fact that $\|\delta_{k-1}\|\cdot\|m-h\|\leq \frac{1}{2}(\|\delta_{k-1}\|^2+\|m-h\|^2)$, equation (<ref>) implies that

\begin{equation}\label{Math22}
\sup_{g\in \mathcal{D}ic_{{\C}_0}}|\langle \delta_{k-1},g\rangle|\geq \frac{\|\delta_{k-1}\|^2-\|m-h\|^2}{2\cdot\|h\|_{EL^1}}.
\end{equation}

It follows from (<ref>) and (<ref>) that

\begin{equation}\label{Math24}
\|\delta_{k-1}\|^2-\E_{\Xi_k|\Xi_{k-1}}\|\delta_k\|^2\geq \frac{\max{(\|\delta_{k-1}\|^2-\|m-h\|^2,0})^2}{4\cdot\left(\|h\|_{EL^1}\cdot{c(p,q)}\right)^2}.
\end{equation}

Taking expectation w.r.t. $\Xi_{k-1}$ on both sides of (<ref>), by the law of iterated expectation and Jensen's inequality we have

\begin{equation}\label{Math25}
\E_{\Xi_{k}}\|\delta_k\|^2\leq  \E_{\Xi_{k-1}}\|\delta_{k-1}\|^2-\frac{ \left(\E_{\Xi_{k-1}}\max{(\|\delta_{k-1}\|^2-\|m-h\|^2,0})\right)^2}{4\cdot\left(\|h\|_{EL^1}\cdot{c(p,q)}\right)^2}.
\end{equation}

Let $\Delta_\iota=\E_{\Xi_{\iota}} \left( \|\delta_{\iota}\|^2-\|m-h\|^2 \right),\ 1\le \iota\le k$. The theorem  is equivalent to

\begin{equation} \label{equiv}
\Delta_k\leq \frac{M}{k+1}
\end{equation}

where $M=4\left(\|h\|_{EL^1}\cdot{c(p,q)}\right)^2$. To prove (<ref>), consider  two scenarios: $\Delta_{k-1}\leq 0$ and $\Delta_{k-1}\ge 0$. For the first scenario, i.e. $\Delta_{k-1}\leq 0$, due to the projection we have $\E_{\Xi_{k}|\Xi_{k-1}}\|\delta_{k}\|^2\leq \|\delta_{k-1}\|^2$ a.s.  It follows that

\begin{equation}\label{Theorem2ogalast}
\E_{\Xi_{k}}(\|\delta_{k}\|^2)=\E_{\Xi_{k-1}}\E_{\Xi_{k}|\Xi_{k-1}}\|\delta_{k}\|^2\leq \E_{\Xi_{k-1}}\|\delta_{k-1}\|^2.
\end{equation}

$$ \Delta_k =  \E_{\Xi_{k}} \left( \|\delta_{k}\|^2-\|m-h\|^2 \right)  \le \E_{\Xi_{k-1}} \left( \|\delta_{k-1}\|^2-\|m-h\|^2 \right) \leq 0.$$

Thus, (<ref>) holds. Now consider the second scenario, $\Delta_{k-1}\geq 0$. Since $\{\Delta_{i}\}_{i=0}^{k-1}$ is a decreasing sequence by a similar argument to (<ref>), we have $\Delta_{\iota}\geq 0$  for each $0\leq \iota\leq k-1$. Therefore, it follows from (<ref>)  that

\begin{equation}\label{Math26}
\Delta_\iota\leq  \Delta_{\iota-1}\left(1-\frac{\Delta_{\iota-1}}{M}\right),\ 1 \leq \iota\leq k.
\end{equation}

If $\Delta_0\leq M$, we know $\Delta_k\leq \frac{M}{k+1}$  by mathematical induction, which is (<ref>); if $\Delta_0 > M$, then  $\Delta_1<0$ by (<ref>), which leads to $\Delta_k<0$ for all $k\geq 1$ and (<ref>) also holds. In conclusion, the result for rPPR.oga is true.

§.§.§ rPPR.aga

The proof for rPPR.aga is similar to that of rPPR.oga.

§.§.§ rPPR.rga
Use the notations above and define the residual by $\tilde{\delta}_k=m-m_k$ for each $k\ge 0$. For any given $h(x)=\sum_{\iota=1}^\infty \beta_{\iota}\cdot v_\iota (\theta_{\iota}^\top x_{\mathcal{C}_\iota })\in \mathcal{H}_{q}$ with $\beta_\iota\ge 0$ and $\|v_\iota (\theta_{\iota}^\top x_{\mathcal{C}_\iota })\|=1$,  construct a  discrete random variable for each index $\iota\ge 1$:

\xi_{k,\iota}=\left\{
\begin{array}{ll}
\frac{2\cdot c(p,q)}{k+2}\cdot \|h\|_{EL^1}, & \mathcal{C}_\iota \subset \A_k,\\
0, & \text{otherwise},\\
\end{array}
\right

where $\A_k$ is a random subset of $ \{1, ..., p\} $ defined in Section <ref> and $\C_\iota $ is a fixed subset. Now, we fix the index $\iota$. Let $g_{\A_k}$ be a random function in the dictionary $\mathcal{D}ic_{\A_k}$:   $g_{\A_k}= g_\iota (\theta_{\iota}^\top x_{\mathcal{C}_\iota })\in \mathcal{D}ic_{\A_k} $ if ${\mathcal{C}_\iota }\subset \A_k$, and $g_{\A_k}=0$ if ${\mathcal{C}_\iota }\not\subset \A_k$. Based on the Algorithm rPPR.rga,
we have

\begin{align}\label{RGAproof2a}
\|\tilde{\delta}_k\|^2&\leq \|m-\alpha_k\cdot m_{k-1}-\xi_{k,\iota}\cdot g_{\A_k}\|^2\\
&=\|\alpha_k\cdot \tilde{\delta}_{k-1}+\bar{\alpha}_k\cdot m-\xi_{k,\iota}\cdot g_{\A_k}\|^2\nonumber\\
&=\alpha_k^2\cdot\|\tilde{\delta}_{k-1}\|^2+2\alpha_k\cdot \langle \tilde{\delta}_{k-1},\bar{\alpha}_k \cdot m-\xi_{k,\iota}\cdot g_{\A_k}\rangle+\|\bar{\alpha}_k\cdot m-\xi_{k,\iota}\cdot g_{\A_k}\|^2\nonumber\\
&=\alpha_k^2\cdot\|\tilde{\delta}_{k-1}\|^2+2\alpha_k\cdot \langle\tilde{\delta}_{k-1},\bar{\alpha}_k m-\xi_{k,\iota}\cdot g_{\A_k}\rangle +\|\bar{\alpha}_k\cdot (m-h)+\bar{\alpha}_k \cdot h-\xi_{k,\iota}\cdot g_{\A_k}\|^2\nonumber\\
&=\alpha_k^2\cdot\|\tilde{\delta}_{k-1}\|^2+2\alpha_k\cdot \langle \tilde{\delta}_{k-1},\bar{\alpha}_k\cdot m-\xi_{k,\iota}\cdot g_{\A_k}\rangle+\|\bar{\alpha}_k\cdot (m-h)\|^2\nonumber\\
&\ \ \ \ +2\bar{\alpha}_k\cdot \langle m-h, \bar{\alpha}_k\cdot h-\xi_{k,\iota}\cdot g_{\A_k}\rangle + \bar{\alpha}_k^2\cdot \|h\|^2 +\xi_{k,\iota}^2-2\bar{\alpha}_k \cdot \langle h,\xi_{k,\iota}\cdot g_{\A_k}\rangle\nonumber
\end{align}

where $\bar{\alpha}_k=1-\alpha_k$ and $k\geq 1$. Noe that $\E_{\Xi_k|\Xi_{k-1}}(\xi_{k,\iota}\cdot g_{\A_k})=\bar{\alpha}_k\cdot \|h\|_{EL^1}\cdot v_\iota (\theta_{\iota}^\top  x_{\C_\iota })$.  Taking conditional expectation w.r.t. $\Xi_{k}|\Xi_{k-1}$ on both sides of equation (<ref>) we have

\begin{multline}\label{RGAproof2b}
\E_{\Xi_k|\Xi_{k-1}}\|\tilde{\delta}_k\|^2\leq \alpha_k^2\cdot \|\tilde{\delta}_{k-1}\|^2+2\alpha_k\cdot \langle \tilde{\delta}_{k-1},\bar{\alpha}_k\cdot m-\bar{\alpha}_k\cdot \|h\|_{EL^1}\cdot v_\iota (\theta_{\iota}^\top x_{\C_\iota })\rangle\\+\bar{\alpha}_k^2\cdot \|m-h\|^2+
2\bar{\alpha}_k^2\cdot \langle m-h, h-\|h\|_{EL^1}\cdot v_\iota (\theta_{\iota}^\top x_{\C_\iota })\rangle\\+\bar{\alpha}_k^2\cdot \|h\|^2+\bar{\alpha}_k^2\cdot c^2(p,q)\cdot \|h\|_{EL^1}^2-2\bar{\alpha}_k^2\cdot \langle h, \|h\|_{EL^1}\cdot v_\iota (\theta_{\iota}^\top  x_{\C_\iota })\rangle,
\end{multline}

where $\|h\|_{EL^1}=\sum_{\iota=1}^\infty{\beta_\iota}$. Since  (<ref>) holds for all $\iota\geq 1$, it also holds on their average with the  weight $\frac{\beta_\iota}{\|h\|_{EL^1}}$ for the $\iota $-th part. Therefore,

\begin{align}\label{RGAproof2c}
\E_{\Xi_k|\Xi_{k-1}}\|\tilde{\delta}_k\|^2
&\leq \alpha_k^2\cdot \|\tilde{\delta}_{k-1}\|^2+ 2\alpha_k\bar{\alpha}_k\cdot \langle \tilde{\delta}_{k-1},m-h\rangle\\
& \ \ \ \ \ +\bar{\alpha}_k^2\cdot \|m-h\|^2+\bar{\alpha}_k^2\cdot c^2(p,q)\cdot \|h\|_{EL^1}^2\nonumber\\
&\leq \alpha_k^2\cdot \|\tilde{\delta}_{k-1}\|^2+ 2\alpha_k\bar{\alpha}_k\cdot \|\tilde{\delta}_{k-1}\|\|m-h\|+\bar{\alpha}_k^2\cdot \|m-h\|^2 \nonumber
\\
& \ \ \ \ +\bar{\alpha}_k^2\cdot c^2(p,q)\cdot \|h\|_{EL^1}^2\nonumber\\
&\leq \alpha_k\cdot \|\tilde{\delta}_{k-1}\|^2+\bar{\alpha}_k\cdot \|m-h\|^2+\bar{\alpha}_k^2\cdot c^2(p,q)\cdot \|h\|_{EL^1}^2\nonumber,
\end{align}

where in the second inequality the Caucy-Schwarz inequality is used. Taking expectation w.r.t. $\Xi_{k-1}$ on both sides of (<ref>) gives

\begin{equation}\label{RGAproof2d}
\E_{\Xi_k}(\|\tilde{\delta}_k\|^2-\|m-h\|^2)\leq \alpha_k\cdot \E_{\Xi_{k-1}}(\|\tilde{\delta}_{k-1}\|^2-\|m-h\|^2)+ \bar{\alpha}_k^2\cdot c^2(p,q)\cdot \|h\|_{EL^1}^2.
\end{equation}

It is easy to check that the result for rPPR.rga  holds obviously for $k=0$. Assuming that the result for rPPR.rga holds for $k-1$ where  $k\ge 1$, by (<ref>) we  also have

\begin{align}
\E_{\Xi_k}(\|\tilde{\delta}_k\|^2-\|m-h\|^2)& \leq \left(1-\frac{2}{k+2}\right)\cdot \frac{4c^2(p,q)\cdot \|h\|_{EL^1}^2}{k}+ \frac{4c^2(p,q)\cdot \|h\|_{EL^1}^2}{(k+2)^2}\nonumber\\
&\leq \frac{4c^2(p,q)\cdot \|h\|_{EL^1}^2}{k+1}\nonumber.
\end{align}

By mathematical induction, we have completed the proof  for rPPR.rga.    \(\Box\)

§.§ Two main technical results for the proofs of Theorems <ref>–<ref>

Let the sequence $t_n\asymp \ln{n}$. Assume
\E(e^{c\cdot Y^2})<\infty
Then, for any $h \in \mathcal{H}_{\infty,q_0,n}$ and and bounded $m(x)$ with $x\in\B^p$, we have

\begin{equation} \label{aaa00}
\E_{\mathcal{D}_n, \Xi_{k_n}} (\I(\hat{m}_{k_n,b}(x))) \leq   \frac{4\left(c\cdot \|h\|_{EL^1}\right)^2}{(k_n+1)}
+ \E\left(m(X)-h(X)\right)^2+ \omega^1_n;
\end{equation}

and for any  $ \tau \in\mathbb{Z}_+$ (i.e. number of iterations)  and $ \nu > 0 $, we have

\begin{equation}\label{Thesecondmainlemma}
\E_{\mathcal{D}_{n},\Xi_{\infty}}\left( \I(\hat{m}_{k_b^*,b}(x))\cdot \mathbb{I}(A_n)\right)\leq  \frac{ 8\left(c \cdot \|h\|_{EL^1}\right)^2}{\tau+1} + c\cdot \E\left(m(X)-h(X)\right)^2+ \omega_n^2,
\end{equation}

where $ \omega^1_n = {c n^{-1} } pk_n J_n\ln(pk_n J_n) \cdot \ln^3 n $ and $ \omega_n^2 = 2\tau\cdot {n^{-1} \ln{n}\cdot (q+J_n^{1+\nu})}
+c\cdot J_n^{-\nu}\cdot\ln^4 n$, and $A_n:=\{\max_{1\leq i\leq n}|Y_i|\leq t_n\}$ and  $\mathbb{I}(A_n)$ is the indicator of  $A_n$.

Proofs for (<ref>) and (<ref>) are similar;
we here only give the detail for former while that for the latter is given in . The proof of Lemma <ref> is based on Theorem 1 in [3] (see also Lemma <ref> in ), which showed

\begin{equation}\label{Mainlemmaformula}
\begin{aligned}
\E_{\mathcal{D}_n, \Xi_{k_n}} (\I(\hat{m}_{k_n,b}(x))) &\leq 2 \E_{\mathcal{D}_{n},\Xi_{k_n}}\left( \|{m}_{k_n, b}(X)-Y\|_n^2-\|m(X)-Y\|_n^2\right)\\
& \ \ \ +\frac{c\ln^2 n}{n}\cdot\sup_{x_1^n}\ln\left(\mathcal{N}(1/(80nt_n),\mathcal{H}^{t_n}_{k_n,p,n},x_1^n)\right),
\end{aligned}
\end{equation}

where $ {m}_{k_n, b}$ is defined in Section <ref> and $\mathcal{N}$ is the covering number of $\mathcal{H}^{t_n}_{k_n,p,n}:=\max\{\min\{\mathcal{H}_{k_n,p,n},t_n\},-t_n\}$. To give a  bound for the first part on the RHS of (<ref>), let

\begin{align*}
\Delta &=\|\tilde{m}_{k_n,b}(X)-Y\|_n^2-\|m(X)-Y\|_n^2 \\
&= (\|Y-\tilde{m}_{k_n,b}\|_{n}^2-\|Y-h\|_{n}^2) + (\|Y-h\|_{n}^2-\|Y-m\|_{n}^2) := \text{I} + \text{II}.
\end{align*}

For $\text{I}$, following the steps in proving Theorem <ref> with mean squared norm $\|\cdot\|$ replaced by the empirical norm $\|\cdot\|_n$, we have

\begin{equation}\label{Sta1}
\E_{\Xi_{k_n}|\mathcal{D}_{n}}(\text{I})\leq 4\left(c\cdot \|h\|_{n,EL^1}\right)^2\cdot(k_n+1)^{-1}.
\end{equation}

Note that, by Lemma 3.4 of [4],  $\E_{\mathcal{D}_{n}}(\|h\|_{n,EL^1}^2)\leq \|h\|_{EL^1}^2$. Thus, by taking expectation on both sides of (<ref>) we have

\begin{equation}\label{Sta2}
\E_{\Xi_{k_n}, \mathcal{D}_{n}}(\text{I})\leq 4\left(c\cdot \|h\|_{EL^1}\right)^2\cdot(k_n+1)^{-1}.
\end{equation}

It is easy to see that $\text{II}$ is independent of $\Xi_{k_n}$, thus we further have

\begin{equation}\label{Sta3}
\E_{\Xi_{k_n}, \mathcal{D}_{n}}(\text{II})=\E\left((Y-h(X))^2- (Y-m(X))^2\right)=\E\left(m(X)-h(X)\right)^2.
\end{equation}

Combination of (<ref>) and (<ref>) yields

\begin{equation}\label{Sta4}
\E_{\Xi_{k_n}, \mathcal{D}_{n}}(\Delta)\leq
4\left(c \cdot \|h\|_{EL^1}\right)^2\cdot(k_n+1)^{-1} + \E\left(m(X)-h(X)\right)^2.
\end{equation}

For the second part on the RHS of (<ref>),  following the proof of (<ref>) we have

\begin{equation}\label{Sta5}
\sup_{ x_1^n}  \mathcal{N}(1/(80n t_n),\mathcal{H}^{t_n}_{k_n,p,n},x_1^n) \leq  3\left( 480en t_n \right)^{c\cdot pk_nJ_n\ln(pk_nJ_n)},
\end{equation}

where ${\mathcal{H}^{t_n}_{k_n,p,n}}=\{\max\{\min\{g,t_n\},-t_n\}: g\in\mathcal{H}_{k_n,p,n}\}$ and $c>0$. It follows from (<ref>) and (<ref>) that

\begin{align}\label{Sta6}
%  \MSE(\hat{m}_{k_n,b}(x))
\E_{\mathcal{D}_n, \Xi_{k_n}} (\I(\hat{m}_{k_n,b}(x)))  &\leq  \frac{c\ln^2 n}{n}pk_nJ_n\ln(pk_nJ_n) \ln\left(n t_n \right) +\E_{\Xi_{k_n}, \mathcal{D}_{n}}(\Delta).
\end{align}

Lemma <ref> follows from (<ref>) and   (<ref>).

§.§ Proof of Theorem  <ref>

By (<ref>), there is a series of functions
$$h_n=\sum_{i=1}^L\sum_{j=1}^{K_n}{S_{n,i,j}(\theta_j^\top  x_{\C_j})}\in \mathcal{H}_{\infty,q_0,n},$$

where $S_{n,i,j}\in {\cal S}_{u_n,s}[-1,1]$ and $K_n=J_n^{q_0-1}$, such that
\|m-h_n\|^2\leq c\cdot LJ_n^{-2r}
$. As $J_n\to\infty$ when $n\to\infty$, for given any $\epsilon>0$ there is $n_1(\epsilon)\in\mathbb{Z}_+$ and $h_{n_1(\epsilon)}\in \mathcal{H}_{\infty,q_0,n}$ such that for any $n\ge n_1(\epsilon)$

\begin{equation*}\label{dfngmdfm}
\text{I}:= \E(m(X) - h_{n_1(\epsilon)}(X))^2 \le c\cdot LJ_{n_1(\epsilon)}^{-2r}\le \epsilon/3.
\end{equation*}

Now fix $n_1(\epsilon)$, thus $\|h_{n_1(\epsilon)}\|_{EL^1}$ is fixed. Since $ k_n \to \infty $, there is $n_2(\epsilon)\in\mathbb{Z}_+$ such that

\text{II}:= 4\left(c\cdot \|h_{n_1(\epsilon)}\|_{EL^1}\right)^2\cdot(k_n+1)^{-1} \le \epsilon/3

for any $n\ge n_2(\epsilon)$. By (<ref>) of Lemma <ref>, there is $n(\epsilon) \ge max(n_1(\epsilon), n_2(\epsilon))$ such that for any $n>n(\epsilon)$

\begin{equation*}
\MSE(\hat{m}_{k_n,b}(x)) \leq   \text{I} + \text{II} + \omega^1_n \le \epsilon
\end{equation*}

which implies

\begin{equation} \label{dfknge}
\MSE(\hat{m}_{k_n,b}(x)) \to 0,\ \text{as}\ n\to\infty.
\end{equation}

Note that, by Jensen's inequality,

\begin{equation}\label{JensenAB}
|\hat{m}^{ens}_{k_n}(x)-m(x)|^2\leq \frac{1}{B}\sum_{b=1}^B{| \hat{m}_{k_n,b}(x) - m(x)|^2}.
\end{equation}

\begin{equation}\label{JensenAB1}
\MSE(\hat{m}^{ens}_{k_n}(x))\leq \frac{1}{B}\sum_{b=1}^B \MSE(\hat{m}_{k_n,b}(x)).
\end{equation}

As a consequence, Theorem  <ref> follows from (<ref>) and (<ref>).  \(\Box\)

§.§ Proof of Theorem  <ref>

Recall that
m(x)=\sum_{\iota=1}^\infty {g_\iota(\theta_\iota^\top  x_{\C_\iota})}
where $\theta_\iota\in\mathbb{R}^{q_0}$ and $ card(\C_\iota)= q_0\leq q $. By (<ref>), there is a function $S_{n,\iota}\in {\cal S}_{u,s}[-1,1]$ such that

\| g_\iota-S_{n,\iota}\|\leq c(r)\cdot \|g_
\iota\|_{H^r(\B^p)}\cdot J_n^{-r}.

Therefore, we can find a series of functions $h_n(x)=\sum_{\iota=1}^\infty {S_{n,\iota}(\theta_\iota^\top x_{\C_\iota})}\in \mathcal{H}_{\infty,q_0,n}$, such that

\|m-h_n\|^2 \leq c\cdot a^2\cdot J_n^{-2r},

where $a=\sum_{\iota=1}^\infty{\|g_\iota\|_{H^r(\B^p)}}<\infty$ by the assumption of the theorem.

By (<ref>) of Lemma <ref> and $\|h_n\|_{EL^1}\le ca\cdot J_n^{-r}+a$ for each $n\in\mathbb{Z}_+$, we have

\begin{equation}\label{Theorem4part1}
\MSE(\hat{m}_{k_n,b}(x)) \leq   c \cdot(k_n+1)^{-1}
+ (ca)^2\cdot J_n^{-2r} + \omega^1_n.
\end{equation}

Properly selecting $k_n, J_n$ such that $k_n^{-1}=J_n^{-2r}=\frac{k_nJ_n}{n}$, namely $J_n=n^{\frac{1}{4r+1}}$, the RHS of (<ref>) is bounded by $c\cdot n^{-\frac{2r}{4r+1}}\cdot\ln^4 n$.

Finally, Theorem  <ref> follows from (<ref>) and (<ref>).   \(\Box\)

§.§ Proof of Theorem <ref>

Recalling $A_n $ in Lemma <ref>, let $ A_n^c $ be the complement. Write

\begin{align}
\I(\hat{m}_{k_b^*,b}(x))
&=\I(\hat{m}_{k_b^*,b}(x))\cdot \mathbb{I}(A_n^c)+ \I(\hat{m}_{k_b^*,b}(x)))\cdot \mathbb{I}(A_n) := \text{I}+\text{II}. \nonumber%\label{Theorem5part1}
\end{align}

For any given $\epsilon>0$, we have

\begin{equation} \label{Theorem5part2}
\begin{aligned}
\P\left( \I(\hat{m}_{k_b^*,b}(x))>\epsilon\right)&\leq \P(\text{I}>\epsilon/2)+\P(\text{II}>\epsilon/2)\\
&\leq \P(A_n^c)+\P(\text{II}>\epsilon/2):=\text{III}+\text{IV}.
\end{aligned}
\end{equation}

First, we consider part III. By the assumption on $ Y $, we have

\begin{align}
\label{subgaussguanjian}    \P(A_n^c)&=1- \P\left(\max_{1\leq i\leq n}|Y_i|\leq t_n\right)
=1-\left[\P(|Y_i|\leq t_n)\right]^n
\leq 1-(1-c\cdot e^{-c\cdot t_n^2})^n \\
&\le 1-e^{n\cdot \ln(1-c\cdot e^{-c\cdot t_n^2})}. \nonumber
\end{align}

\text{III} = \P(A_n^c) \to 0.

Next, consider part IV using (<ref>) of Lemma <ref>.

Choose the same $h_n$ in the proof of Theorem <ref>. For any given $\eta>0$, there is $n_1(\eta)\in\mathbb{Z}_+$ such that for any $n\ge n_1(\eta)$,

\begin{equation}\label{adfgarohr}
\text{V}:=\E\left(m(X)-h_{n_1(\eta)}(X)\right)^2 \leq c\cdot LJ_{n_1(\eta)}^{-2r} \le \eta/2.
\end{equation}

Now, fix $ n_1(\eta) $ for the given $ \eta $. Choose $\tau_n=n^\frac{1}{2}\cdot {J_n^{-\frac{1+\nu}{2}}}$ which goes to infinity as $n\to\infty$, i.e.
\omega_n^2 \to 0.

Therefore, there is $n_2(\eta)\in\mathbb{Z}_+$ such that for any $n\ge n_2(\eta)$

\begin{equation}\label{adfgarohr2}
\text{VI}:=  8\left(c\cdot \|h_{n_1(\eta)}\|_{EL^1}\right)^2\cdot(\tau_n+1)^{-1}+ \omega_n^2\le \eta/2.
\end{equation}

Combination of (<ref>) and (<ref>) leads to

\begin{equation*}
\E_{\mathcal{D}_{n},\Xi_{\infty}}\left( \I(\hat{m}_{k_b^*,b}(x)) \cdot \mathbb{I}(A_n)\right) \leq  \text{V}+\text{VI}\leq \eta
\end{equation*}

for any $n\ge \max\{n_1(\eta),n_2(\eta)\}$, which implies that for every $ b $,

\begin{equation}\label{Theorem5lastone}
\lim_{n\to\infty}{\E_{\mathcal{D}_{n},\Xi_{\infty}} \left(\I(\hat{m}_{k_b^*,b}(x))\cdot\mathcal{I}(A_n)\right)}=0.
\end{equation}

By  Chebyshev inequality, equation (<ref>) implies

\text{IV}=\P(\text{II}>\epsilon/2)\to 0.

In conclusion, we have proved that the LHS of (<ref>) converges to $0$. On the other hand, by Jensen's inequality, we have

\begin{equation}\label{JensenCD}
|\hat{m}^{ens}_{*}(x)-m(x)|^2\leq \frac{1}{B}\sum_{b=1}^B{| \hat{m}_{k_b^*,b}(x) - m(x)|^2}.
\end{equation}

Theorem <ref> follows from  (<ref>) and (<ref>).   \(\Box\)

§.§ Proof of Theorem <ref>

Let $A_n$ and $A_n^c $ be defined in the proof of Theorem <ref> and $a_n=n^{-\frac{2r'}{4r'+1+\nu}}\cdot \ln^4 n$ for each $n\in\mathbb{Z}_+$. Write

\begin{equation}\label{Theorem7part1}
\frac{1}{a_n}\cdot\I(\hat{m}_{k_b^*,b}(x))=\frac{1}{a_n}\cdot\I(\hat{m}_{k_b^*,b}(x)\cdot \mathbb{I}(A_n^c)
+ \frac{1}{a_n}\cdot\I(\hat{m}_{k_b^*,b}(x))\cdot \mathbb{I}(A_n) := \text{I}+\text{II}.
\end{equation}

For any given $\delta>0$  and $M(\delta)>0$, we have

\begin{equation}\label{Theorem7part2}
\begin{aligned}
\P\left(\frac{1}{a_n}\cdot\I(\hat{m}_{k_b^*,b}(x))>M(\delta)\right)&\leq\P(\text{I}>M(\delta)/2)
&\leq \P(A_n^c)+\P(\text{II}>M(\delta)/2):=\text{III}+\text{IV}.
\end{aligned}
\end{equation}

For III, we can use the same argument used in the proof of Theorem  <ref> to show that there is $n_1(\delta)\in\mathbb{Z}_+$ such that  $\text{III} < \delta/2$ for all $n\geq n_1(\delta)$.

Next, we consider part IV. Choose the same $h_n$ in the proof of Theorem <ref> satisfying

\|m-h_n\|^2 \leq c\cdot J_n^{-2r}.

By (<ref>) of Lemma <ref> and the fact that $\|h_n\|_{EL^1}\le ca\cdot J_n^{-r}+a$ for each $n\in\mathbb{Z}_+$, we have

\begin{equation}\label{Theorem7part1}
\E_{\mathcal{D}_{n},\Xi_{\infty}}\left( \I(\hat{m}_{k_b^*,b}(x)\cdot \mathbb{I}(A_n)\right)\leq   c\cdot(\tau_n+1)^{-1}
+ c\cdot J_n^{-2r}+ \omega_n^2
\end{equation}

for any $\tau_n \in\mathbb{Z}_+$. Properly selecting $\tau_n,J_n$ such that $\tau_n^{-1}=J_n^{-2r}=\tau_n\cdot J_n^{1+\nu}\cdot n^{-1}$, namely $J_n=n^{\frac{1}{4r+1+\nu}}$ and $\tau_n=n^\frac{1}{2}\cdot J_n^{-\frac{1+\nu}{2}}\to\infty$ as long as $ 2r\leq \nu$, the LHS of  (<ref>) is bounded by the product of a constant and $a_n=n^{-\frac{2r'}{4r'+1+\nu}}\cdot \ln^4 n$. By Chebyshev inequality again, there is $M_1(\delta)>0$ such that

\text{IV}:=\P(\text{II}>M_1(\delta)/2)\leq \frac{2\cdot \E_{\mathcal{D}_{n},\Xi_{\infty}}(\text{II})}{M_1(\delta)}\leq \delta/2

holds for all $n\in\mathbb{Z}_+$.

Therefore, by (<ref>) we have shown that

\begin{equation}\label{Aldfng}
\P\left(\frac{1}{a_n}\cdot\I(\hat{m}_{k_b^*,b}(x))>M_1(\delta)\right)\le \text{III}+\text{IV}\le \delta
\end{equation}

holds for each $n\ge n_1(\delta)$. Theorem <ref> follows from (<ref>) and (<ref>).  \(\Box\)

§ DISCUSSION

As commented by [32], “PPR does represent an important intellectual advance, one that has blossomed in its reincarnation in the field of neural networks”. [32] attributed the low popularity of PPR to the computational issue, which is part of the motivation for this paper. This paper studied three greedy algorithms of PPR and proposed an ensemble procedure, ePPR, based on “feature boosting”. Extensive numerical study suggests that ePPR can generally produce obviously more accurate predictions than the popular methods, including RF and its variants, for data with either continuous response or categorical response.

One possible reason why ePPR has an advantage over RF is in their respective underlying models, PPR and CART. The former can approximate any functions as shown in this paper, while the latter cannot because of its splitting is based on marginal variables. Moreover, PPR is continuous but CART is a step function which may cause problems in prediction ([52]). The main problem with ePPR is the computational burden. However, the ensemble procedure can be easily reduce the calculation burden by using parallel computation. On the other hand, we believe there are more connections between ePPR and RF or ANN, which need further investigations.

SUPPLEMENTARY MATERIAL

Supplement to: "Ensemble Projection Pursuit  for General Nonparametric Regression"
provides the auxiliary and other relevant results that are used in the proofs in Section <ref>.

[1]
Adams, Robert AR. A. Fournier, John JFJ. J.
Sobolev spaces.

[2]
Athey, SusanS.,
Tibshirani, JulieJ. Wager, StefanS.
Generalized random forests.
The Annals of Statistics

[3]
Bagirov, Adil MA. M.,
Clausen, ConnyC. Kohler, MichaelM.
Estimation of a regression function by maxima of minima of linear
IEEE Transactions on Information Theory

[4]
Barron, Andrew RA. R.,
Cohen, AlbertA.,
Dahmen, WolfgangW. DeVore, Ronald AR. A.
Approximation and learning by greedy algorithms.
The annals of statistics

[5]
Bergmeir, ChristophC. Benítez, José M.J. M.
Neural Networks in R Using the Stuttgart Neural Network Simulator:
Journal of Statistical Software

[6]
Braun, AlinaA.,
Kohler, MichaelM. Walk, HarroH.
On the rate of convergence of a neural network regression estimate
learned by gradient descent.
arXiv preprint arXiv:1912.03921.

[7]
Breiman, L.L.
Classification and Regression Trees Regression Trees.
Wadsworth/Thomson Learning.

[8]
Breiman, LeoL.
Statistical Modeling: The Two Cultures (with comments and a rejoinder
by the author).
Statistical Science
199 – 231.

[9]
Breiman, LeoL.
Random forests.
Machine learning

[10]
Breiman, LeoL.,
Friedman, JeromeJ.,
Stone, Charles J.C. J. Olshen, R. A.R. A.
Classification and Regression Trees.
Chapman and Hall/CRC.

[11]
Brenner, Susanne CS. C.,
Scott, L RidgwayL. R. Scott, L RidgwayL. R.
The mathematical theory of finite element methods

[12]
Brown, G.G.,
Wyatt, J.J.,
Harris, R.R. Yao, XX.
Diversity creation methods: a survey and categorisation.
Information Fusion

[13]
Caruana, RichR.,
Karampatziakis, NikosN. Yessenalina, AinurA.
An Empirical Evaluation of Supervised Learning in High Dimensions.
In Proceedings of the 25th International Conference on Machine
ICML '08
Association for Computing Machinery, New York, NY, USA.

[14]
Caruana, RichR. Niculescu-Mizil, AlexandruA.
An Empirical Comparison of Supervised Learning Algorithms.
In Proceedings of the 23rd international conference on Machine

[15]
Chen, HungH.
Estimation of a projection-pursuit type regression model.
Annals of Statistics

[16]
Chen, TianqiT.,
He, TongT.,
Benesty, MichaelM.,
Khotilovich, VadimV.,
Tang, YuanY.,
Cho, HyunsuH.,
Chen, KailongK. et al.
Xgboost: extreme gradient boosting.
R package version 0.4-2

[17]
Chi, Chien-MingC.-M.,
Vossler, PatrickP.,
Fan, YingyingY. Lv, JinchiJ.
Asymptotic Properties of High-Dimensional Random Forests.
arXiv preprint arXiv:2004.13953.

[18]
Cortes, CorinnaC. Vapnik, VladimirV.
Support-vector networks.
Machine learning

[19]
Culp, MarkM.,
Johnson, KjellK. Michailides, GeorgeG.
ada: An r package for stochastic boosting.
Journal of statistical software

[20]
Davis, Geoffrey MG. M.,
Mallat, Stephane GS. G. Zhang, ZhifengZ.
Adaptive time-frequency decompositions with matching pursuit.
In Wavelet Applications

[21]
DeVore, Ronald AR. A. Lorentz, George GG. G.
Constructive approximation
Springer Science & Business Media.

[22]
Dimitriadou, EvgeniaE.,
Hornik, KurtK.,
Leisch, FriedrichF.,
Meyer, DavidD. Weingessel, AndreasA.
Misc functions of the Department of Statistics (e1071), TU Wien.
R Package

[23]
Friedberg, RinaR.,
Tibshirani, JulieJ.,
Athey, SusanS. Wager, StefanS.
Local linear forests.
Journal of Computational and Graphical Statistics

[24]
Friedman, JeromeJ.,
Hastie, TrevorT. Tibshirani, RobertR.
Regularization Paths for Generalized Linear Models via Coordinate
Journal of Statistical Software

[25]
Friedman, Jerome HJ. H.
Multivariate adaptive regression splines.
Annals of Statistics

[26]
Friedman, Jerome HJ. H. Stuetzle, WernerW.
Projection pursuit regression.
Journal of the American statistical Association

[27]
Friedman, J. H.J. H. Tukey, J. W.J. W.
A projection pursuit algorithm for exploratory data analysis.
IEEE Transactions on computers

[28]
Genuer, R.R. Poggi, J. M.J. M.
Random Forests with R.
Use R!
Springer International Publishing.

[29]
Giné, E.E. Nickl, R.R.
Mathematical Foundations of Infinite-Dimensional Statistical Models.
Cambridge Series in Statistical and Probabilistic Mathematics.
Cambridge University Press.

[30]
Györfi, LászlóL.,
Kohler, MichaelM., Krzyżak, AdamA. Walk, HarroH.
A distribution-free theory of nonparametric regression

[31]
Härdle, WolfgangW.,
Hall, PeterP. Ichimura, HidehikoH.
Optimal Smoothing in Single-Index Models.
The Annals of Statistics
157 – 178.

[32]
Hastie, TrevorT.,
Tibshirani, RobertR. Friedman, JeromeJ.
The Elements of Statistical Learning.
Springer Series in Statistics.
Springer New York Inc., New York, NY, USA.

[33]
Hastie, Trevor JT. J. Tibshirani, Robert JR. J.
Generalized additive models.
Chapman & Hall/CRC.

[34]
Ho, Tin KamT. K.
Random decision forests.
In Document analysis and recognition, 1995., proceedings of the
third international conference on

[35]
Hristache, MarianM.,
Juditsky, AnatoliA. Spokoiny, VladimirV.
Direct estimation of the index coefficient in a single-index model.
The Annals of Statistics
593 – 623.

[36]
Huber, Peter JP. J.
Projection Pursuit.
Annals of Statistics

[37]
Jones, Lee KL. K.
On a conjecture of Huber concerning the convergence of projection
pursuit regression.
The annals of statistics

[38]
Jones, Lee K.L. K.
A Simple Lemma on Greedy Approximation in Hilbert Space and Convergence
Rates for Projection Pursuit Regression and Neural Network Training.
The Annals of Statistics

[39]
Kidger, PatrickP. Lyons, TerryT.
Universal approximation with deep narrow networks.
In Conference on learning theory

[40]
Klusowski, Jason MJ. M.
Universal Consistency of Decision Trees in High Dimensions.
arXiv preprint arXiv:2104.13881.

[41]
Kruskal, J. B.J. B.
Toward a practical method which helps uncover the structure of a set of
multivariate observations by finding the linear transformation which
optimizes a new index of condensation.
In Statistical Computation.
Academic Press.

[42]
Kuncheva, L.L. Whitaker, C.C.
Measures of diversity in classifier ensembles.
Machine Learning

[43]
Liaw, AndyA.,
Wiener, MatthewM. et al.
Classification and regression by randomForest.
R news

[44]
Pati, Yagyensh ChandraY. C.,
Rezaiifar, RaminR. Krishnaprasad, Perinkulam SambamurthyP. S.
Orthogonal matching pursuit: Recursive function approximation with
applications to wavelet decomposition.
In Proceedings of 27th Asilomar conference on signals, systems and

[45]
Petrushev, Pencho PP. P.
Approximation by ridge functions and neural networks.
SIAM Journal on Mathematical Analysis

[46]
Rashid, A. AlievA. A. Vugar, E. IsmailovE. I.
A representation problem for smooth sums of ridge functions.
Journal of Approximation Theory

[47]
Sande, EspenE.,
Manni, CarlaC. Speleers, HendrikH.
Explicit error estimates for spline approximation of arbitrary
smoothness in isogeometric analysis.
Numerische Mathematik

[48]
Scornet, ErwanE.,
Biau, GérardG. Vert, Jean-PhilippeJ.-P.
Consistency of random forests.
The Annals of Statistics

[49]
Shalev-Shwartz, S.S. Ben-David, S.S.
Understanding Machine Learning: From Theory to Algorithms.
Understanding Machine Learning: From Theory to Algorithms.
Cambridge University Press.

[50]
Stone, Charles J.C. J.
Optimal Global Rates of Convergence for Nonparametric Regression.
The Annals of Statistics
1040 – 1053.

[51]
Tan, XinX. Zhan, HaoranH.
Estimation of Projection Pursuit Regression via Alternating

[52]
Wager, StefanS. Walther, GuentherG.
Adaptive concentration of regression trees, with application to random
arXiv preprint arXiv:1503.06388.

[53]
Xia, YingcunY.,
Tong, HowellH. Li, Wai KeungW. K.
On extended partially linear single-index models.

[54]
Zhang, HepingH. Wang, MinghuiM.
Search for the smallest random forest.
Statistics and its interface

[55]
Zhu, RuoqingR.,
Zeng, DonglinD. Kosorok, Michael R.M. R.
Reinforcement Learning Trees.
Journal of the American Statistical Association