Additive Models with Trend Filtering

TL;DR

This paper introduces additive trend filtering, a method for modeling additive functions with components that are piecewise polynomial, offering theoretical optimality, computational efficiency, and a novel parallelizable backfitting algorithm.

Contribution

It develops a new additive trend filtering approach with proven minimax optimal error rates and introduces a parallelizable backfitting algorithm for efficient computation.

Findings

Achieves minimax optimal error rates for additive functions with bounded variation derivatives.

Demonstrates superior theoretical and computational properties over additive smoothing splines.

Introduces a novel parallel backfitting algorithm for additive trend filtering.

Abstract

We study additive models built with trend filtering, i.e., additive models whose components are each regularized by the (discrete) total variation of their $k$th (discrete) derivative, for a chosen integer $k \geq 0$. This results in $k$th degree piecewise polynomial components, (e.g., $k=0$ gives piecewise constant components, $k=1$ gives piecewise linear, $k=2$ gives piecewise quadratic, etc.). Analogous to its advantages in the univariate case, additive trend filtering has favorable theoretical and computational properties, thanks in large part to the localized nature of the (discrete) total variation regularizer that it uses. On the theory side, we derive fast error rates for additive trend filtering estimates, and show these rates are minimax optimal when the underlying function is additive and has component functions whose derivatives are of bounded variation. We also show that these rates are unattainable by additive smoothing splines (and by additive models built from linear smoothers, in general). On the computational side, as per the standard in additive models, backfitting is an appealing method for optimization, but it is particularly appealing for additive trend filtering because we can leverage a few highly efficient univariate trend filtering solvers. Going one step further, we describe a new backfitting algorithm whose iterations can be run in parallel, which (as far as we know) is the first of its kind. Lastly, we present experiments to examine the empirical performance of additive trend filtering.

Figures

Equations

$$Y^i=\mu +\sum _{j=1}^d{f}_{0j}(X_j^i)+{ϵ}^i,i=1,\dots ,n,$$

$$\displaystyle$\displaystyle\min_{\theta_{1},\ldots,\theta_{d} \in\mathbb{R}^{n}}$\hfil\ignorespaces$$

$$\displaystyle$\displaystyle\mathop{\mathrm{subject\,\,to}}$\hfil\ignorespaces$$

$${\theta }_j^i=f_j(X_j^i),i=1,\dots ,n,j=1,\dots ,d.$$

$$\theta_{j}=(I+\lambda Q_{j})^{-1}\bigg{(}Y-\bar{Y}\mathbbm{1}-\sum_{\ell\not=j}\theta_{\ell}\bigg{)},\;\;\;j=1,\ldots,d,$$

$$\theta_{j}=S_{j}\bigg{(}Y-\bar{Y}\mathbbm{1}-\sum_{\ell\not=j}\theta_{\ell}\bigg{)},\;\;\;j=1,\ldots,d,$$

$${\mathcal{F}}_k(C)=\{f:\mathrm{TV}(f^{(k)})\le C\},$$

$$\underset{\theta \in {\mathbb{R}}^n}{\min }\frac{1}{2}\Vert Y-\theta {\Vert }_2^2+\lambda \Vert D^{(X,k+1)}\theta {\Vert }_1,$$

$$D^{(X,1)}$$

$$D^{(X,k+1)}$$

$$\displaystyle$\displaystyle\min_{\theta_{1},\ldots,\theta_{d} \in\mathbb{R}^{n}}$\hfil\ignorespaces$$

$$\displaystyle$\displaystyle\mathop{\mathrm{subject\,\,to}}$\hfil\ignorespaces$$

$$S_j{X}_j=(X_j^{(1)},X_j^{(2)},\dots ,X_j^{(n)}),j=1,\dots ,d.$$

$$f_{01}(t)=\min(t,1-t)^{0.2}\sin\bigg{(}\frac{2.85\pi}{0.3+\min(t,1-t)}\bigg{)},\;\;\;f_{02}(t)=e^{3t}\sin(4\pi t),\;\;\;f_{03}(t)=-(t-1/2)^{2},$$

$$\sum_{j=1}^{d}\lambda_{j}\theta_{j}^{T}Q_{j}\theta_{j}\;\;\;\text{and}\;\;\;\sum_{j=1}^{d}\lambda_{j}\big{\|}D^{(X_{j},k+1)}S_{j}\theta_{j}\big{\|}_{1},$$

$$\begin{array}{c}{h}_i(t)=\prod _{\ell =1}^{i-1}(t-t^{\ell }),i=1,\dots ,k+1,\\ h_{i+k+1}(t)=\prod _{\ell =1}^k(t-t^{i+\ell })\cdot 1\{t>t^{i+k}\},i=1,\dots ,n-k-1.\end{array}$$

$$g=\sum _{i=1}^n{\alpha }^i{h}_i$$

$$\displaystyle$\displaystyle\min_{\alpha_{1},\ldots,\alpha_{d} \in\mathbb{R}^{n}}$\hfil\ignorespaces$$

$$\displaystyle$\displaystyle\mathop{\mathrm{subject\,\,to}}$\hfil\ignorespaces$$

$${\hat{\theta }}_j^i=\sum _{\ell =1}^n{\hat{\alpha }}_j^{\ell }{h}_{\ell }^{(X_j)}(X_j^i),i=1,\dots ,n,j=1,\dots ,d.$$

$$\displaystyle$\displaystyle\min_{f_{j} \in\mathcal{H}_{j}, \, j=1,\ldots,d}$\hfil\ignorespaces$$

$$\displaystyle$\displaystyle\mathop{\mathrm{subject\,\,to}}$\hfil\ignorespaces$$

$${\hat{\theta }}_j^i={\hat{f}}_j(X_j^i),i=1,\dots ,n,j=1,\dots ,d.$$

$$H_{i\ell }^{(X_j,k)}=h_{\ell }^{(X_j)}(X_j^i),i=1,\dots ,n,\ell =1,\dots ,n.$$

$$(S_{j}H^{(X_{j},k)})^{-1}=\left[\begin{array}[]{c}C^{(X_{j},k+1)}\\ \frac{1}{k!}D^{(X_{j},k+1)}\end{array}\right],$$

$$(H^{(X_{j},k)})^{-1}=\left[\begin{array}[]{c}C^{(X_{j},k+1)}\\ \frac{1}{k!}D^{(X_{j},k+1)}\end{array}\right]S_{j}.$$

$$f_j^{(k)}(t)=k!+k!\sum _{\ell =k+2}^n{\alpha }_j^{\ell }\cdot 1\{t>X_j^{\ell -1}\},$$

$${\hat{f}}_j(x_j)=\sum _{\ell =1}^n{\hat{\alpha }}_j^{\ell }{h}_{\ell }^{(X_j)}(x_j),x_j\in \mathbb{R}.$$

$${\hat{a}}_j=C^{(X_j,k+1)}{S}_j{\hat{\theta }}_j,$$

$${\hat{b}}_j=\frac{1}{k!}{D}^{(X_j,k+1)}{S}_j{\hat{\theta }}_j.$$

$$\tilde{P}=M\big{[}\,\tilde{P}^{(X_{1},k)}\;\;\ldots\;\;\tilde{P}^{(X_{d},k)}\,\big{]}\in\mathbb{R}^{n\times dk},$$

$$\tilde{K}=U^{T}M\big{[}\,K^{(X_{1},k)}\;\;\ldots\;\;K^{(X_{d},k)}\,\big{]}\in\mathbb{R}^{(n-kd-1)\times(n-k-1)d},$$

$${\hat{\theta }}_j={\hat{\theta }}_j^{\mathrm{poly}}+{\hat{\theta }}_j^{\mathrm{knot}},j=1,\dots ,d,$$

Full text

1 Introduction

As the dimension of the input space grows large, nonparametric regression turns into a notoriously difficult problem. In this work, we adopt the stance taken by many others, and consider an additive model for responses $Y^{i}\in\mathbb{R}$, $i=1,\ldots,n$ and corresponding input points $X^{i}=(X^{i}_{1},\ldots,X^{i}_{d})\in\mathbb{R}^{d}$, $i=1,\ldots,n$, of the form

[TABLE]

where $\mu\in\mathbb{R}$ is an overall mean parameter, each $f_{0j}$ is a univariate function with $\sum_{i=1}^{n}f_{0j}(X^{i}_{j})=0$ for identifiability, $j=1,\ldots,d$, and the errors $\epsilon^{i}$, $i=1,\ldots,n$ are i.i.d. with mean zero. A comment on notation: here and throughout, when indexing over the $n$ samples we use superscripts, and when indexing over the $d$ dimensions we use subscripts, so that, e.g., $X^{i}_{j}$ denotes the $j$th component of the $i$th input point. (Exceptions will occasionally be made, but the role of the index should be clear from the context.)

Additive models are a special case of the more general projection pursuit regression model of Friedman and Stuetzle (1981). Additive models for the Cox regression and logistic regression settings were studied in Tibshirani (1983) and Hastie (1983), respectively. Some of the first asymptotic theory for additive models was developed in Stone (1985). Two algorithms closely related to (backfitting for) additive models are the alternating least squares and alternating conditional expectations methods, from van der Burg and de Leeuw (1983) and Breiman and Friedman (1985), respectively. The work of Buja et al. (1989) advocates for the use of additive models in combination with linear smoothers, a surprisingly simple combination that gives rise to flexible and scalable multidimensional regression tools. The book by Hastie and Tibshirani (1990) is the definitive practical guide for additive models for exponential family data distributions, i.e., generalized additive models.

More recent work on additive models is focused on high-dimensional nonparametric estimation, and here the natural goal is to induce sparsity in the component functions, so that only a few select dimensions of the input space are used in the fitted additive model. Some nice contributions are given in Lin and Zhang (2006); Ravikumar et al. (2009); Meier et al. (2009), all primarily focused on fitting splines for component functions and achieving sparsity through a group lasso type penalty. In other even more recent and interesting work sparse additive models, Lou et al. (2016) consider a semiparametric (partially linear) additive model, and Petersen et al. (2016) consider a formulation that uses fused lasso (i.e., total variation) penalization applied to the component functions.

The literature on additive models (and by now, sparse additive models) is vast and the above is far form a complete list of references. In this paper, we examine a method for estimating additive models wherein each component is fit in a way that is locally adaptive to the underlying smoothness along its associated dimension of the input space. The literature on this line of work, as far as we can tell, is much less extensive. First, we review linear smoothers in additive models, motivate our general goal of local adaptivity, and then describe our specific proposal.

1.1 Review: additive models and linear smoothers

The influential paper by Buja et al. (1989) studies additive minimization problems of the form

[TABLE]

where $Y=(Y^{1},\ldots,Y^{n})\in\mathbb{R}^{n}$ denotes the vector of responses, and $Y-\bar{Y}\mathbbm{1}$ is its centered version, with $\bar{Y}=\frac{1}{n}\sum_{i=1}^{n}Y^{i}$ denoting the sample mean of $Y$, and $\mathbbm{1}=(1,\ldots,1)\in\mathbb{R}^{n}$ the vector of all 1s. Each vector $\theta_{j}=(\theta^{1}_{j},\ldots,\theta^{n}_{j})\in\mathbb{R}^{n}$ represents the evaluations of the $j$th component function $f_{j}$ in our model, i.e., tied together by the relationship

[TABLE]

In the problem (1), $\lambda\geq 0$ is a regularization parameter and $Q_{j}$, $j=1,\ldots,d$ are penalty matrices. As a typical example, we might consider $Q_{j}$ to be the Reinsch penalty matrix for smoothing splines along the $j$th dimension of the input space, for $j=1,\ldots,d$. Under this choice, a backfitting (block coordinate descent) routine for (1) would repeatedly cycle through the updates

[TABLE]

where the $j$th update fits a smoothing spline to the $j$th partial residual, over the $j$th dimension of the input points, denoted by $X_{j}=(X^{1}_{j},X^{2}_{j},\ldots X^{n}_{j})\in\mathbb{R}^{n}$. At convergence, we arrive at an additive smoothing spline estimate, which solves (1).

Modeling the component functions as smoothing splines is arguably the most common formulation for additive models, and it is the standard in several statistical software packages like the R package gam. As Buja et al. (1989) explain, the backfitting steps in (2) suggest that a more algorithmic approach to additive modeling can be taken. Instead of starting with a particular criterion in mind, as in (2), one can instead envision repeatedly cycling through updates

[TABLE]

where each $S_{j}$ is a particular (user-chosen) linear smoother, meaning, a linear map that performs a univariate smoothing across the $j$th dimension of inputs $X_{j}$. The linear smoothers $S_{j}$, $j=1,\ldots,d$ could correspond to, e.g., smoothing splines, regression splines (regression using a spline basis with given knots), kernel smoothing, local polynomial smoothing, or a combination of these, across the input dimensions. The convergence point of the iterations (3) solves a problem of the form (1) with $\lambda Q_{j}=S_{j}^{+}-I$, where $S_{j}^{+}$ is the Moore-Penrose pseudoinverse of $S_{j}$, for $j=1,\ldots,d$.

The class of linear smoothers is broad enough to offer fairly flexible, interesting mechanisms for smoothing, and simple enough to understand precisely. Buja et al. (1989) provide a unified analysis of additive models with linear smoothers, in which they derive the effective degrees of freedom of these estimators and a generalized cross-validation routine for tuning; they also study fundamental properties such as uniqueness of the component fits, and convergence of the backfitting steps.

Much of the work following Buja et al. (1989) remains in keeping with the idea of using linear smoothers in combination with additive models. Studying high-dimensional additive models, Lin and Zhang (2006); Ravikumar et al. (2009); Meier et al. (2009); Koltchinskii and Yuan (2010); Raskutti et al. (2012) all essentially build their methods off of linear smoothers, with modifications to induce sparsity in the estimated component functions. Ravikumar et al. (2009) consider a sparsified version of backfitting in (3), while the others consider penalized versions of the additive criterion in (1).

1.2 The limitations of linear smoothers

The beauty of linear smoothers lies in their simplicity. However, with this simplicity comes serious limitations, in terms of their ability to adapt to varying local levels of smoothness. In the univariate setting, the seminal theoretical work by Donoho and Johnstone (1998) makes this idea precise. With $d=1$, suppose that underlying regression function $f_{0}$ lies in the univariate function class

[TABLE]

for a constant $C>0$, where $\mathrm{TV}(\cdot)$ is the total variation operator, and $f^{(k)}$ the $k$th weak derivative of $f$. The class in (4) allows for greater fluctuation in the local level of smoothness of $f_{0}$ than, say, more typical function classes like Holder and Sobolev spaces. The results of Donoho and Johnstone (1998) (see also Section 5.1 of Tibshirani (2014)) imply that the minimax error rate for estimation over $\mathcal{F}_{k}(C)$ is $n^{-(2k+2)/(2k+3)}$, but the minimax error rate when we consider only linear smoothers (linear transformations of $Y$) is $n^{-(2k+1)/(2k+2)}$. This difference is highly nontrivial, e.g., for $k=0$ this is a difference of $n^{-2/3}$ (optimal) versus $n^{-1/2}$ (optimal among linear smoothers) for estimating a function $f_{0}$ of bounded variation.

It is important to emphasize that this shortcoming is not just a theoretical one; it is also clearly noticeable in basic practical examples. This does not bode well for additive models built from linear smoothers, when estimating component functions $f_{0j}$, $j=1,\ldots,d$ that display locally heterogeneous degrees of smoothness. Just as linear smoothers will struggle in the univariate setting, an additive estimate based on linear smoothers will not be able to efficiently track local changes in smoothness, across any of the input dimensions. This could lead to a loss in accuracy even if only some (or one) of the components $f_{0j}$, $j=1,\ldots,d$ possesses heterogeneous smoothness across its domain.

Two well-studied univariate estimators that are locally adaptive, i.e., that attain the minimax error rate over the $k$th order total variation class in (4), are wavelet smoothing and locally adaptive regression splines, as developed by Donoho and Johnstone (1998) and Mammen and van de Geer (1997), respectively. There is a substantial literature on these methods in the univariate case (especially for wavelets), but fewer authors have considered using these locally adaptive estimators in the additive models context. Some notable exceptions are Zhang and Wong (2003); Sardy and Tseng (2004), who study additive models built from wavelets, and Petersen et al. (2016), who study sparse additive models with components given by 0th order locally adaptive regression splines (i.e., the components are regularized via fused lasso penalties or total variation penalties). The latter work is especially related to our focus in this paper.

1.3 Additive trend filtering

We consider additive models that are constructed using trend filtering (instead of linear smoothers, wavelets, or locally adaptive regression splines) as their componentwise smoother. Proposed independently by Steidl et al. (2006) and Kim et al. (2009), trend filtering is a relatively new approach to univariate nonparametric regression. As explained in Tibshirani (2014), it can be seen as a discrete-time analog of the locally adaptive regression spline estimator. Denoting by $X=(X^{1},\ldots,X^{n})\in\mathbb{R}^{n}$ the vector of univariate input points, where we assume $X^{1}<\ldots<X^{n}$, the trend filtering estimate of order $k\geq 0$ is defined as the solution of the optimization problem

[TABLE]

where $\lambda\geq 0$ is a tuning parameter, and $D^{(X,k+1)}\in\mathbb{R}^{(n-k-1)\times n}$ is a $k$th order difference operator, constructed based on $X$. These difference operators can be defined recursively, as in

[TABLE]

(The leading matrix $D^{(X,1)}$ in (11) is the $(n-k-1)\times(n-k)$ version of the difference operator in (10).) Intuitively, the interpretation is that the problem (5) penalizes the sum of absolute $(k+1)$st order discrete derivatives of $\theta^{1},\ldots,\theta^{n}$ across the input points $X^{1},\ldots,X^{n}$. Thus, at optimality, the coordinates of the trend filtering solution $\hat{\theta}^{1},\ldots,\hat{\theta}^{n}$ obey a $k$th order piecewise polynomial form.

This intuition is formalized in Tibshirani (2014) and Wang et al. (2014), where it is shown that the components of the $k$th order trend filtering estimate $\hat{\theta}$ are precisely the evaluations of a fitted $k$th order piecewise polynomial function across the inputs, and that the trend filtering and locally adaptive regression spline estimates of the same order $k$ are asymptotically equivalent. When $k=0$ or $k=1$, in fact, there is no need for asymptotics, and the equivalence between trend filtering and locally adaptive regression spline estimates is exact in finite samples. It is also worth pointing out that when $k=0$, the trend filtering estimate reduces to the 1d fused lasso estimate (Tibshirani et al., 2005), which is known as 1d total variation denoising in signal processing (Rudin et al., 1992).

Over the $k$th order total variation function class defined in (4), Tibshirani (2014); Wang et al. (2014) prove that $k$th order trend filtering achieves the minimax optimal $n^{-(2k+2)/(2k+3)}$ error rate, just like $k$th order locally adaptive regression splines. Another important property, as developed by Kim et al. (2009); Tibshirani (2014); Ramdas and Tibshirani (2016), is that trend filtering estimates are relatively cheap to compute—much cheaper than locally adaptive regression spline estimates—owing to the bandedness of the difference operators in (10), (11), which means that specially implemented convex programming routines can solve (5) in an efficient manner.

It is this computational efficiency, along with its capacity for local adaptivity, that makes trend filtering a particularly desirable candidate to extend to the additive model setting. Specifically, we consider the additive trend filtering estimate of order $k\geq 0$, defined as a solution in the problem

[TABLE]

As before, $Y-\bar{Y}\mathbbm{1}$ is the centered response vector, and $\lambda\geq 0$ is a regularization parameter. Not to be confused with the notation for linear smoothers from a previous subsection, $S_{j}\in\mathbb{R}^{n\times n}$ in (12) is a permutation matrix that sorts the $j$th component of inputs $X_{j}=(X^{1}_{j},X^{2}_{j},\ldots X^{n}_{j})$ into increasing order, i.e.,

[TABLE]

Also, $D^{(X_{j},k+1)}$ in (12) is the $(k+1)$st order difference operator, as in (10), (11), but defined over the sorted $j$th dimension of inputs $S_{j}X_{j}$, for $j=1,\ldots,d$. With backfitting (block coordinate descent), computation of a solution in (12) is still quite efficient, since we can leverage the efficient routines for univariate trend filtering.

1.4 A motivating example

Figure 1 shows a simulated example that compares the additive trend filtering estimates in (12) (of quadratic order, $k=2$), to the additive smoothing spline estimates in (1) (of cubic order). In the simulation, we used $n=3000$ and $d=3$. We drew input points $X^{i}\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Unif}[0,1]^{3}$, $i=1,\ldots,3000$, and drew responses $Y^{i}\overset{\mathrm{i.i.d.}}{\sim}N(\sum_{j=1}^{3}f_{0j}(X^{i}_{j}),\sigma^{2})$, $i=1,\ldots,3000$, where $\sigma=1.72$ was set to give a signal-to-noise ratio of about 1. The underlying component functions were defined as

[TABLE]

so that $f_{01},f_{02},f_{03}$ possess different levels of smoothness ($f_{03}$ being the smoothest, $f_{02}$ less smooth, and $f_{01}$ the least smooth), and so that $f_{01}$ itself has heteregeneous smoothness across its domain.

The first row of Figure 1 shows the estimated component functions from additive trend filtering, at a value of $\lambda$ that minimizes the mean squared error (MSE), computed over 20 repetitions. The second row shows the estimates from additive smoothing splines, also at a value of $\lambda$ that minimizes the MSE. We see that the trend filtering fits adapt well to the varying levels of smoothness, but the smoothing spline fits are undersmoothed, for the most part. In terms of effective degrees of freedom (df), the additive smoothing spline estimate is much more complex, having about 85 df (computed via Monte Carlo over the 20 repetitions); the additive trend filtering has only about 42 df. The third row of the figure shows the estimates from additive smoothing splines, when $\lambda$ is chosen so that the resulting df is roughly matches that of additive trend filtering in the first row. Now we see that the first component fit is oversmoothed, yet the third is still undersmoothed.

Figure 2 displays the MSE curves from additive trend filtering, as a function of df. We see that trend filtering achieves a lower MSE, and moreover, its MSE curve is optimized at a lower df (i.e., less complex model) than that for smoothing splines. This is analogous to what is typically seen in the univariate setting (Tibshirani, 2014).

1.5 Multiple tuning parameters

In problems (1) and (12), consider generalizing the penalties in the criterions by

[TABLE]

respectively, which means we would now have $d$ tuning parameters $\lambda_{j}$, $j=1,\ldots,d$. In practice, selecting multiple tuning parameters is significantly more challenging than selecting a single one, as is needed in (1) and (12). However, as pointed out by a referee of this article, there has been a considerable amount of work dedicated to this problem by authors studying additive models built from splines (or other linear smoothers), e.g., Gu and Wahba (1991); Wood (2000); Fahrmeir and Lang (2001); Ruppert et al. (2003); Wood (2004); Kim and Gu (2004); Rue et al. (2009); Wood (2011); Wood et al. (2015, 2016). Many of these papers use an efficient computational approach based on restricted maximum likelihood (REML) for selecting $\lambda_{j}$, $j=1,\ldots,d$; see also Wood (2017) for a nice introduction and description of this approach. Unfortunately, as far as we see it, REML does not easily apply to additive trend filtering.

In this paper, we focus on a single tuning parameter $\lambda$ as in (1) and (12) mainly for simplicity; we are not suggesting that this is always the preferred formulation in practice. Many of our results (the basis formulation, and the uniqueness and degrees of freedom results in Section 2) carry over immediately to the multiple tuning parameter case. Others (the error bounds in Section 3) do not, though extending them may certainly be possible. Furthermore, the motivating example of the last subsection is meant to elucidate the differences in what additive smoothing splines and additive trend filtering can do with a single tuning parameter; a serious applied statistician, in just $d=3$ dimensions, would likely use REML or some related technique to fit a multiple tuning parameter smoothing spline model, which would bring it closer to additive trend filtering in performance here (it would be able to adjust to the variable smoothness across the components, though still not that within $f_{01}$). Of course, as $d$ grows larger, a separate tuning parameter per component will generally become infeasible, and the single tuning parameter comparisons will become more meaningful.

1.6 Summary of contributions

A summary of our contributions, and an outline for the rest of this paper, are given below.

• •

In Section 2, we investigate basic properties of the additive trend filtering model: an equivalent continuous-time formulation, a condition for uniqueness of component function estimates, and a simple formula for the effective degrees of freedom of the additive fit. We also introduce two estimators related to additive trend filtering, based on splines, that facilitate theoretical analysis (and are perhaps of interest in their own right).

• •

In Section 3, we derive error bounds for additive trend filtering. Assuming that the underlying regression function is additive, denoted by $f_{0}=\sum_{j=1}^{d}f_{0j}$, and that $\mathrm{TV}(f_{0j}^{(k)})$ is bounded, for $j=1,\ldots,d$, we prove that the $k$th order additive trend filtering estimator converges to $f_{0}$ at the rate $n^{-(2k+2)/(2k+3)}$ when the dimension $d$ is fixed (under weak assumptions), and at the rate $dn^{-(2k+2)/(2k+3)}$ when $d$ is growing (under stronger assumptions). We prove that these rates are optimal in a minimax sense, and also show that additive smoothing splines (or more generally, additive models built from linear smoothers of any kind) are suboptimal over such a class of functions $f_{0}$.

• •

In Section 4, we study the backfitting algorithm for additive trend filtering models, and give a connection between backfitting and an alternating projections scheme in the additive trend filtering dual problem. This inspires a new parallelized backfitting algorithm for additive trend filtering.

• •

In Section 5, we present empirical experiments and comparisons, and we also investigate the use of cross-validation for tuning parameter selection, and multiple tuning parameter models as in (13). In Section 6, we give a brief discussion.

2 Basic properties

In this section, we derive a number of basic properties of additive trend filtering estimates, starting with a representation for the estimates as continuous functions over $\mathbb{R}^{d}$ (rather than simply discrete fitted values at the input points).

2.1 Falling factorial representation

We may describe additive trend filtering in (12) as an estimation problem written in analysis form. The components are modeled directly by the parameters $\theta_{j}$, $j=1,\ldots,d$, and the desired structure is established by regularizing the discrete derivatives of these parameters, through the penalty terms $\|D^{(X_{j},k+1)}S_{j}\theta_{j}\|_{1}$, $j=1,\ldots,d$. Here, we present an alternative representation for (12) in basis form, where each component is expressed as a linear combination of basis functions, and regularization is applied to the coefficients in this expansion.

Before we derive the basis formulation that underlies additive trend filtering, we first recall the falling factorial basis (Tibshirani, 2014; Wang et al., 2014). Given knot points $t^{1}<\ldots<t^{n}\in\mathbb{R}$, the $k$th order falling factorial basis functions $h_{1},\ldots,h_{n}$ are defined by

[TABLE]

We denote $1\{t>a\}=1$ when $t>a$, and 0 otherwise. (Also, our convention is to define the empty product to be 1, so that $h_{1}(t)=1$.) The functions $h_{1},\ldots,h_{n}$ are piecewise polynomial functions of order $k$, and appear very similar in form to the $k$th order truncated power basis functions. In fact, when $k=0$ or $k=1$, the two bases are exactly equivalent (meaning that they have the same span). Similar to an expansion in the truncated power basis, an expansion in the falling factorial basis,

[TABLE]

is a continuous piecewise polynomial function, having a global polynomial structure determined by $\alpha^{1},\ldots,\alpha^{k+1}$, and exhibiting a knot—i.e., a change in its $k$th derivative—at the location $t^{i+k}$ when $\alpha^{i+k+1}\not=0$. But, unlike the truncated power functions, the falling factorial functions in (14) are not splines, and when $g$ (as defined above) has a knot at a particular location, it displays a change not only in its $k$th derivative at this location, but also in all lower order derivatives (i.e., all derivatives of orders $1,\ldots,k-1$).

Tibshirani (2014); Wang et al. (2014) establish a connection between univariate trend filtering and the falling factorial functions, and show that the trend filtering problem can be interpreted as a sparse basis regression problem using these functions. As we show next, the analogous result holds for additive trend filtering.

Lemma 1** (Falling factorial representation).**

For $j=1,\ldots,d$, let $h_{1}^{(X_{j})},\ldots,h_{n}^{(X_{j})}$ be the falling factorial basis in (14) with knots $(t^{1},\ldots,t^{n})=S_{j}X_{j}$, the $j$th dimension of the input points, properly sorted. Then the additive trend filtering problem (12) is equivalent to the problem

[TABLE]

in that, at any solutions in (12), (15), we have

[TABLE]

An alternative way of expressing problem (15) is

[TABLE]

where $\mathcal{H}_{j}=\mathrm{span}\{h^{(X_{j})}_{1},\ldots,h^{(X_{j})}_{n}\}$ is the span of the falling factorial basis over the $j$th dimension, and $f_{j}^{(k)}$ is the $k$th weak derivative of $f_{j}$, $j=1,\ldots,d$. In this form, at any solutions in (12), (16),

[TABLE]

Proof.

For $j=1,\ldots,d$, define the $k$th order falling factorial basis matrix $H^{(X_{j},k)}\in\mathbb{R}^{n\times n}$ by

[TABLE]

Note that the columns of $H^{(X_{j},k)}$ follow the order of the sorted inputs $S_{j}X_{j}$, but the rows do not; however, for $S_{j}H^{(X_{j},k)}$, both its rows and columns of follow the order of $S_{j}X_{j}$. From Wang et al. (2014), we know that

[TABLE]

for some matrix $C^{(X_{j},k+1)}\in\mathbb{R}^{(k+1)\times n}$, i.e.,

[TABLE]

Problem (15) is given by reparameterizing (12) according to $\theta_{j}=H^{(X_{j},k)}\alpha_{j}$, for $j=1,\ldots,d$. As for (16), the equivalence between this and (15) follows by noting that for $f_{j}=\sum_{\ell=1}^{n}\alpha^{\ell}_{j}h_{\ell}^{(X_{j})}$, we have

[TABLE]

and so $\mathrm{TV}(f_{j}^{(k)})=k!\sum_{\ell=k+2}^{n}|\alpha^{\ell}_{j}|$, for each $j=1,\ldots,d$. ∎

This lemma not only provides an interesting reformulation for additive trend filtering, it is also practically useful in that it allows us to perform interpolation or extrapolation using the additive trend filtering model. That is, from a solution $\hat{\theta}=(\hat{\theta}_{1},\ldots,\hat{\theta}_{d})$ in (12), we can extend each component fit $\hat{\theta}_{j}$ to the real line, by forming an appropriate linear combination of falling factorial functions:

[TABLE]

The coefficients above are determined by the relationship $\hat{\alpha}_{j}=(H^{(X_{j},k)})^{-1}\hat{\theta}_{j}$, and are easily computable given the highly structured form of $(H^{(X_{j},k)})^{-1}$, as revealed in (18). Writing the coefficients in block form, as in $\hat{\alpha}_{j}=(\hat{a}_{j},\hat{b}_{j})\in\mathbb{R}^{(k+1)}\times\mathbb{R}^{(n-k-1)}$, we have

[TABLE]

The first $k+1$ coefficients $\hat{a}_{j}$ index the pure polynomial functions $h^{(X_{j})}_{1},\ldots,h^{(X_{j})}_{k+1}$. These coefficients will be generically dense (the form of $C^{(X_{j},k+1)}$ is not important here, so we omit it for simplicity, but details are given in Appendix A.1). The last $n-k-1$ coefficients $\hat{b}_{j}$ index the knot-producing functions $h^{(X_{j})}_{k+2},\ldots,h^{(X_{j})}_{n}$, and when $(\hat{b}_{j})_{\ell}=\frac{1}{k!}(D^{(X_{j},k+1)}S_{j}\hat{\theta}_{j})_{\ell}\not=0$, the fitted function $\hat{f}_{j}$ exhibits a knot at the $(\ell+k)$th sorted input point among $S_{j}X_{j}$, i.e., at $X_{j}^{(\ell+k)}$. Figure 3 gives an example.

We note that the coefficients $\hat{\alpha}_{j}=(\hat{a}_{j},\hat{b}_{j})$ in (20), (21) can be computed in $O(n)$ operations and $O(1)$ memory. This makes extrapolation of the $j$th fitted function $\hat{f}_{j}$ in (19) highly efficient. Details are given in Appendix A.1.

2.2 Uniqueness of component fits

It is easy to see that, for the problem (12), the additive fit $\sum_{j=1}^{d}\hat{\theta}_{j}$ is always uniquely determined: denoting $\sum_{j=1}^{d}\theta_{j}=T\theta$ for a linear operator $T$ and $\theta=(\theta_{1},\ldots,\theta_{d})\in\mathbb{R}^{nd}$, the loss term $\|y-T\theta\|_{2}^{2}$ is strictly convex in the variable $T\theta$, and this, along with the convexity of the problem (12), implies a unique additive fit $T\hat{\theta}$, no matter the choice of solution $\hat{\theta}=(\hat{\theta}_{1},\ldots,\hat{\theta}_{d})\in\mathbb{R}^{nd}$.

On the other hand, when $d>1$, the criterion in (12) is not strictly convex in $\theta$, and hence there need not be a unique solution $\hat{\theta}$, i.e., the individual components fits $\hat{\theta}_{j}$, $j=1,\ldots,d$ need not be uniquely determined. We show next that uniqueness of the component fits can be guaranteed under some conditions on the input matrix $X=[X_{1}\,\cdots\,X_{d}]\in\mathbb{R}^{n\times d}$. We will rely on the falling factorial representation for additive trend filtering, introduced in the previous subsection, and on the notion of general position: a matrix $A\in\mathbb{R}^{m\times p}$ is said to have columns in general position provided that, for any $\ell<\min\{m,p\}$, subset of $\ell+1$ columns denoted $A_{i_{1}},\ldots,A_{i_{\ell+1}}$, and signs $s_{1},\ldots,s_{\ell+1}\in\{-1,1\}$, the affine span of $\{s_{1}A_{i_{1}},\ldots,s_{\ell+1}A_{i_{\ell+1}}\}$ does not contain any element of $\{\pm A_{i}:i\not=i_{1},\ldots,i_{\ell+1}\}$. Informally, if the columns of $A$ are not in general position, then there must be some small subset of columns that are affinely dependent.

Lemma 2** (Uniqueness).**

For $j=1,\ldots,d$, let $H^{(X_{j},k)}\in\mathbb{R}^{n\times n}$ be the falling factorial basis matrix constructed over the sorted $j$th dimension of inputs $S_{j}X_{j}\in\mathbb{R}^{n}$, as in (17). Decompose $H^{(X_{j},k)}$ into its first $k+1$ columns $P^{(X_{j},k)}\in\mathbb{R}^{n\times(k+1)}$, and its last $n-k-1$ columns $K^{(X_{j},k)}\in\mathbb{R}^{n\times(n-k-1)}$. The former contains evaluations of the pure polynomials $h^{(X_{j})}_{1},\ldots,h^{(X_{j})}_{k+1}$; the latter contains evaluations of the knot-producing functions $h^{(X_{j})}_{k+2},\ldots,h^{(X_{j})}_{n}$. Also, let $\tilde{P}^{(X_{j},k)}$ denote the matrix $P^{(X_{j},k)}$ with its first column removed, for $j=1,\ldots,d$, and $M=I-\mathbbm{1}\mathbbm{1}^{T}/n$. Define

[TABLE]

the product of $M$ and the columnwise concatenation of $\tilde{P}^{(X_{j},k)}$, $j=1,\ldots,d$. Let $UU^{T}$ denote the projection operator onto the space orthogonal to the column span of $\tilde{P}$, where $U\in\mathbb{R}^{n\times(n-kd-1)}$ has orthonormal columns, and define

[TABLE]

the product of $U^{T}M$ and the columnwise concatenation of $K^{(X_{j},k)}$, $j=1,\ldots,d$. A sufficient condition for uniqueness of the additive trend filtering solution in (12) can now be given in two parts.

If $\tilde{K}$ has columns in general position, then the knot-producing parts of all component fits are uniquely determined, i.e., for each $j=1,\ldots,d$, the projection of $\hat{\theta}_{j}$ onto the column space of $K^{(X_{j},k)}$ is unique. 2. 2.

If in addition $\tilde{P}$ has full column rank, then the polynomial parts of component fits are uniquely determined, i.e., for each $j=1,\ldots,d$, the projection of $\hat{\theta}_{j}$ onto the column space of $P^{(X_{j},k)}$ is unique, and thus the component fits $\hat{\theta}_{j}$, $j=1,\ldots,d$ are all unique.

The proof is deferred to Appendix A.2. To rephrase, the above lemma decomposes each component of the additive trend filtering solution according to

[TABLE]

where $\hat{\theta}_{j}^{\mathrm{poly}}$ exhibits a purely polynomial trend over $S_{j}X_{j}$, and $\hat{\theta}_{j}^{\mathrm{knot}}$ exhibits a piecewise polynomial trend over $S_{j}X_{j}$, and hence determines the knot locations, for $j=1,\ldots,d$. The lemma shows that the knot-producing parts $\hat{\theta}_{j}^{\mathrm{knot}}$, $j=1,\ldots,d$ are uniquely determined when the columns of $\tilde{K}$ are in general position, and the polynomial parts $\hat{\theta}_{j}^{\mathrm{knot}}$, $j=1,\ldots,d$ are unique when the columns of $\tilde{K}$ are in general position, and the columns of $\tilde{P}$ are linearly independent.

The conditions placed on $\tilde{P},\tilde{K}$ in Lemma 2 are not strong. When $n>kd$, and the elements of input matrix $X$ are drawn from a density over $\mathbb{R}^{nd}$, it is not hard to show that $\tilde{P}$ has full column rank with probability 1. We conjecture that, under the same conditions, $\tilde{K}$ will also have columns in general position with probability 1, but do not pursue a proof.

Remark 1** (Relationship to concurvity).**

It is interesting to draw a connection to Buja et al. (1989). In the language used by these authors, when $\tilde{P}$ has linearly dependent columns, we say that the predictor variables display concurvity, i.e., linear dependence after nonlinear (here, polynomial) transformations are applied. Buja et al. (1989) establish that the components in the additive model (1), built with quadratic penalties, are unique provided there is no concurvity between variables. In comparison, Lemma 2 establishes uniqueness of the additive trend filtering components when there is no concurvity between variables, and additionally, the columns of $\tilde{K}$ are in general position. The latter two conditions together can be seen as requiring no generalized concurvity—if $\tilde{K}$ were to fail the general position assumption, then there would be a small subset of the variables that are linearly dependent after nonlinear (here, piecewise polynomial) transformations are applied.

2.3 Dual problem

Let us abbreviate $D_{j}=D^{(X_{j},k+1)}$, $j=1,\ldots,d$ for the penalty matrices in the additive trend filtering problem (12). Basic arguments in convex analysis, deferred to Appendix A.3, show that the dual of problem (12) can be expressed as:

[TABLE]

and that primal and dual solutions in (12), (24) are related by:

[TABLE]

From the form of (24), it is clear that we can write the (unique) dual solution as $\hat{u}=\Pi_{U}(Y-\bar{Y}\mathbbm{1})$, where $\Pi_{U}$ is the (Euclidean) projection operator onto $U$. Moreover, using (25), we can express the additive fit as $\sum_{j=1}^{d}\hat{\theta}_{j}=(\mathrm{Id}-\Pi_{U})(Y-\bar{Y}\mathbbm{1})$, where $\mathrm{Id}-\Pi_{U}$ is the operator that gives the residual from projecting onto $U$. These relationships will be revisited in Section 4, where we return to the dual perspective, and argue that the backfitting algorithm for the additive trend filtering problem (12) can be seen as a type of alternating projections algorithm for its dual problem (24).

2.4 Degrees of freedom

In general, given data $Y\in\mathbb{R}^{n}$ with $\mathbb{E}(Y)=\eta$, $\mathrm{Cov}(Y)=\sigma^{2}I$, and an estimator $\hat{\eta}$ of $\eta$, recall that we define the effective degrees of freedom of $\hat{\eta}$ as (Efron, 1986; Hastie and Tibshirani, 1990):

[TABLE]

where $\hat{\eta}(Y)=(\hat{\eta}^{1}(y),\ldots,\hat{\eta}^{n}(Y))$. Roughly speaking, the above definition sums the influence of the $i$th component $Y^{i}$ on its corresponding fitted value $\hat{\eta}^{i}(Y)$, across $i=1,\ldots,n$. A precise understanding of degrees of freedom is useful for model comparisons (recall the x-axis in Figure 2), and other reasons. For linear smoothers, in which $\hat{\eta}(Y)=SY$ for some $S\in\mathbb{R}^{n\times n}$, it is clear that $\mathrm{df}(\hat{\eta})=\mathrm{tr}(S)$, the trace of $S$. (This also covers additive models whose components are built from univariate linear smoothers, because in total these are still just linear smoothers: the additive fit is still just a linear function of $Y$.)

Of course, additive trend filtering is a not a linear smoother; however, it is a particular type of generalized lasso estimator, and degrees of freedom for such a class of estimators is well-understood (Tibshirani and Taylor, 2011, 2012). The next result is an consequence of existing generalized lasso theory, proved in Appendix A.4.

Lemma 3** (Degrees of freedom).**

Assume the conditions of Lemma 2, i.e., that the matrix $\tilde{P}$ in (22) has full column rank, and the matrix $\tilde{K}$ in (23) is in general position. Assume also that the response is Gaussian, $Y\sim N(\eta,\sigma^{2}I)$, and treat the input points $X^{i}\in\mathbb{R}^{d}$, $i=1,\ldots,n$ as fixed and arbitrary, as well as the tuning parameter value $\lambda\geq 0$. Then the additive trend filtering fit from (12) has degrees of freedom

[TABLE]

Remark 2** (The effect of shrinkage).**

Lemma 3 says that for an unbiased estimate of the degrees of freedom of the additive trend filtering fit, we count the number of knots in each component fit $\hat{\theta}_{j}$ (recall that this is the number of nonzeros in $D^{(X_{j},k+1)}\hat{\theta}_{j}$, i.e., the number of changes in the discrete $(k+1)$st derivative), add them up over $j=1,\ldots,d$, and add $kd$. This may seem surprising, as these knot locations are chosen adaptively based on the data $Y$. But, such adaptivity is counterbalanced by the shrinkage induced by the $\ell_{1}$ penalty in (12) (i.e., for each component fit $\hat{\theta}_{j}$, there is shrinkage in the differences between the attained $k$th derivatives on either side of a selected knot). See Tibshirani (2015) for a study of this phenomenon.

Remark 3** (Easy unbiased degrees of freedom estimation).**

It is worth emphasizing that an unbiased estimate from Lemma 3 for the degrees of freedom of the total fit in additive trend filtering is very easy to calculate: we scan the individual component fits and add up the number of knots that appear in each one. The same cannot be said for additive smoothing splines, or additive models built from univariate linear smoothers, in general. Although computing the fit itself is typically cheaper with additive linear smoothers than with additive trend filtering, computing the degrees of freedom is more challenging. For example, for the additive model in (1) built with quadratic penalties, we have

[TABLE]

where $F\in\mathbb{R}^{n\times nd}$ has $d$ copies of the centering matrix $M=I-\mathbbm{1}\mathbbm{1}^{T}/n\in\mathbb{R}^{n\times n}$ stacked across its columns, $Q\in\mathbb{R}^{nd\times nd}$ is a block diagonal matrix with blocks $MQ_{j}M$, $j=1,\ldots,d$, and $A^{+}$ denotes the Moore-Penrose pseudoinverse of a matrix $A$. The above formula does not obviously decompose into a sum of quantities across components, and is nontrivial to compute post optimization of (1), specifically when a backfitting algorithm as in (2) has been used to compute a solution.

2.5 Two related additive spline estimators

From its equivalent formulation in (16), additive trend filtering is seen to be closely related to two other additive spline estimators, which we introduce here. Consider, for univariate function classes $\mathcal{S}_{j}$, $j=1,\ldots,d$, the problem

[TABLE]

When each $\mathcal{S}_{j}$, $j=1,\ldots,d$ is the set of $k$ times weakly differentiable functions, we call the solution in (26) the additive locally adaptive regression spline of order $k\geq 0$, as it is the natural extension of the univariate estimator considered in Mammen and van de Geer (1997). Denote by $\hat{f}_{j}$, $j=1,\ldots,d$ this solution; the representation arguments used by these authors apply immediately to the additive setting, and imply that each $\hat{f}_{j}$, $j=1,\ldots,d$ is indeed a spline of degree $k$ (justifying the choice of name). The same arguments show that, for $k=0$ or $k=1$, the knots of the spline $\hat{f}_{j}$ lie among the $j$th dimension of the input points $X^{1}_{j},\ldots,X^{n}_{j}$, for $j=1,\ldots,d$, but for $k\geq 2$, this need not be true, and in general the components will be splines with knots at locations other than the inputs.

We can facilitate computation by taking $\mathcal{S}_{j}=\mathcal{G}_{j}$, where $\mathcal{G}_{j}$ is the set of splines of degree $k$ with knots lying among the $j$th dimension of inputs $X^{1}_{j},\ldots,X^{n}_{j}$, for $j=1,\ldots,d$. We call the resulting solution the restricted additive locally adaptive regression spline of order $k\geq 0$. More precisely, we require that the splines in $\mathcal{G}_{j}$ have knots in a set $T_{j}$, which, writing $t_{j}=S_{j}X_{j}$ for the sorted inputs along the $j$th dimension, is defined by

[TABLE]

i.e., defined by removing $k+1$ input points at the boundaries, for $j=1,\ldots,d$. Setting $\mathcal{S}_{j}=\mathcal{G}_{j}$, $j=1,\ldots,d$ makes (26) a finite-dimensional problem, just like (16). When $k=0$ or $k=1$, as claimed in Section 2.1 (and shown in Tibshirani (2014)), the falling factorial functions are simply splines, which means that $\mathcal{H}_{j}=\mathcal{G}_{j}$ for $j=1,\ldots,d$, hence additive trend filtering and restricted additive locally adaptive regression splines are the same estimator. When $k\geq 2$, this is no longer true, and they are not the same. Additive trend filtering will be much easier to compute, since $\mathrm{TV}(g^{(k)})$ does not admit a nice representation in terms of discrete derivatives for a $k$th order spline (and yet it does for a $k$th order falling factorial function, as seen in (12)).

To summarize, additive locally adaptive splines, restricted additive locally adaptive splines, and additive trend filtering all solve a problem of the form (26) for different choices of function classes $\mathcal{S}_{j}$, $j=1,\ldots,d$. For $k=0$ or $k=1$, these three estimators are equivalent. For $k\geq 2$, they will be generically different, though our intuition tells us that their differences should not be too large: the unrestricted problem admits a solution that is a spline in each component; the restricted problem simply forces these splines to have knots at the input points; and the trend filtering problem swaps splines for falling factorial functions, which are highly similar in form. Next, we give theory that confirms this intuition, in large samples.

3 Error bounds

We derive error bounds for additive trend filtering and additive locally adaptive regression splines (both the unrestricted and restricted variants), when the underlying regression function is additive, and has components whose derivatives are of bounded variation. These results are actually special cases of a more general result we prove in this section, on a generic roughness-regularized additive estimator, where we assume a certain decay for the entropy of the unit ball in the roughness operator. We treat separately the settings in which the dimension $d$ of the input space is fixed and growing. We also complement our error rates with minimax lower bounds. We start by introducing helpful notation.

3.1 Notation

Given a distribution $Q$ supported on a set $D$, and i.i.d. samples $X^{i}$, $i=1,\ldots,n$ from $Q$, denote by $Q_{n}$ the associated empirical distribution. We define the $L_{2}(Q)$ and $L_{2}(Q_{n})$ inner products, denoted $\langle\cdot,\cdot\rangle_{L_{2}(Q)}$ and $\langle\cdot,\cdot\rangle_{L_{2}(Q_{n})}$, respectively, over functions $m,r:D\to\mathbb{R}$

[TABLE]

Definitions for the corresponding $L_{2}(Q)$ and $L_{2}(Q_{n})$ norms, denoted $\|\cdot\|_{L_{2}(Q)}$ and $\|\cdot\|_{L_{2}(Q_{n})}$, respectively, arise naturally from these inner products, defined by

[TABLE]

Henceforth, we will abbreviate subscripts when using these norms and inner products, writing $\|\cdot\|_{2}$ and $\|\cdot\|_{n}$ for the $L_{2}(Q)$ and $L_{2}(Q_{n})$ norms, respectively, and similarly for the inner products. This abbreviated notation omits the underlying distribution $Q$; thus, unless explicitly stated otherwise, the underlying distribution should always be interpreted as the distribution of the input points. We will often call $\|\cdot\|_{2}$ the $L_{2}$ norm and $\|\cdot\|_{n}$ the empirical norm, and similarly for inner products.

In what follows, of particular interest will be the case when $D=[0,1]^{d}$, and $m:[0,1]^{d}\to\mathbb{R}$ is an additive function, of the form

[TABLE]

which we write to mean $m(x)=\sum_{j=1}^{d}m_{j}(x_{j})$. In a slight abuse of notation (overload of notation), for each $j=1,\ldots,d$, we will abbreviate the $L_{2}(Q_{j})$ norm by $\|\cdot\|_{2}$, where $Q_{j}$ is the $j$th marginal of $Q$, and will also abbreviate $L_{2}(Q_{jn})$ norm by $\|\cdot\|_{n}$, where $Q_{jn}$ is the empirical distribution of $X^{i}_{j}$, $i=1,\ldots,n$. We will use similar abbreviations for the inner products.

A few more general definitions are in order. We denote the $L_{\infty}$ norm, also called the sup norm, of a function $f:D\to\mathbb{R}$ by $\|f\|_{\infty}=\mathop{\mathrm{ess\,\,sup}}_{z\in D}|f(z)|$. For a functional $\nu$, acting on functions from $D$ to $\mathbb{R}$, we write $B_{\nu}(\delta)$ for the $\nu$-ball of radius $\delta>0$, i.e., $B_{\nu}(\delta)=\{f:\nu(f)\leq\delta\}$. We abbreviate $B_{n}(\delta)$ for the $\|\cdot\|_{n}$-ball of radius $\delta$, $B_{2}(\delta)$ for the $\|\cdot\|_{2}$-ball of radius $\delta$, and $B_{\infty}(\delta)$ for the $\|\cdot\|_{\infty}$-ball of radius $\delta$. We will use these concepts fluidly, without explicit reference to the domain $D$ (or its dimensionality), as the meaning should be clear from the context.

Lastly, for a set $S$ and norm $\|\cdot\|$, we define the covering number $N(\delta,\|\cdot\|,S)$ to be the smallest number of $\|\cdot\|$-balls of radius $\delta$ to cover $S$, and the packing number $M(\delta,\|\cdot\|,S)$ to be the largest number of disjoint $\|\cdot\|$-balls of radius $\delta$ that are contained in $S$. We call $\log N(\delta,\|\cdot\|,S)$ the entropy number.

3.2 Error bounds for a fixed dimension $d$

We consider error bounds for the generic roughness-penalized estimator defined as a solution of

[TABLE]

where $\mathcal{S}_{j}$, $j=1,\ldots,d$ are univariate function spaces, and $J$ is a regularizer that acts on univariate functions. We assume in this subsection that the dimension $d$ of the input space is fixed, i.e., it does not grow with $n$. Before stating our main result in this setting, we list and briefly discuss our other assumptions. First, we give our assumptions on the data generation process.

Assumption A1 (Continuous input distribution)****.

The input points $X^{i}$, $i=1,\ldots,n$ are i.i.d. from a continuous distribution $Q$ supported on $[0,1]^{d}$.

**Assumption B1 (Generic regression model, with sub-Gaussian

errors)****.**

The responses $Y^{i}$, $i=1,\ldots,n$ follow the model

[TABLE]

with overall mean $\mu\in\mathbb{R}$, where $\sum_{i=1}^{n}f_{0}(X^{i})=0$ for identifiability. The errors $\epsilon^{i}$, $i=1,\ldots,n$ are uniformly sub-Gaussian and have mean zero, i.e.,

[TABLE]

for a constant $\sigma>0$. The errors and input points are independent.

Assumptions Assumption A1 (Continuous input distribution) and Assumption B1 (Generic regression model, with sub-Gaussian errors) are very weak. We do not place any specific smoothness or additivity conditions on the underlying regression function $f_{0}$, as our error bound in Theorem 1 will involve the error of the closest additive function to $f_{0}$, whose components lie in the given function spaces $\mathcal{S}_{j}$, $j=1,\ldots,d$.

Next, we present our assumptions on the regularizer $J$. We write $\|\cdot\|_{Z_{n}}$ for the empirical norm defined over a set of univariate points $Z_{n}=\{z^{1},\ldots,z^{n}\}\subseteq[0,1]$, i.e., $\|g\|_{Z_{n}}^{2}=\frac{1}{n}\sum_{i=1}^{n}g^{2}(z^{i})$.

**Assumption C1 (Seminorm regularizer, null space of

polynomials)****.**

The regularizer $J$ is a seminorm, and its domain is contained in the space of $k$ times weakly differentiable functions, for an integer $k\geq 0$. Furthermore, its null space contains all $k$th order polynomials, i.e., $J(g)=0$ for all $g(t)=t^{\ell}$, $\ell=0,\ldots,k$.

Assumption C2 (Relative boundedness of derivatives)****.

There is a constant $L>0$ such that $\mathop{\mathrm{ess\,\,sup}}_{t\in[0,1]}g^{(k)}(t)-\mathop{\mathrm{ess\,\,inf}}_{t\in[0,1]}g^{(k)}(t)\leq L$ for $g\in B_{J}(1)$ (with $g^{(k)}$ the $k$th weak derivative of $g$).

Assumption C3 (Entropy bound)****.

There are constants $0<w<2$ and $K>0$ such that

[TABLE]

Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C3 (Entropy bound) on the regularizer $J$ are not strong, and are satisfied by various common regularizers, e.g., $J(g)=[\int_{0}^{1}(g^{(k+1)}(t))^{2}\,dt]^{1/2}$ or $J(g)=\mathrm{TV}(g^{(k)})$ (the latter studied shortly).

We now state our main result in the fixed $d$ case, which is proved in Appendix A.5, A.6.

Theorem 1**.**

Assume Assumption A1 (Continuous input distribution), Assumption B1 (Generic regression model, with sub-Gaussian errors) on the data distribution, and assume Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C3 (Entropy bound) on the seminorm $J$. Also, assume that the dimension $d$ of the input space is fixed. Let $C_{n}\geq 1$ be an arbitrary sequence. There exist constants $c_{1},c_{2},c_{3},n_{0}>0$, that depend only on $d,\sigma,k,L,K,w$, such that for all $c\geq c_{1}$, $n\geq n_{0}$, and tuning parameter values $\lambda\geq cn^{w/(2+w)}C_{n}^{-(2-w)/(2+w)}$, any solution in (28) satisfies

[TABLE]

with probability at least $1-\exp(-c_{2}c)-\exp(-c_{3}\sqrt{n})$, simultaneously over all $\tilde{f}=\sum_{j=1}^{d}\tilde{f}_{j}$, feasible for the problem (28), such that $\|\tilde{f}-f_{0}\|_{n}\leq\max\{C_{n},\sum_{j=1}^{d}J(\tilde{f}_{j})\}$.

Remark 4** (Error bound for additive, $J$-smooth $f_{0}$).**

Assume $f_{0}=\sum_{j=1}^{d}f_{0j}$, where $f_{0j}\in\mathcal{S}_{j}$, $j=1,\ldots,d$, and $\sum_{j=1}^{d}J(f_{0j})\leq C_{n}$. Letting $\tilde{f}=f_{0}$, the approximation error term in (29) (the first term on the right-hand side) is zero, and for $\lambda=cn^{w/(2+w)}C_{n}^{-(2-w)/(2+w)}$, the result in the theorem reads

[TABLE]

with probability at least $1-\exp(-c_{2}c)-\exp(-c_{3}\sqrt{n})$. As we will see in the minimax lower bound in Theorem 3 (plugging in $c_{n}=C_{n}/d$, and taking $d$ to be a constant), the rate $n^{-2/(2+w)}C_{n}^{2w/(2+w)}$ is optimal for such a class of functions.

Remark 5** **(**Distance to best additive, $J$-smooth approximation of

$f_{0}$).**

The arguments used to establish the oracle-type inequality (29) also imply a result on the empirical norm error between $\hat{f}$ and the best additive approximation of $f_{0}$. To be precise, let $(f^{\mathrm{best}}_{1},\ldots,f^{\mathrm{best}}_{d})$ denote a solution in the population-level problem

[TABLE]

We note that (31) has “half” of the regularization of problem (28), as it uses a penalty parameter of $\lambda/2$ versus $\lambda$. We can think of $f^{\mathrm{best}}=\sum_{j=1}^{d}f^{\mathrm{best}}_{j}$ as the best additive, $J$-smooth approximation of $f_{0}$, where $\lambda$ as usual controls the level of smoothness. The following is a consequence of the proof of Theorem 1, verified in Appendix A.7: assume that $\|f^{\mathrm{best}}-f_{0}\|_{n}\leq\max\{C_{n},\sum_{j=1}^{d}J(f^{\mathrm{best}}_{j})\}$ almost surely (with respect to $Q$), for sufficiently large $\lambda$; then any solution in (28) satisfies for all $c\geq c_{1}$, $n\geq n_{0}$, and $\lambda\geq cn^{w/(2+w)}C_{n}^{-(2-w)/(2+w)}$,

[TABLE]

with probability at least $1-\exp(-c_{2}c)-\exp(-c_{3}\sqrt{n})$, where as before $c_{1},c_{2},c_{3},n_{0}>0$ are constants that depend only on $d,\sigma,k,L,K,w$. Notably, the right-hand side in the bound (32) does not depend on the approximation error; in particular, we do not even require $\|f^{\mathrm{best}}-f_{0}\|_{n}$ to converge to zero. This is analogous to classical results from Stone (1985).

We examine a special case of the generic problem (28) when the regularizer is $J(g)=\mathrm{TV}(g^{(k)})$, and derive implications of the above Theorem 1 for additive locally regression adaptive splines and additive trend filtering, corresponding to different choices of the function classes $\mathcal{S}_{j}$, $j=1,\ldots,d$ in (28). We must introduce an additional (weak) assumption on the input distribution, for the results on restricted locally adaptive regression splines and trend filtering.

Assumption A2 (Bounded input density)****.

The density of the input distribution $Q$ is bounded below by a constant $b_{0}>0$.

Here is our result for additive locally adaptive splines and additive trend filtering. The proof is given in Appendix A.8, A.9.

Corollary 1**.**

Assume Assumption A1 (Continuous input distribution), Assumption B1 (Generic regression model, with sub-Gaussian errors) on the data distribution. Also, assume that the dimension $d$ of the input space is fixed, and that the underlying regression function is additive, $f_{0}=\sum_{j=1}^{d}f_{0j}$, where the components $f_{0j}$, $j=1,\ldots,d$ are $k$ times weakly differentiable, such that $\sum_{j=1}^{d}\mathrm{TV}(f_{0j}^{(k)})\leq C_{n}$ for a sequence $C_{n}\geq 1$. For $J(g)=\mathrm{TV}(g^{(k)})$, Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C3 (Entropy bound) hold with $L=1$ and $w=1/(k+1)$. Furthermore, the following is true of the estimator defined by problem (28).

• (a)

Let $\mathcal{S}_{j}$ be the set of all $k$ times weakly differentiable functions, for each $j=1,\ldots,d$. There are constants $c_{1},c_{2},c_{3},n_{0}>0$, depending only on $d,\sigma,k$, such that for all $c\geq c_{1}$ and $n\geq n_{0}$, any solution in the additive locally adaptive regression spline problem (28), with tuning parameter value $\lambda=cn^{1/(2k+3)}C_{n}^{-(2k+1)/(2k+3)}$, satisfies

[TABLE]

with probability at least $1-\exp(-c_{2}c)-\exp(-c_{3}\sqrt{n})$.

• (b)

*Let $\mathcal{S}_{j}=\mathcal{G}_{j}$, the set of *$k$th degree splines with knots in the set $T_{j}$ in (27), for $j=1,\ldots,d$, and assume Assumption A2 (Bounded input density) on the input density. Then there are constants $c_{1},c_{2},c_{3},n_{0}>0$, that depend only on $d,b_{0},\sigma,k$, such that for all $c\geq c_{1}$ and $n(\log n)^{-(1+1/k)}\geq n_{0}C_{n}^{(2k+2)/(2k^{2}+2k-1)}$, any solution in the restricted additive locally adaptive spline problem (28), with $\lambda=cn^{1/(2k+3)}C_{n}^{-(2k+1)/(2k+3)}$, satisfies the same result in (33), with probability at least $1-\exp(-c_{2}c)-c_{3}/n$.

• (c)

*Let $\mathcal{S}_{j}=\mathcal{H}_{j}$, the set of *$k$*th degree falling factorial functions defined over $X_{j}$ (the *$j$th dimension of inputs), for $j=1,\ldots,d$, and assume Assumption A2 (Bounded input density). Then there exist constants $c_{1},c_{2},c_{3},n_{0}>0$, that depend only on $d,b_{0},\sigma,k$, such that for all $c\geq c_{1}$ and $n(\log n)^{-(2k+3)}\geq n_{0}C_{n}^{4k+4}$, any solution in the additive trend filtering problem (28), with $\lambda=cn^{1/(2k+3)}C_{n}^{-(2k+1)/(2k+3)}$, satisfies (33), with probability at least $1-\exp(-c_{2}c)-c_{3}/n$.

Remark 6** (Spline and falling factorial approximants).**

For part (a) of the corollary, the approximation error (the first term on the right-hand side) in (30) is zero by definition, and we need only verify Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C3 (Entropy bound) for the regularizer $J(g)=\mathrm{TV}(g^{(k)})$. Parts (b) and (c) require control over the approximation error, because the underlying regression function $f_{0}=\sum_{j=1}^{d}f_{0j}$ need not have components that lie in the chosen function spaces $\mathcal{S}_{j}$, $j=1,\ldots,d$. To be clear: for $k=0$ or $k=1$, as discussed in Section 2.5, all three problems considered in parts (a), (b), (c) are equivalent; hence parts (b) and (c) really only concern the case $k\geq 2$. For both of these parts, we control the approximation error by controlling the univariate approximation error and then applying the triangle inequality. For part (b), we use a special spline quasi-interpolant from Proposition 7 in Mammen and van de Geer (1997) (who in turn construct this using results from de Boor (1978)); for part (c), we develop a new falling factorial approximant that may be of independent interest.

3.3 Error bounds for a growing dimension $d$

In this subsection, we allow the input dimension $d$ to grow with the sample size $n$. To keep our analysis as clean as possible, we consider a constrained version of the problem (28), namely

[TABLE]

for a tuning parameter $\delta>0$. (The penalized problem (28) can also be analyzed in the setting of growing $d$, but we find that the analysis is messier and requires more assumptions in order to obtain the same results.) Instead of Assumption A1 (Continuous input distribution), we now use the following assumption in the input distribution.

Assumption A3 (Product input distribution)****.

The input points $X^{i}$, $i=1,\ldots,n$ are i.i.d. from a continuous distribution $Q$ supported on $[0,1]^{d}$, that decomposes as $Q=Q_{1}\times\cdots\times Q_{d}$, where the density of each $Q_{j}$ is lower and upper bounded by constants $b_{1},b_{2}>0$, for $j=1,\ldots,d$.

Assumption Assumption A3 (Product input distribution) is fairly restrictive, since it requires the input distribution $Q$ to be independent across dimensions of the input space. The reason we use this assumption: when $Q=Q_{1}\times\cdots\times Q_{d}$, additive functions enjoy a key decomposability property in terms of the (squared) $L_{2}$ norm defined with respect to $Q$. In particular, if $m=\sum_{j=1}^{d}m_{j}$ has components with $L_{2}$ mean zero, denoted by $\bar{m}_{j}=\int_{0}^{1}m_{j}(x_{j})\,dQ_{j}(x_{j})=0$, $j=1,\ldots,d$, then we have

[TABLE]

This is explained by the fact that each pair of components $m_{j}$, $m_{\ell}$ with $j\not=\ell$ are orthogonal with respect to the $L_{2}$ inner product, since

[TABLE]

The above orthogonality, and thus the decomposability property in (35), is only true because of the product form $Q=Q_{1}\times\cdots\times Q_{d}$. Such decomposability is not generally possible with the empirical norm (the inner products between components do not vanish even if all empirical means are zero). In the proof of Theorem 2, we move from considering the empirical norm of the error vector to the $L_{2}$ norm, in order to leverage the property in (35), which eventually leads to an error rate that has a linear dependence on the dimension $d$. In the absence of $L_{2}$ decomposability, the same error rate can be achieved with a weaker incoherence bound, as in (40); see Remark 9 after the theorem.

We now state our main result in the growing $d$ case, whose proof is in Appendix A.10, A.11.

Theorem 2**.**

Assume Assumption A3 (Product input distribution), Assumption B1 (Generic regression model, with sub-Gaussian errors) on the data distribution, and assume Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C3 (Entropy bound) on the seminorm $J$. Let $\delta\geq 1$ be arbitrary. There are constants $c_{1},c_{2},c_{3},n_{0}>0$, that depend only on $b_{1},b_{2},\sigma,k,L,K,w$, such that for all $c\geq c_{1}$ and $n\geq n_{0}(d\delta)^{1+w/2}$, any solution in (34) satisfies both

[TABLE]

with probability at least $1-\exp(-c_{2}c)-c_{3}/n$, simultaneously over all functions $\tilde{f}=\sum_{j=1}^{d}\tilde{f}_{j}$, feasible for the problem (34).

Remark 7** (Error bound for additive, $J$-smooth $f_{0}$).**

Assume $f_{0}=\sum_{j=1}^{d}f_{0j}$, where $f_{0j}\in\mathcal{S}_{j}$ and $J(f_{0j})\leq c_{n}$, $j=1,\ldots,d$, for a sequence $c_{n}\geq 1$. Letting $\tilde{f}=f_{0}$, and $\delta=c_{n}$, the results in (36), (37) translate to

[TABLE]

*with probability at least $1-\exp(-c_{2}c)-c_{3}/n$, provided that $n\geq n_{0}(dc_{n})^{1+w/2}$. From the minimax lower bound in Theorem 3, we can see that the optimal rate for such a class of functions is in fact $dn^{-2/(2+w)}c_{n}^{2w/(2+w)}$, which reveals that the rates in (38) are tight when $c_{n}$ is a constant, but not when $c_{n}$ grows with $n$. It is worth noting that the dependence of the bounds on $c_{n}$ in Theorem 2 (and hence in (38)) can be improved to have the optimal scaling of $c_{n}^{2w/(2+w)}$ by assuming that $f_{0}$ is sup norm bounded, and additionally placing a sup norm bound on the components in (34). This feels like an unnecessary restriction, so we prefer to present results without it, as in Theorem 2 (and (38)). *

Remark 8** **(**Distance to best additive, $J$-smooth approximation of

$f_{0}$).**

A consequence of the proof of (36) is a bound on the empirical norm error between $\hat{f}$ and the best additive approximation of $f_{0}$. To be precise, let $f^{\mathrm{best}}=\sum_{j=1}^{d}f^{\mathrm{best}}_{j}$ minimize $\|\sum_{j=1}^{d}\tilde{f}_{j}-f_{0}\|_{n}^{2}$ over all additive functions $\tilde{f}=\sum_{j=1}^{d}\tilde{f}_{j}$ feasible for problem (34). Then following directly from (85) in the proof of Theorem 2, we have for all $c\geq c_{1}$ and $n\geq n_{0}(d\delta)^{1+w/2}$,

[TABLE]

with probability at least $1-\exp(-c_{2}c)-c_{3}/n$, where again $c_{1},c_{2},c_{3},n_{0}>0$ are constants that depend on $b_{1},b_{2},\sigma,k,L,K,w$. Just as we saw in fixed $d$ case, the right-hand side in (39) does not depend on the approximation error $\|f^{\mathrm{best}}-f_{0}\|_{n}$, which is analogous to classical results from Stone (1985).

Remark 9** ($L_{2}$ decomposability and incoherence**).

The decomposability property in (35) is critical in obtaining the sharp (linear) dependence on $d$ in the error rates (36), (37). However, it is worth noting that all that is needed in the proof is in fact a lower bound of the form

[TABLE]

for a constant $\phi_{0}>0$, rather than an exact equality, as in (35). The condition (40) is an incoherence condition that can hold for nonproduct distributions $Q$, over an appropriate class of functions (additive functions with smooth components), provided that the correlations between components of $Q$ are not too large. See Meier et al. (2009); van de Geer (2014) for similar incoherence conditions.

Next we present our results for additive locally adaptive regression splines (both unresricted and restricted variants) and additive trend filtering. The proof is in Appendix A.12.

Corollary 2**.**

Assume Assumption A3 (Product input distribution), Assumption B1 (Generic regression model, with sub-Gaussian errors) on the data distribution. Also, assume that the underlying regression function is additive, $f_{0}=\sum_{j=1}^{d}f_{0j}$, where the components $f_{0j}$, $j=1,\ldots,d$ are $k$ times weakly differentiable, such that $\mathrm{TV}(f_{0j}^{(k)})\leq c_{n}$, $j=1,\ldots,d$, for a sequence $c_{n}\geq 1$. Then for $J(g)=\mathrm{TV}(g^{(k)})$, the following is true of the estimator defined by problem (34).

• (a)

Let $\mathcal{S}_{j}$ be the space of all $k$ times weakly differentiable functions, for each $j=1,\ldots,d$. There exist constants $c_{1},c_{2},c_{3},n_{0}>0$, that depend only on $b_{1},b_{2},\sigma,k$, such that for all $c\geq c_{1}$ and $n\geq n_{0}(dc_{n})^{(2k+3)/(2k+2)}$, any solution in the constrained-form additive locally adaptive spline problem (34), with tuning parameter $\delta=c_{n}$, satisfies

[TABLE]

with probability at least $1-\exp(-c_{2}c)-c_{3}/n$.

• (b)

*Let $\mathcal{S}_{j}=\mathcal{G}_{j}$, the set of *$k$th degree splines with knots in the set $T_{j}$ in (27), for $j=1,\ldots,d$. There exist constants $c_{1},c_{2},c_{3},n_{0}>0$, that depend only on $b_{1},b_{2},\sigma,k$, such that for $c\geq c_{1}$ and $n\geq(dc_{n})^{(2k+3)/(2k+2)}$, any solution in the constrained-form restricted additive locally adaptive spline problem (34), with tuning parameter $\delta=a_{k}c_{n}$, where $a_{k}\geq 1$ is a constant that depends only on $k$, satisfies (41), with probability at least $1-\exp(-c_{2}c)-c_{3}d/n$.

• (c)

*Let $\mathcal{S}_{j}=\mathcal{H}_{j}$, the set of *$k$*th degree falling factorial functions defined over $X_{j}$ (the *$j$th dimension of input points), for $j=1,\ldots,d$. Then there are constants $c_{1},c_{2},c_{3},n_{0}>0$, depending only on $b_{1},b_{2},\sigma,k$, such that for all $c\geq c_{1}$ and $n\geq n_{0}(dc_{n})^{(2k+3)/(2k+2)}$, any solution in the constrained-form additive trend filtering problem (34), with tuning parameter $\delta=a_{k}c_{n}$, where $a_{k}\geq 1$ is a constant depending only on $k$, satisfies (41), with probability at least $1-\exp(-c_{2}c)-c_{3}d/n$.

3.4 Minimax lower bounds

We consider minimax lower bounds for estimation over the class of additive functions whose components are smooth with respect to the seminorm $J$. We allow the dimension $d$ to grow with $n$. As for the data distribution, we will use the following assumptions in place of Assumption A1 (Continuous input distribution), Assumption A2 (Bounded input density), Assumption A3 (Product input distribution), Assumption B1 (Generic regression model, with sub-Gaussian errors).

Assumption A4 (Uniform input distribution)****.

The inputs $X^{i}$, $i=1,\ldots,n$ are i.i.d. from the uniform distribution on $[0,1]^{d}$.

Assumption B2 (Additive model, Gaussian errors)****.

The responses $Y^{i}$, $i=1,\ldots,n$ follow

[TABLE]

with mean $\mu\in\mathbb{R}$, where $\int_{[0,1]^{d}}f_{0}(x)\,dx=0$ for identifiability. The errors $\epsilon^{i}$, $i=1,\ldots,n$ are i.i.d. $N(0,\sigma^{2})$, for some constant $\sigma>0$. The errors and input points are independent.

For the regularizer $J$, assumed to satisfy Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), we will replace Assumption Assumption C3 (Entropy bound) by the following assumption, on the log packing and log covering (entropy) numbers.

**Assumption C4 (Matching packing and

covering number bounds)****.**

There exist constants $0<w<2$ and $K_{1},K_{2}>0$ such that

[TABLE]

(To be clear, here $\|\cdot\|_{2}$ is the $L_{2}$ norm defined with respect to the uniform distribution on $[0,1]$.)

Let us introduce the notation

[TABLE]

Now we state our main minimax lower bound. The proof is given in Appendix A.13, A.14.

Theorem 3**.**

Assume Assumption A4 (Uniform input distribution), Assumption B2 (Additive model, Gaussian errors) on the data distribution, and Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C4 (Matching packing and covering number bounds) on the seminorm $J$. Then there exist constants $c_{0},n_{0}>0$, that depend only on $\sigma,k,L,K_{1},K_{2},w$, such that for all $c_{n}\geq 1$ and $n\geq n_{0}d^{1+w/2}c_{n}^{w}$, we have

[TABLE]

When we choose $J(g)=\mathrm{TV}(g^{(k)})$ as our regularizer, the additive function class $B_{J}^{d}(\delta)$ becomes

[TABLE]

and Theorem 3 implies the following result, whose proof is in Appendix A.15.

Corollary 3**.**

Assume Assumption A4 (Uniform input distribution), Assumption B2 (Additive model, Gaussian errors) on the data distribution. Assume further that $f_{0j}$, $j=1,\ldots,d$ are $k$ times weakly differentiable. Then there are constants $c_{0},n_{0}>0$, that depend only on $\sigma,k$, such that for all $c_{n}\geq 1$ and and $n\geq n_{0}d^{(2k+3)/(2k+2)}c_{n}^{1/(k+1)}$,

[TABLE]

Remark 10** **(Optimality for a fixed dimension $d$).

For a fixed $d$, the estimator defined by (28) is minimax rate optimal over the class of additive functions $f_{0}$ such that $\sum_{j=1}^{d}J(f_{0j})\leq C_{n}$. To see this, note that such a class of functions contains $B_{J}^{d}(C_{n}/d)$, therefore plugging $c_{n}=C_{n}/d$ into the right-hand side in (42) yields a lower bound rate of $n^{-2/(2+w)}C_{n}^{2w/(2+w)}$, which matches the upper bound rate in (30).

Furthermore, when $J(g)=\mathrm{TV}(g^{(k)})$, the lower bound rate given by plugging $c_{n}=C_{n}/d$ into the right-hand side in (43) is $n^{-(2k+2)/(2k+3)}C_{n}^{2/(2k+3)}$, matching the upper bound rate in (33). Hence additive locally adaptive regression splines, restricted additive locally adaptive regression splines, and additive trend filtering all achieve the minimax rate over the space of additive functions $f_{0}$ such that $\sum_{j=1}^{d}\mathrm{TV}(f_{0j}^{(k)})\leq C_{n}$.

Remark 11** **(Optimality for a growing dimension $d$).

For growing $d$, the estimator defined by (34) is minimax rate optimal over the class of additive functions $f_{0}$ such that $J(f_{0j})\leq c$, $j=1,\ldots,d$, where $c>0$ is a constant. This is verified by noting that the lower bound rate of $dn^{-2/(2+w)}$ in (42) matches the upper bound rates in (36), (37).

When $J(g)=\mathrm{TV}(g^{(k)})$, and again, $c_{n}=c$ (a constant), the lower bound rate of $dn^{-(2k+2)/(2k+3)}$ in (43) matches the upper bound rates in (41). Thus additive locally adaptive regression splines, restricted additive locally adaptive regression splines, and additive trend filtering all attain the minimax rate over the space of additive functions $f_{0}$ with $\mathrm{TV}(f_{0j}^{(k)})\leq c$, $j=1,\ldots,d$.

For growing $c_{n}$, we note that the upper bounds in (38) and (41) have an inflated dependence on $c_{n}$, compared to (42) and (43). It turns out that the latter (lower bounds) are tight, and the former (upper bounds) are loose. The upper bounds can be tightened under further boundedness assumptions (see Remark 7).

Remark 12** **(Suboptimality of additive linear smoothers).

Seminal theory from Donoho and Johnstone (1998) on minimax linear rates over Besov spaces shows that, under Assumption Assumption B2 (Additive model, Gaussian errors), and with the inputs $X^{i}$, $i=1,\ldots,n$ being now nonrandom and occurring over the regular $d$-dimensional lattice $\{1/N,2/N,\ldots,1\}^{d}\subseteq[0,1]^{d}$ with $N=n^{1/d}$, we have

[TABLE]

for all $n\geq n_{0}$, where $c_{0},n_{0}>0$ are constants, depending only on $\sigma,k$. On the left-hand side in (44) the infimum is taken over all additive linear smoothers, i.e., estimators $\hat{f}=\sum_{j=1}^{d}\hat{f}_{j}$ such that each component $\hat{f}_{j}$ is a linear smoother, for $j=1,\ldots,d$. The additive linear smoother lower bound (44) is verified in Appendix A.16.

For a fixed $d$, we can see that all additive linear smoothers—e.g., additive smoothing splines, additive kernel smoothing estimators, additive RKHS estimators, etc.—are suboptimal over the class of additive functions $f_{0}$ with $\sum_{j=1}^{d}\mathrm{TV}(f_{0j}^{(k)})\leq C_{n}$, as the optimal linear rate in (44) (set $c_{n}=C_{n}/d$) is $n^{-(2k+1)/(2k+2)}C_{n}^{2/(2k+2)}$, slower than the optimal rate $n^{-(2k+2)/(2k+3)}C_{n}^{2/(2k+2)}$ of additive locally adaptive regression splines and additive trend filtering in (33).

For growing $d$, and $c_{n}=c$ (a constant), we also see that additive linear smoothers are suboptimal over the class of additive functions $f_{0}$ such that $\mathrm{TV}(f_{0j}^{(k)})\leq c$, $j=1,\ldots,d$, as the optimal linear rate in (44) is $dn^{-(2k+1)/(2k+2)}$, slower than the optimal rate $dn^{-(2k+2)/(2k+3)}$ of additive locally adaptive regression splines and additive trend filtering in (41).

4 Backfitting and the dual

We now examine computational approaches for the additive trend filtering problem (12). This is a convex optimization problem, and many standard approaches can be applied. For its simplicity and its ubiquity in additive modeling, we focus on the backfitting algorithm in particular.

4.1 Backfitting

The backfitting approach for problem (12) is described in Algorithm 1. We write $\mathrm{TF}_{\lambda}(r,X_{j})$ for the univariate trend filtering fit, with a tuning parameter $\lambda>0$, to a response vector $r=(r^{1},\ldots,r^{n})\in\mathbb{R}^{n}$ over an input vector $X_{j}=(X^{1}_{j},\ldots,X^{n}_{j})\in\mathbb{R}^{n}$. In words, the algorithm cycles over $j=1,\ldots,d$, and at each step updates the estimate for component $j$ by applying univariate trend filtering to the $j$th partial residual (i.e., the current residual excluding component $j$). Centering in Step 2b part (ii) is optional, because the fit $\mathrm{TF}_{\lambda}(r,X_{j})$ will have mean zero whenever $r$ has mean zero, but centering can still be performed for numerical stability. In general, the efficiency of backfitting hinges on the efficiency of the univariate smoother employed; to implement Algorithm 1 in practice we can use fast interior point methods (Kim et al., 2009) or fast operator splitting methods (Ramdas and Tibshirani, 2016) for univariate trend filtering, both of which result in efficient empirical performance.

Algorithm 1 is equivalent to block coordinate descent (BCD), also called exact blockwise minimization, applied to problem (12) over the coordinate blocks $\theta_{j}$, $j=1,\ldots,d$. A general treatment of BCD is given in Tseng (2001), who shows that for a convex criterion that decomposes into smooth plus separable terms, as does that in (12), all limit points of the sequence of iterates produced by BCD are optimal solutions. A recent wave of work from the optimization community gives refined convergence analyses for coordinate descent (or its variants) in particular settings. We do not pursue the implications of this work for our problem; our interest here is primarily in developing a connection between BCD for problem (12) and alternating projections in its dual problem (24), which is the topic of the next subsection.

4.2 Dual alternating projections

Using the additive trend filtering problem (12) and its dual (24), related by the transformation (25), we see that for any dimension $j=1,\ldots,d$, the univariate trend filtering fit with response vector $r=(r^{1},\ldots,r^{n})$ and input vector $X_{j}=(X^{1}_{j},\ldots,X^{n}_{j})$ can be written as

[TABLE]

where $U_{j}=\{S_{j}D_{j}^{T}v_{j}:\|u\|_{\infty}\leq\lambda\}$, and recall, we abbreviate $D_{j}=D^{(X_{j},k+1)}$. (This follows from (45) specialized to the univariate trend filtering problem.) The backfitting approach in Algorithm 1 can be viewed (ignoring the optional centering step) as performing the updates, for $t=1,2,3,\ldots$,

[TABLE]

or, reparametrized in terms of the primal-dual relationship $u=Y-\bar{Y}\mathbbm{1}-\sum_{j=1}^{d}\theta_{j}$ in (25),

[TABLE]

Thus the backfitting algorithm for (12), as expressed above in (47), is seen to be a particular type of alternating projections method applied to the dual problem (24), cycling through projections onto $U_{j}$, $j=1,\ldots,d$. Interestingly, as opposed to the classical alternating projections approach, which would repeatedly project the current iterate $u_{j-1}^{(t)}$ onto $U_{j}$, $j=1,\ldots,d$, the steps in (47) repeatedly project an “offset” version $u_{j-1}^{(t)}+\theta_{j}^{(t-1)}$ of the current iterate, for $j=1,\ldots,d$ (this corresponds to running univariate trend filtering on the current residual, in the iterations (46)).

There is a considerable literature on alternating projections in optimization, see, e.g., Bauschke and Borwein (1996) for a review. Many alternating projections algorithms can be derived from the perspective of an operator splitting technique, e.g., the alternative direction method of multipliers (ADMM). Indeed, the steps in (47) appear very similar to those from an ADMM algorithm applied to the dual (24), if we think of the “offset” variables $\theta_{j}$, $j=1\ldots,d$ in the iterations (47) as dual variables in the dual problem (24) (i.e., if we think of the primal variables $\theta_{j}$, $j=1,\ldots,d$ as dual variables in the dual problem (24)). This connection inspires a new parallel version of backfitting, presented in the next subsection.

4.3 Parallelized backfitting

We have seen that backfitting is a special type of alternating projections algorithm, applied to the dual problem (24). For set intersection problems (where we seek a point in the intersection of given closed, convex sets), the optimization literature offers a variety of parallel projections methods (in contrast to alternating projections methods) that are provably convergent. One such method can be derived using ADMM (e.g., see Section 5.1 of Boyd et al. (2011)), and a similar construction can be used for the dual problem (24). We first rewrite this problem as

[TABLE]

where we write $I_{S}$ for the indicator function of a set $S$ (equal to 0 on $S$, and $\infty$ otherwise). Then we define the augmented Lagrangian, for an arbitrary $\rho>0$, as

[TABLE]

The ADMM steps for problem (48) are now given by repeating, for $t=1,2,3,\ldots$,

[TABLE]

Now compare (49) to (47)—the key difference is that in (49), the updates to $u_{j}$, $j=1,\ldots,d$, i.e., the projections onto $U_{j}$, $j=1,\ldots,d$, completely decouple and can hence be performed in parallel. Run properly, this could provide a large speedup over the sequential projections in (47).

Of course, for our current study, the dual problem (48) is really only interesting insofar as it is connected to the additive trend filtering problem (12). Fortunately, the parallel projections algorithm (49) maintains a very useful connection to the primal problem (12): $\rho\hat{\gamma}_{j}$, $j=1,\ldots,d$ (i.e., the scaled iterates $\rho\gamma_{j}^{(t)}$, $j=1,\ldots,d$ at convergence) are optimal for the additive trend filtering problem (12). This is simply because the dual of the dual problem (48) is the additive trend filtering problem (12) (therefore $\rho\gamma_{j}$, $j=1,\ldots,d$, which are dual to the constraints in (48), are equivalent to the primal parameters $\theta_{j}$, $j=1,\dots d$ in (12)). We state this next as a theorem, and transcribe the iterations in (49) into an equivalent primal form, in Algorithm 2. For details, see Appendix A.17.

Theorem 4**.**

Initialized arbitrarily, the ADMM steps (49) produce parameters $\hat{\gamma}_{j}$, $j=1,\ldots,d$ (i.e., the iterates $\gamma_{j}^{(t)}$, $j=1,\ldots,d$ at convergence) such that the scaled parameters $\rho\hat{\gamma}_{j}$, $j=1,\ldots,d$ solve additive trend filtering (12). Further, the outputs $\hat{\theta}_{j}$, $j=1,\ldots,d$ of Algorithm 2 solve additive trend filtering (12).

Written in primal form, we see that the the parallel backfitting approach in Algorithm 2 differs from what may be considered the “naive” approach to parallelizing the usual backfitting iterations in Algorithm 1. Consider $\rho=1$. If we were to replace Step 2a in Algorithm 2 with $u_{0}^{(t)}=r^{(t-1)}$, the full residual

[TABLE]

then the update steps for $\theta_{j}^{(t)}$, $j=1,\ldots,d$ that follow would be just given by applying univariate trend filtering to each partial residual (without sequentially updating the partial residuals between trend filtering runs). This naive parallel method has no convergence guarantees, and can fail even in simple practical examples to produce optimal solutions. Importantly, Algorithm 2 does not take $u_{0}^{(t)}$ to be the full residual, but as Step 2a shows, uses a less greedy choice: it basically takes $u_{0}^{(t)}$ to be a convex combination of the residual $r^{(t-1)}$ and its previous value $u_{0}^{(t-1)}$, with higher weight on the latter. The subsequent parallel updates for $\theta_{j}^{(t)}$, $j=1,\ldots,d$ are still given by univariate trend filtering fits, and though these steps do not exactly use partial residuals (since $u_{0}^{(t)}$ is not exactly the full residual), they are guaranteed to produce additive trend filtering solutions upon convergence (as per Theorem 4). An example of cyclic versus parallelized backfitting is given in Appendix A.18.

5 Experiments

Through empirical experiments, we examine the performance of additive trend filtering relative to additive smoothing splines. We also examine the efficacy of cross-validation for choosing the tuning parameter $\lambda$, as well as the use of multiple tuning parameters (i.e., a separate parameter $\lambda_{j}$ for each component $j=1,\ldots,d$). All experiments were performed in R. For the univariate trend filtering solver, we used the trendfilter function in the glmgen package, which is an implementation of the fast ADMM algorithm given in Ramdas and Tibshirani (2016). For the univariate smoothing spline solver, we used the smooth.spline function in base R.

5.1 Simulated heterogeneously-smooth data

We sampled $n=2500$ input points in $d=10$ dimensions, by assigning the inputs along each dimension $X_{j}=(X^{1}_{j},\ldots,X^{n}_{j})$ to be a different permutation of the equally spaced points $(1/n,2/n,\ldots,1)$, for $j=1,\ldots,10$. For the componentwise trends, we examined sinusoids with Doppler-like spatially-varying frequencies:

[TABLE]

We then defined the component functions as $f_{0j}=a_{j}g_{0j}-b_{j}$, $j=1,\ldots,d$, where $a_{j},b_{j}$ were chosen so that $f_{0j}$ had empirical mean zero and empirical norm $\|f_{0j}\|_{n}=1$, for $j=1,\ldots,d$. The responses were generated according to $Y^{i}\overset{\mathrm{i.i.d.}}{\sim}N(\sum_{j=1}^{d}f_{0j}(X^{i}_{j}),\sigma^{2})$, $i=1,\ldots,2500$. By construction, in this setup, there is considerable heterogeneity in the levels of smoothness both within and between the component functions.

The left panel of Figure 4 shows a comparison of the MSE curves from additive trend filtering in (12) (of quadratic order, $k=2$) and additive smoothing splines in (1) (of cubic order). We set $\sigma^{2}$ in the generation of the responses so that the signal-to-noise ratio (SNR) was $\|f_{0}\|_{n}^{2}/\sigma^{2}=4$, where $f_{0}=\sum_{j=1}^{d}f_{0j}$. The two methods (additive trend filtering and additive smoothing splines) were each allowed their own sequence of tuning parameter values, and results were averaged over 10 repetitions from the simulation setup described above. As we can see, additive trend filtering achieves a better minimum MSE along its regularization path, and does so at a less complex model (lower df).

The right panel of Figure 4 shows the best-case MSEs for additive trend filtering and additive smoothing splines (i.e., the minimum MSE over their regularization paths) as the noise level $\sigma^{2}$ is varied so that the SNR ranges from 0.7 to 1.6, in equally spaced values on the log scale. The results were again averaged over 10 repetitions of data drawn from a simulation setup essentially the same as the one described above, except that we considered a smaller problem size, with $n=1000$ and $d=6$. The plot reveals that additive trend filtering performs increasingly well (in comparison to additive smoothing splines) as the SNR grows, which is not surprising, since for high SNR levels it is able to better capture the heterogeneity in the component functions.

Lastly, in Appendix A.19, we present results from an experimental setup mimicking that in this subsection, except with the component functions $f_{0j}$, $j=1,\ldots,d$ having homogeneous smoothness throughout. The results show that additive trend filtering and additive smoothing splines perform nearly exactly the same.

5.2 Cross-validation for tuning parameter selection

Sticking to the simulation setup from the last subsection, but at the smaller problem size, $n=1000$ and $d=6$ (used to produce the right panel of Figure 4), we study in the left panel of Figure 5 the use of 5-fold cross-validation (CV) to select the tuning parameter $\lambda$ for additive trend filtering and additive smoothing splines. Displayed are the resulting MSE curves as the SNR varies from 0.7 to 16. Also shown on the same plot are the oracle MSE curves (which are the same as those the right panel of Figure 4), in which $\lambda$ has been chosen to minimize the MSE for each method. We can see that the performance of each method degrades using CV, but not by much.

In the right panel of the figure, we examine the use of multiple tuning parameters for additive smoothing splines and additive trend filtering, i.e., replacing the penalties in (1) and (12) by those in (13), respectively, so we now have $d$ tuning parameters $\lambda_{j}$, $j=1,\ldots,d$. When the function we are estimating has different amounts of smoothness along different dimensions, we have argued (and seen through examples) that additive trend filtering—using only a single tuning parameter $\lambda$—can accomodate these differences, at least somewhat, thanks to its locally adaptive nature. But, when the difference in smoothness across dimensions is drastic enough, it may be worthwhile to move to multiple tuning parameters $\lambda_{j}$, $j=1,\ldots,d$. Certainly, moving to multiple tuning parameters will help additive smoothing splines address heterogeneity in smoothness across dimensions, which will have greater difficulty in accommodating such heterogeneity using a single tuning parameter, compared to additive trend filtering.

When $d$ is moderate (even just for $d=6$), cross-validation over a $d$-dimensional grid of values for $\lambda_{j}$, $j=1,\ldots,d$ can be prohibitive. As discussed in Section 1.5, there is a lot of literature dedicated to an alternative approach based on restricted maximum likelihood (REML), mostly surrounding additive models built from linear smoothers. Unfortunately, as far as we understand, REML does not easily apply to additive trend filtering. We thus use the following simple approach for multiple tuning parameter selection: within each backfitting loop, for each component $j=1,\ldots,d$, we use (univariate) CV to choose $\lambda_{j}$. While this does not solve a particular convex optimization problem, and is not guaranteed to converge in general, we have found it to work quite well in practice. (We have helpful to use CV to find the best choice of single tuning parameter for additive trend filtering, and then use this estimate to initialize the backfitting routine.) The right panel of Figure 5 compares the performance of this so-called backfit-CV tuning to the oracle, that chooses just a single tuning parameter. Both additive trend filtering and additive smoothing splines are seen to improve with $d$ tuning parameters, tuned by backfit-CV, in comparison to the oracle choice of tuning parameter. Interestingly, we also see that additive smoothing splines with $d$ tuning parameters performs on par with additive trend filtering with the oracle choice of tuning parameter. (In this example, REML tuning for additive smoothing splines—as implemented by the mgcv R package—performed a bit worse than backfit-CV tuning, and so we only show results from the latter.)

6 Discussion

We have studied additive models built around the univariate trend filtering estimator, i.e., defined by penalizing according to the sum of $\ell_{1}$ norms of discrete derivatives of the component functions. We examined basic properties of these additive models, such as extrapolation of the fitted values to a $d$-dimensional surface, uniqueness of the component fits, and characterization of effective degrees of freedom of the fit. When the underlying regression function is additive, with components whose $k$th derivatives are of bounded variation, we derived error rates for $k$th order additive trend filtering: $n^{-(2k+2)/(2k+3)}$ for a fixed input dimension $d$ (under weak assumptions), and $dn^{-(2k+2)/(2k+3)}$ for a growing dimension $d$ (under stronger assumptions). We showed these rates are sharp by establishing matching minimax lower bounds, and showed that additive linear smoothers (e.g., additive smoothing splines) can at best achieve a rate of $n^{-(2k+1)/(2k+2)}$ for fixed $d$, and $dn^{-(2k+1)/(2k+2)}$ for growing $d$, over the same function class. Finally, on the computational side, we devised a provably convergent parallel backfitting algorithm for additive trend filtering.

It is worth noting that our parallel backfitting method is not specific to additive trend filtering, but it can be embedded in a more general parallel coordinate descent framework (Tibshirani, 2017). It would be interesting to follow up on this framework to see what practical benefits it can provide for additive models in particular.

A natural natural extension of our work is to consider the high-dimensional case, where $d$ is comparable or possibly even much larger than $n$, and we fit a sparse additive model by employing an additional sparsity penalty in problem (12). Another natural extension is to consider responses $Y^{i}|X^{i}$, $i=1,\ldots,n$ from an exponential family distribution, and we fit a generalized additive model by changing the loss in (12). After we completed an initial version of this paper, both extensions have been pursued: Tan and Zhang (2017) develop a suite of error bounds for sparse additive models, with various form of penalties (which include total variation on derivatives of components); and Haris et al. (2018) give comprehensive theory for sparse generalized additive models, with various types of penalties (which again include total variation on derivatives of components).

Acknowledgements

We are very thankful to Garvesh Raskutti for his generous help and insight on various issues, and Martin Wainwright for generously sharing his unpublished book with us. We are also grateful to an anonymous referee whose thoughtful comments improved our paper.

Appendix A Appendix

A.1 Fast extrapolation

We discuss extrapolation using the fitted functions $\hat{f}_{j}$, $j=1,\ldots,d$ from additive trend filtering (16), as in (19). We must compute the coefficients $\hat{\alpha}_{j}=(\hat{a}_{j},\hat{b}_{j})$ whose block form is given in (20), (21). Clearly, the computation of $\hat{b}_{j}$ in (21) requires $O(n)$ operations (owing to the bandedness of $D^{(X_{j},k+1)}$, and treating $k$ as a constant). As for $\hat{a}_{j}$ in (20), it can be seen from the structure of $C^{(X_{j},k+1)}$ as described in Wang et al. (2014) that

[TABLE]

which takes only $O(1)$ operations (again treating $k$ as constant, and now using the bandedness of each $D_{j}^{(X_{j},\ell-1)}$, $\ell=2,\ldots,k+1$). In total then, computing the coefficients $\hat{\alpha}_{j}=(\hat{a}_{j},\hat{b}_{j})$ requires $O(n)$ operations, and computing $\hat{\alpha}=(\hat{\alpha}_{1},\ldots,\alpha_{d})$ requires $O(nd)$ operations.

After having computed $\hat{\alpha}=(\hat{\alpha}_{1},\ldots,\alpha_{d})$, which only needs to be done once, a prediction at a new point $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ with the additive trend filtering fit $\hat{f}$ is given by

[TABLE]

This requires $O(d+\sum_{j=1}^{d}\sum_{\ell=k+2}^{n}1\{\hat{\alpha}^{\ell}_{j}\neq 0\})$ operations, utilizing the sparsity of the components in $\hat{\alpha}$ not associated with the polynomial basis functions.

A.2 Proof of Lemma 2

We begin by eliminating the constraint in the additive trend filtering problem (12), rewriting it as

[TABLE]

where $M=I-\mathbbm{1}\mathbbm{1}^{T}/n$. Noting that $D^{(X_{j},k+1)}\mathbbm{1}=0$ for $j=1,\ldots,d$, we can replace the penalty term above by $\sum_{j=1}^{d}\|D^{(X_{j},k+1)}S_{j}\theta_{j}\|_{1}$. Reparametrizing using the falling factorial basis, as in Lemma 1, yields the problem

[TABLE]

where we have used the abbreviation $P_{j}=P^{(X_{j},k)}$ and $K_{j}=K^{(X_{j},k)}$, as well as the block representation $\alpha_{j}=(a_{j},b_{j})\in\mathbb{R}^{(k+1)}\times\mathbb{R}^{(n-k-1)}$, for $j=1,\ldots,d$. Since each $P_{j}$, $j=1,\ldots d$ has $\mathbbm{1}$ for its first column, the above problem is equivalent to

[TABLE]

where $\tilde{P}_{j}$ denotes $P_{j}$ with the first column removed, for $j=1,\ldots,d$. To be clear, solutions in the above problem and the original trend filtering formulation (12) are related by

[TABLE]

Furthermore, we can see that $\hat{a}=(\hat{a}_{1},\ldots\hat{a}_{d})$ solves

[TABLE]

where $\tilde{P}$ is as defined in (22), and $\hat{b}=(\hat{b}_{1},\ldots,\hat{b}_{d})$ solves

[TABLE]

where $UU^{T}$ is the projection orthogonal to the column space of $\tilde{P}$, i.e., it solves

[TABLE]

where $\tilde{K}$ is as in (23). Since problem (51) is a standard lasso problem, existing results on the lasso (e.g., Tibshirani (2013)) imply that the solution $\hat{b}$ is unique whenever $\tilde{K}$ has columns in general position. This proves the first part of the lemma. For the second part of the lemma, note that the solution $\hat{a}$ in the least squares problem (50) is just given by the regression of $MY-M\sum_{j=1}^{d}K_{j}\hat{b}_{j}$ onto $\tilde{P}$, which is unique whenever $\tilde{P}$ has full column rank. This completes the proof.

A.3 Derivation of additive trend filtering dual

As in the proof of Lemma 2, we begin by rewriting the problem (12) as

[TABLE]

where $M=I-\mathbbm{1}\mathbbm{1}^{T}/n$. Then, we reparametrize the above problem,

[TABLE]

and form the Lagrangian

[TABLE]

Minimizing the Lagrangian $L$ over all $\theta,z$ yields the dual problem

[TABLE]

The claimed dual problem (24) is just the above, rewritten in an equivalent form.

A.4 Proof of Lemma 3

We first eliminate the equality constraint in (12), rewriting this problem, as was done in the proof of Lemma 2, as

[TABLE]

where $M=I-\mathbbm{1}\mathbbm{1}^{T}/n$, and $D_{j}=D^{(X_{j},k+1)}$, $j=1,\ldots,d$. This is a generalized lasso problem with a design matrix $T\in\mathbb{R}^{n\times nd}$ that has $d$ copies of $M$ stacked along its columns, and a penalty matrix $D\in\mathbb{R}^{nd\times nd}$ that is block diagonal in the blocks $D_{j}$, $j=1,\ldots,d$. Applying Theorem 3 of Tibshirani and Taylor (2012), we see that

[TABLE]

where $A=\mathrm{supp}(D\hat{\theta})$, and where $D_{-A}$ denotes the matrix $D$ with rows removed that correspond to the set $A$. The conditions for uniqueness in the lemma now precisely imply that

[TABLE]

where $A_{j}$ denotes the subset of $A$ corresponding to the block of rows occupied by $D_{j}$, and $|A_{j}|$ its cardinality, for $j=1,\ldots,d$. This can be verified by transforming to the basis perspective as utilized in the proofs of Lemmas 1 and 2. The desired result is obtained by noting that, for $j=1,\ldots,d$, the component $\hat{\theta}_{j}$ exhibits a knot for each element in $A_{j}$.

A.5 Preliminaries for the proof of Theorem 1

Before the proof of Theorem 1, we collect important preliminary results. We start with a result on orthonormal polynomials. We thank Dejan Slepcev for his help with the next lemma.

Lemma 4**.**

Given an integer $\kappa\geq 0$, and a continuous measure $\Lambda$ on $[0,1]$, whose Radon-Nikodym derivative $\lambda$ is bounded below and above by constants $b_{1},b_{2}>0$, respectively. Denote by $\phi_{0},\ldots,\phi_{\kappa}$ an orthonormal basis for the space of polynomials of degree $\kappa$ on $[0,1]$, given by running the Gram-Schmidt procedure on the polynomials $1,t,\ldots,t^{\kappa}$, with respect to the $L_{2}(\Lambda)$ inner product. Hence, for $\ell=0,\ldots,\kappa$, $\phi_{\ell}$ is an $\ell$th degree polynomial, orthogonal (in $L_{2}(\Lambda)$) to all polynomials of degree less than $\ell$, and we denote its leading coefficient by $a_{\ell}>0$. Now define, for $t\in[0,1]$:

[TABLE]

Then the following two relations hold:

[TABLE]

and

[TABLE]

Proof.

First, we use induction to show that for $t\in[0,1]$,

[TABLE]

This statement holds for $\ell=1$ by definition of $\Phi_{\kappa,0},\Phi_{\kappa,1}$. Assume it holds at some $\ell>1$. Then

[TABLE]

where we used inductive hypothesis in the first line and Fubini’s theorem in the second line, which completes the inductive proof.

Now, the relation in (54) shows that $\Phi_{\kappa,\ell}(1)$ is the $L_{2}(\Lambda)$ inner product of $\phi_{\kappa}$ and an $(\ell-1)$st degree polynomial, for $\ell=1,\ldots,\kappa$. As $\phi_{\kappa}$ is orthogonal to all polynomials of degree less than $\kappa$, we have $\Phi_{\kappa,\ell}(1)=0$, $\ell=1,\ldots,\kappa$. For $\ell=\kappa+1$, note that this same orthogonality along with (54) also shows

[TABLE]

where $\langle\cdot,\cdot\rangle_{2}$ is the $L_{2}(\Lambda)$ inner product. This establishes the statement in (52).

As for (53), note that if $\kappa=0$ then the statement holds, because the left-hand side is 1 and the right-hand side is always larger than 1. Hence consider $\kappa\geq 1$. From (54), we have, for any $t\in[0,1]$,

[TABLE]

where in the second line we used Cauchy-Schwartz, and in the third line we used the fact that $\phi_{\kappa}$ has unit norm, and the upper bound $b_{2}$ on $\lambda$. Next we bound $a_{\kappa}$. Let $p$ be the projection of $x^{\kappa}$ onto the space of polynomials of degree $\kappa-1$, with respect to the $L_{2}(\Lambda)$ inner product. Then we have $\phi_{\kappa}=(x^{\kappa}-p)/\|x^{\kappa}-p\|_{2}$, thus its leading coefficient is $a_{\kappa}=1/\|x^{\kappa}-p\|_{2}$, where $\|\cdot\|_{2}$ is the $L_{2}(\Lambda)$ norm. Consider

[TABLE]

In the first line we used the lower bound $b_{1}$ on $\lambda$. In the second we used the fact the the Legendre polynomial $P_{\kappa}$ of degree $\kappa$, shifted to $[0,1]$ but unnormalized, is the result from projecting out $1,t,\ldots,t^{\kappa-1}$ from $t^{\kappa}$, with respect to the uniform measure. In the last step we used the fact that $P_{\kappa}$ has norm $1/(\sqrt{2\kappa+1}{2\kappa\choose\kappa})$. Combining (55) and (A.5) gives the result (53). ∎

Remark 13** **(**Special case: uniform measure and Rodrigues’

formula**).

In the case of $\Lambda=U$, the uniform measure on $[0,1]$, we can just take $\phi_{0},\ldots,\phi_{\kappa}$ to be the Legendre polynomials, shifted to $[0,1]$ and normalized appropriately. Invoking Rodrigues’ formula to express these functions,

[TABLE]

*the results in Lemma 4 can be directly verified. *

We use Lemma 4 to construct a sup norm bound on functions in $B_{J}(1)$ that are orthogonal (in $L_{2}(\Lambda)$) to all polynomials of degree $k$. We again Dejan Slepcev for his help with the next lemma.

Lemma 5**.**

Given an integer $k\geq 0$, and a continuous measure $\Lambda$ on $[0,1]$, whose Radon-Nikodym derivative $\lambda$ is bounded below and above by constants $b_{1},b_{2}>0$, respectively. Let $J$ be a functional satisfying Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials) and Assumption C2 (Relative boundedness of derivatives), for a constant $L>0$. There is a constant $R_{0}>0$, that depends only on $k,b_{1},b_{2},L$, such that

[TABLE]

where $\langle\cdot,\cdot\rangle_{2}$ denotes the $L_{2}(\Lambda)$ inner product.

Proof.

Fix an arbitrary function $g\in B_{J}(1)$, orthogonal (in $L_{2}(\Lambda)$) to all polynomials of degree $k$. Using integration by parts, and repeated application of Lemma 4, we have

[TABLE]

where $w_{\ell}(t)=(-1)^{\ell}a_{\ell}\ell!\Phi_{\ell,\ell}(t)$, $\ell=0,\ldots,k$, and by properties (52), (53) of Lemma 4,

[TABLE]

Now, we will prove the following by induction:

[TABLE]

Starting at $\ell=k$, the statement holds because, using (57), for almost every $t\in[0,1]$,

[TABLE]

where in the second line we used the fact that the weight function integrates to 1 from (58), and in the third we used Assumption Assumption C2 (Relative boundedness of derivatives) and the upper bound on the integrated absolute weights from (58). Assume the statement holds at some $\ell<k$. Then again by (57), (58), for almost every $t\in[0,1]$,

[TABLE]

the last line using $\mathop{\mathrm{ess\,\,sup}}_{0\leq u<v\leq 1}|\int_{u}^{v}g^{(\ell)}(s)\,ds|\leq\|g^{(\ell)}\|_{\infty}$ and the inductive hypothesis. This verifies (59). Taking $\ell=0$ in (59) and defining $R_{0}=L(b_{2}/b_{1})^{(k+1)/2}\prod_{i=0}^{k}{2i\choose i}$ proves the lemma. ∎

We study the minimum eigenvalue of the (uncentered) empirical covariance matrix of a certain basis for additive $k$th degree polynomials in $\mathbb{R}^{kd}$. We thank Mathias Drton for his help with part (a) of the next lemma.

Lemma 6**.**

Let $X^{i}$, $i=1,\ldots,n$ denote an i.i.d. sample from a continuous distribution $Q$ on $[0,1]^{d}$. For an integer $k\geq 0$, let $V\in\mathbb{R}^{n\times kd}$ be a matrix whose $i$th row is given by

[TABLE]

for $i=1,\ldots,n$. Let

[TABLE]

where $\lambda_{\min}(\cdot)$ denotes the minimum eigenvalue of its argument. Assuming that $n\geq kd$, the following properties hold:

• (a)

$\nu_{n}>0$, almost surely with respect to $Q$;

• (b)

$\nu_{0}>0$;

• (c)

for any $0\leq t\leq 1$, $\mathbb{P}(\nu_{n}^{2}>t\nu_{0}^{2})$ with probability at least \displaystyle 1-d\exp\Big{(}-\frac{(1-t)^{2}\nu_{0}n}{2(kd)^{2}}\Big{)}.

Proof.

For part (a), if the claim holds for $n=kd$, then it holds for all $n>kd$, so we may assume without a loss of generality that $n=kd$. Note that the determinant of $V\in\mathbb{R}^{n\times kd}$ is a polynomial function, call it $q(X)$, of the elements $X^{i}_{j}$, $i=1,\ldots,n$, $j=1,\ldots,d$. By Lemma 1 of Okamoto (1973), the roots of any polynomial—that is not identically zero—form a set of Lebesgue measure zero. To check that the polynomial $q$ in question is not identically zero, it suffices to show that it is nonzero at a single realization of $X$. To this end, consider an input matrix defined by

[TABLE]

the rowwise concatenation of $\alpha_{\ell}I\in\mathbb{R}^{d\times d}$, $\ell=1,\ldots,k$. By the blockwise Vandermonde structure of the corresponding basis matrix $V$, we have that $q(X)\not=0$ provided the coefficients $\alpha_{\ell}$, $\ell=1,\ldots,k$ are all distinct. Therefore $q$ is not identically zero, and with respect to the continuous distribution $Q$, the determinant of $V$ is nonzero, i.e., $\nu_{n}>0$, almost surely.

For part (b), given any $a\in\mathbb{R}^{kd}$ with $a\not=0$, we know that $Va\not=0$ almost surely, since $\nu_{n}>0$ almost surely, by part (a). Thus

[TABLE]

which proves that $\nu_{0}>0$.

Part (c) is an application of a matrix Chernoff bound from Tropp (2012). In order to apply this result, we must obtain an almost sure upper bound $R$ on $\lambda_{\max}(V^{i}(V^{i})^{T})$, with $V^{i}$ as in (60) and $\lambda_{\max}(\cdot)$ denoting the maximum eigenvalue of its argument. This follows as

[TABLE]

as each component of $V^{i}$ has absolute magnitude at most 1 (recalling that $Q$ is supported on $[0,1]^{d}$). Taking $R=(kd)^{2}$ and applying Corollary 5.2 of Tropp (2012) (to be specific, applying the form of the Chernoff bound given in Remark 5.3 of this paper) gives the result. ∎

The next lemma pertains to the additive function space

[TABLE]

We give a sup norm bound on the components of functions in $\mathcal{M}_{n}(1)\cap B_{n}(\rho)$. The proof combines Lemmas 5 and 6, and uses a general strategy that follows the arguments given in Example 2.1(ii) of van de Geer (1990).

Lemma 7**.**

Let $X^{i}$, $i=1,\ldots,n$ denote an i.i.d. sample from a continuous distribution $Q$ on $[0,1]^{d}$, and let $J$ be a seminorm satisfying Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials) and Assumption C2 (Relative boundedness of derivatives). There are constants $R_{1},R_{2},c_{0},n_{0}>0$, depending only on $d,k,L$, such that for all $\rho>0$ and $n\geq n_{0}$,

[TABLE]

with probabilty at least $1-\exp(-c_{0}n)$, where $\mathcal{M}_{n}(1)$ is the function space in (61).

Proof.

Fix an arbitrary $m=\sum_{j=1}^{d}m_{j}\in\mathcal{M}_{n}(1)\cap B_{n}(\rho)$. For each $j=1,\ldots,d$, decompose

[TABLE]

where $p_{j}$ is a polynomial of degree $k$ such that $\langle p_{j},1\rangle_{n}=0$, and $g_{j}$ is orthogonal to all polynomials of degree $k$, with respect to the $L_{2}(U)$ inner product, with $U$ the uniform distribution on $[0,1]$; in fact, by definition of $\mathcal{M}_{n}(1)$, we know that $\langle m_{j},1\rangle_{n}=0$ so

[TABLE]

By the triangle inequality and Lemma 5 applied to the measure $\Lambda=U$ (whose density is of course lower and upper bounded with $b_{1}=b_{2}=1$), we have, for each $j=1,\ldots,d$,

[TABLE]

where $R_{0}>0$ is the constant from Lemma 5, and we have used $J(m_{j})=J(g_{j})$, for $j=1,\ldots,d$, as the null space of $J$ contains $k$th degree polynomials.

The triangle inequality and (62) now imply

[TABLE]

Write

[TABLE]

for some coefficients $\alpha_{j\ell}$, $j=1,\ldots$, $\ell=1,\ldots,k$. For $V\in\mathbb{R}^{n\times kd}$ the basis matrix as in Lemma 6, and $\alpha=(\alpha_{11},\ldots,\alpha_{1k},\ldots,\alpha_{d1},\ldots,\alpha_{dk})\in\mathbb{R}^{kd}$, we have

[TABLE]

Furthermore, noting

[TABLE]

we have

[TABLE]

where $\nu_{n}^{2}=\lambda_{\min}(V^{T}V/n)$, as in Lemma 6, and we have used the upper bound in (63). Using part (c) of Lemma 6, with $t=1/2$, we have

[TABLE]

with probability at least $1-d\exp(-\nu_{0}n/(8(kd)^{2}))$, where $\nu_{0}^{2}=\lambda_{\min}(\mathbb{E}[V^{T}V]/n)$, as in Lemma 6. Therefore, using the triangle inequality and the fact that $Q$ is supported on $[0,1]^{d}$, we have for each $j=1,\ldots,d$, and any $x_{j}\in[0,1]$,

[TABLE]

with probability at least $1-d\exp(-\nu_{0}n/(8(kd)^{2}))$. Finally, for each $j=1,\ldots,d$, using the triangle inequality, and the sup norm bound from Lemma 5 once again,

[TABLE]

with probability $1-d\exp(-\nu_{0}n/(8(kd)^{2}))$, completing the proof. ∎

We give a simple bound on the entropy of an arbitrary sum of sets in terms of the entropies of the original sets.

Lemma 8**.**

Given sets $S_{1},\ldots,S_{d}$ and a norm $\|\cdot\|$, it holds that

[TABLE]

Proof.

For $j=1,\ldots,d$, suppose that $S_{j}$ can be covered in $N_{j}$ balls of radius $\delta/d$, with centers at $s_{j}^{1},\ldots,s_{j}^{N_{j}}$. Take an arbitrary $s\in S_{1}+\cdots+S_{d}$, and write $s=\sum_{j=1}^{d}s_{j}$, with $s_{j}\in S_{j}$, $j=1,\ldots,d$. For each $j=1,\ldots,d$, there is some $s_{j}^{\ell_{j}}$ such that $\|s_{j}-s_{j}^{\ell}\|\leq\delta/d$, and so by the triangle inequality

[TABLE]

That is, we have shown that $\prod_{j=1}^{d}N_{j}$ balls of radius $\delta$ with centers at

[TABLE]

cover $S$. This completes the proof. ∎

The next result represents our main tool from empirical process theory that will be used in the proof of Theorem 1. It is essentially an application of Lemma 3.5 in van de Geer (1990) (see also van de Geer (2000)).

Lemma 9**.**

Let $X^{i}$, $i=1,\ldots,n$ denote an i.i.d. sample from a continuous distribution $Q$ on $[0,1]^{d}$. Let $\epsilon^{i}$, $i=1,\ldots,n$ be uniformly sub-Gaussian random variates that have variance proxy $\sigma^{2}>0$ and are independent of $X^{i}$, $i=1,\ldots,n$. Let $J$ be a seminorm satisfying Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C3 (Entropy bound), and let $\rho>0$ be arbitrary. Then there are constants $c_{1},c_{2},c_{3},n_{0}>0$, depending only on $d,\sigma,k,L,K,w,\rho$, such that for all $c\geq c_{1}$ and $n\geq n_{0}$,

[TABLE]

with probabilty at least $1-\exp(-c_{2}c^{2})-\exp(-c_{3}n)$.

Proof.

Let $\Omega_{1}$ denote the event on which the conclusion in Lemma 7 holds, which has probability at least $1-\exp(-c_{3}n)$ for $n\geq n_{1}$, for constants $c_{3},n_{1}>0$. Also let $R_{0}=R_{1}\rho+R_{2}$, where $R_{1},R_{2}>0$ are the constants from the lemma. Denote

[TABLE]

On $\Omega_{1}$, consider

[TABLE]

The first inequality (64) uses the sup norm bound from Lemma 7; the second inequality (65) uses

[TABLE]

and applies Lemma 8 to the space on the right-hand side above. The third inequality (66) just uses the fact we may assume $R_{0}\geq 1$, without a loss of generality; and the last inequality (67) uses Assumption Assumption C3 (Entropy bound). The entropy bound established in (67) allows us to apply Lemma 3.5 van de Geer (1990) (see also Lemma 8.4 in van de Geer (2000)), which gives constants $c_{1},c_{2},n_{2}>0$, depending only on $d,\sigma,k,R_{0},K,w$, such that for all $c\geq c_{1}$ and $n\geq n_{1}$,

[TABLE]

on an event $\Omega_{2}$ with probability at least $1-\exp(-c_{2}c^{2})$. The desired result in the lemma therefore holds for all $c\geq c_{1}$ and $n\geq n_{0}=\max\{n_{1},n_{2}\}$, on $\Omega_{1}\cap\Omega_{2}$. ∎

We finish with two simple results, on shifting around exponents in sums and products.

Lemma 10**.**

For any $a,b\geq 0$, and any $0<q<1$,

[TABLE]

Proof.

The function $f(t)=(1+t)^{q}-(1+t^{q})$ has derivative $f^{\prime}(t)=q(1+t)^{q-1}-qt^{q-1}<0$ for all $t>0$, and so $f(t)<f(0)=0$ for all $t>0$. Plugging in $t=a/b$ and rearranging gives the claim. ∎

Lemma 11**.**

For any $a,b\geq 0$, and any $w$,

[TABLE]

Proof.

Note that either $ab^{1-w/2}\leq a^{1/(1+w/2)}b$ or $ab^{1-w/2}\geq a^{1/(1+w/2)}b$, and in the latter case we get $b\leq a^{1/(1+w/2)}$, so $ab^{1-w/2}\leq a^{2/(1+w/2)}$. ∎

A.6 Proof of Theorem 1

This proof roughly follows the ideas in the proof of Theorem 9 in Mammen and van de Geer (1997), though it differs in a few key ways. We use $c>0$ to denote a constant that will multiply our final estimation error bound; it will also control the probability with which our final result holds. Some steps will only hold for sufficiently large $n$, but we do not always make this explicit. Lastly, we will occasionally abuse our notation for the empirical norms and empirical inner products by using them with vector arguments, to be interpreted in the appropriate sense (e.g., $\langle m,v\rangle_{n}=\frac{1}{n}\sum_{i=1}^{n}v^{i}m(X^{i})$ for a function $m$ and vector $v\in\mathbb{R}^{n}$).

We break down the presentation of our proof into mini sections for readability.

Deriving a basic inequality. Denote by $\hat{f}=\sum_{j=1}^{d}\hat{f}_{j}$ the total additive fit in (28). Let $\mathcal{S}$ denote feasible set in (28). For any $f\in\mathcal{S}$, note that by orthogonality,

[TABLE]

where $\bar{\epsilon}=\frac{1}{n}\sum_{i=1}^{n}\epsilon^{i}$. Therefore $\hat{f}$ must also be optimal for the problem

[TABLE]

where $W^{i}=f_{0}(X^{i})+\epsilon^{i}-\bar{\epsilon}$, $i=1,\ldots,n$, and we denote $\lambda_{n}=\lambda/n$ and $J_{d}(f)=\sum_{j=1}^{d}J(f_{j})$. Standard arguments (from first-order optimality) show that any solution $\hat{f}$ in the above satisfies

[TABLE]

for any feasible $\tilde{f}=\sum_{j=1}^{d}\tilde{f}_{j}\in\mathcal{S}$. Expanding the definition of $W$ and rearranging gives

[TABLE]

Using the polarization identity $\langle a,b\rangle=\frac{1}{2}(\|a\|^{2}+\|b\|^{2}-\|a-b\|_{2}^{2})$ for an inner product $\langle\cdot,\cdot\rangle$ and its corresponding norm $\|\cdot\|$,

[TABLE]

Abbreviating $\hat{\Delta}=\hat{f}-\tilde{f}$, $\hat{J}=J_{d}(\hat{f})$, and $\tilde{J}=J_{d}(\tilde{f})$, and using $\langle\bar{\epsilon}\mathbbm{1},\hat{\Delta}\rangle=0$, this becomes

[TABLE]

which is our basic inequality. In what follows, we will restrict our attention to feasible $\tilde{f}$ such that $\|\tilde{f}-f_{0}\|_{n}\leq\max\{C_{n},\tilde{J}\}$.

Localizing the error vector. We prove that $\hat{\Delta}$ is appropriately bounded in the empirical norm. By the tail bound for quadratic forms of sub-Gaussian random variates in Theorem 2.1 of Hsu et al. (2012), for all $t>0$,

[TABLE]

and hence taking $t=\sqrt{n}$,

[TABLE]

on an event $\Omega_{1}$ with probability at least $1-\exp(-\sqrt{n})$. Thus returning to the basic inequality (68), using the Cauchy-Schwartz inequality, and the above bound, we have

[TABLE]

on $\Omega_{1}$. This is a quadratic inequality of the form $x^{2}\leq bx+c$ in $x=\|\hat{\Delta}\|_{n}$, so we can upper bound $x$ by the larger of the two roots, $x\leq(b+\sqrt{b^{2}+4c})/2\leq b+\sqrt{c}$, i.e.,

[TABLE]

on $\Omega_{1}$. Abbreviating $J^{*}=\max\{C_{n},\tilde{J}\}$, and using $J^{*}\geq 1$ (as $C_{n}\geq 1$ by assumption),

[TABLE]

on $\Omega_{1}$. Recalling $\|\tilde{f}-f_{0}\|_{n}\leq J^{*}$, and using the fact that $\lambda_{n}=o(1)$ for our eventual choice of $\lambda_{n}$, we have that for sufficiently large $n$,

[TABLE]

on $\Omega_{1}$.

Bounding the sub-Gaussian complexity term. We focus on the first term on the right-hand side in (68), i.e., the sub-Gaussian complexity term. Let $m=\hat{\Delta}/(\hat{J}+J^{*})$. By construction, we have $J(m)\leq 1$, and from (69), we have $\|m\|_{n}\leq 2\sqrt{5}\sigma+\sqrt{2}$ on $\Omega_{1}$. Then, applying Lemma 9, with the choice $\rho=2\sqrt{5}\sigma+\sqrt{2}$, we see that there are constants $c_{1},c_{2},c_{3}>0$ such that for all $c\geq c_{1}$,

[TABLE]

on $\Omega_{1}\cap\Omega_{2}$, where $\Omega_{2}$ is an event with probability at least $1-\exp(-c_{2}c^{2})-\exp(-c_{3}n)$. Plugging this into (68) gives

[TABLE]

on $\Omega_{1}\cap\Omega_{2}$. By the inequality in Lemma 11, applied to the first term on the right-hand side above, with $a=n^{-1/2}(\hat{J}+J^{*})^{w/2}$ and $b=\|\hat{\Delta}\|_{n}$,

[TABLE]

on $\Omega_{1}\cap\Omega_{2}$, where we abbreviate $r_{n}=n^{-1/(2+w)}$. Applying the simple inequality $2ab\leq a^{2}+b^{2}$ to the first term on the right-hand side, with $a=cr_{n}(\hat{J}+J^{*})^{w/(2+w)}$ and $b=\|\hat{\Delta}\|_{n}$, and subtracting $\|\hat{\Delta}\|_{n}^{2}/2$ from both sides,

[TABLE]

on $\Omega_{1}\cap\Omega_{2}$ (where we have assumed without a loss of generality that $c\geq 1$).

Controlling the effect of the penalty terms. Now we handle the appearances of the achieved penalty term $\hat{J}$. First, set $\lambda_{n}\geq(3/4)c^{2}r_{n}^{2}/C_{n}^{(2-w)/(2+w)}$, and denote

[TABLE]

Consider the case $\hat{J}\geq C_{n}$. Then $-1/C_{n}^{(2-w)/(2+w)}\geq-1/\hat{J}^{(2-w)/(2+w)}$, and

[TABLE]

thus, using the simple inequality in Lemma 10, we have $a\leq 4\lambda_{n}J^{*}$. In the case $\hat{J}<C_{n}$, we have by Lemma 10 again,

[TABLE]

Therefore, altogether, we conclude that $a\leq 6\lambda_{n}J^{*}$, and plugging this into (70) gives

[TABLE]

on $\Omega_{1}\cap\Omega_{2}$. The statement (29) as made in the theorem follows by dropping the nonnegative term $\|\hat{\Delta}\|_{n}^{2}/2$ from the left-hand side, and adjusting the constants $c,c_{1},c_{2},c_{3}>0$ as needed.

A.7 Proof of the best additive approximation bound in

(32)

We follow the exact same arguments as in the proof of Theorem 1, up until the last part, in which we control the achieved penalty terms $\hat{J}$. Now we deviate from the previous arguments, slightly. Set $\lambda_{n}\geq(3/2)c^{2}r_{n}^{2}/C_{n}^{(2-w)/(2+w)}$, and denote

[TABLE]

By the same logic as in the proof of Theorem 1, we have $a\leq 3\lambda_{n}J^{*}$. Plugging this into (70) gives

[TABLE]

on $\Omega_{1}\cap\Omega_{2}$. Rearranging,

[TABLE]

on $\Omega_{1}\cap\Omega_{2}$. But, setting $\tilde{f}=f^{\mathrm{best}}$, the bracketed term on the right-hand side above is nonpositive (by definition of $f^{\mathrm{best}}$ in (31)). This leads to (32), after adjusting $c,c_{1},c_{2},c_{3}>0$ as needed.

A.8 Preliminaries for the proof of Corollary 1

The following two lemmas will be helpful for the proof of Corollary 1.

Lemma 12**.**

Given $f=\sum_{j=1}^{d}f_{j}$, whose component functions are each $k$ times weakly differentiable, there exists an additive spline approximant $\check{f}=\sum_{j=1}^{d}\check{f}_{j}$, where $\check{f}_{j}\in\mathcal{G}_{j}$, the set of $k$th order splines with knots in the set $T_{j}$ defined in (27), for $j=1,\ldots,d$, such that

• (i)

$\mathrm{TV}(\check{f}_{j}^{(k)})\leq a_{k}\mathrm{TV}(f_{j}^{(k)})$, for $j=1,\ldots,d$; and

• (ii)

$\|\check{f}_{j}-f_{j}\|_{\infty}\leq a_{k}W_{\max}^{k}\mathrm{TV}(f_{j}^{(k)})$, for $j=1,\ldots,d$.

Above, $a_{k}\geq 1$ is a constant depending only on $k$, and we define $W_{\max}=\max_{j=1,\ldots,d}W_{j}$, where

[TABLE]

When the input points are drawn from a distribution $Q$ that satisfies Assumptions Assumption A1 (Continuous input distribution), Assumption A2 (Bounded input density), there are universal constants $c_{0},n_{0}>0$ such that for $n\geq n_{0}$, we have $W_{\max}\leq(c_{0}/b_{0})\log n/n$ with probability at least $1-2b_{0}d/n$, and so the bounds in (ii) become

[TABLE]

with probability at least $1-2b_{0}d/n$.

Proof.

Parts (i) and (ii) are simply a componentwise application of Proposition 7 in Mammen and van de Geer (1997). In particular, from their result, we know that for $j=1,\ldots,d$, there is a $k$th degree spline function $\check{f}_{j}$ whose knots lie in $T_{j}$ in (27), with $\mathrm{TV}(\check{f}_{j}^{(k)})\leq a_{k}\mathrm{TV}(f_{j}^{(k)})$ and

[TABLE]

where $a_{k}\geq 1$ depends only on $k$. (This result follows from strong quasi-interpolating properties of spline functions, from de Boor (1978).) This proves parts (i) and (ii).

When we consider random inputs drawn from a distribution $Q$ satisfying Assumptions Assumption A1 (Continuous input distribution), Assumption A2 (Bounded input density), the densities of the marginals $Q_{j}$, $j=1,\ldots,d$ will be bounded below by $b_{0}>0$, and thus there are universal constants $c_{0},n_{0}>0$ such that for $n\geq n_{0}$, we have $W_{j}\leq(c_{0}/b_{0})\log n/n$ with probability at least $1-2b_{0}/n$ (see, e.g., Lemma 5 in Wang et al. (2014)), for $j=1,\ldots,d$, and hence applying a union bound gives the result for $W_{\max}$. ∎

Lemma 13**.**

Given $f=\sum_{j=1}^{d}f_{j}$, whose component functions are each $k$ times weakly differentiable, there is an additive falling factorial approximant $\check{f}=\sum_{j=1}^{d}\check{f}_{j}$, where $\check{f}_{j}\in\mathcal{H}_{j}$, the set of $k$th order falling factorial functions defined over $X_{j}^{1},\ldots,X_{j}^{n}$, for each $j=1,\ldots,d$, such that

• (i)

$\mathrm{TV}(\check{f}_{j}^{(k)})\leq a_{k}\mathrm{TV}(f_{j}^{(k)})$, for $j=1,\ldots,d$; and

• (ii)

$\|\check{f}_{j}-f_{j}\|_{\infty}\leq a_{k}(W_{\max}^{k}+2k^{2}W_{\max})\mathrm{TV}(f_{j}^{(k)})$, for $j=1,\ldots,d$.

Again, $a_{k}\geq 1$ is a constant depending only on $k$, and $W_{\max}$ is as defined in Lemma 12. When the inputs are drawn from a distribution $Q$ satisfying Assumptions Assumption A1 (Continuous input distribution), Assumption A2 (Bounded input density), the bound in (ii) become

[TABLE]

with probability at least $1-2b_{0}d/n$.

Proof.

First we apply Lemma 12 to produce an additive spline approximant, call it $f^{*}=\sum_{j=1}^{d}f_{j}^{*}$, to the given $f=\sum_{j=1}^{d}f_{j}$. Next, we parametrize the spline component functions in a helpful way:

[TABLE]

where $\alpha_{j}^{1},\ldots,\alpha_{j}^{n}\in\mathbb{R}$ are coefficients and $g_{j1}\ldots,g_{jn}$ are the truncated power basis functions over the knot set $T_{j}$ defined in (27), and we write $g_{j\ell}(t)=t^{\ell-1}$, $\ell=1,\ldots,k$ without a loss of generality, for $j=1,\ldots,d$. It is not hard to check that $\mathrm{TV}((f_{j}^{*})^{(k)})=\sum_{\ell=k+2}^{n}|\alpha^{\ell}_{j}|$, for $j=1,\ldots,d$.

We now define $\check{f}=\sum_{j=1}^{d}\check{f}_{j}$, our falling factorial approximant, to have component functions

[TABLE]

where $h_{j1},\ldots,h_{jn}$ are the falling factorial basis functions defined over $X^{1}_{j},\ldots,X^{n}_{j}$, for $j=1,\ldots,d$. (Note that $\check{f}_{j}$ preserves the polynomial part of $f^{*}_{j}$ exactly, for $j=1,\ldots,d$.) Again, it is straightforward to check that $\mathrm{TV}(\check{f}_{j}^{(k)})=\sum_{\ell=k+2}^{n}|\alpha^{\ell}_{j}|$, for $j=1,\ldots,d$, i.e.,

[TABLE]

the inequality coming from part (i) of of Lemma 12. This verifies part (i) of the current lemma. As for part (ii), we note that Lemma 4 of Wang et al. (2014) shows that

[TABLE]

where recall $W_{j}$ is the maximum gap between sorted input points along the $j$th dimension, $j=1,\ldots,d$, as defined in Lemma 12. In fact, a straightforward modification of their proof can be used to strengthen this result to

[TABLE]

which means that by Holder’s inequality,

[TABLE]

Then, by the triangle inequality,

[TABLE]

where we have used part (ii) of Lemma 12. This verifies part (ii) of the current lemma.

Lastly, for random inputs drawn from a distribution $Q$ satisfying Assumptions Assumption A1 (Continuous input distribution), Assumption A2 (Bounded input density), the proof of (72) follows the same arguments as the proof of (71). ∎

A.9 Proof of Corollary 1

We consider first the statement in part (a). We must check that Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C3 (Entropy bound) hold for our choice of regularizer $J(g)=\mathrm{TV}(g^{(k)})$, and then we can apply Theorem 1. Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives) are immediate. As for Assumption Assumption C3 (Entropy bound), consider the univariate function class

[TABLE]

The results in Birman and Solomyak (1967) imply that for any set $Z_{n}=\{z^{1},\ldots,z^{n}\}\subseteq[0,1]$,

[TABLE]

for a universal constant $K>0$. As explained in Mammen (1991); Mammen and van de Geer (1997), this confirms that Assumption Assumption C3 (Entropy bound) holds for our choice of regularizer, with $w=1/(k+1)$. Applying Theorem 1, with $\tilde{f}=f_{0}$, gives the result in (33).

For the statement in part (b), note first that we can consider $k\geq 2$ without a loss of generality, as pointed out in Remark 6 following the corollary. Using Lemma 12, can choose an additive spline approximant $\check{f}$ to $f_{0}$, with components $\check{f}_{j}\in\mathcal{G}_{j}$, $j=1,\ldots,d$. Define $\tilde{f}_{j}$ to be the centered version of $\check{f}_{j}$, with zero empirical mean, $j=1,\ldots,d$. By the fact that centering does not change the penalty, and part (i) of the lemma, $\sum_{j=1}^{d}\mathrm{TV}(\tilde{f}_{j}^{(k)})\leq a_{k}\sum_{j=1}^{d}\mathrm{TV}(f_{0j}^{(k)})$. Also, using the fact that centering cannot increase the empirical norm, the triangle inequality, and (71), we get that with probability least $1-2b_{0}d/n$,

[TABLE]

When $\sum_{j=1}^{d}\mathrm{TV}(f_{0j}^{(k)})\leq C_{n}$, we see that $\|\tilde{f}-f_{0}\|_{n}$ is bounded by $C_{n}$ for large enough $n$. This meets required condition for Theorem 1, by the above display, the approximation error in (29) satisfies

[TABLE]

But when $n/(\log n)^{1+1/k}\geq n_{0}C_{n}^{(2k+2)/(2k^{2}+2k-1)}$, the right-hand side above is upper bounded by $a_{0}n^{-(2k+2)/(2k+3)}C_{n}^{2/(2k+3)}$, for a constant $a_{0}>0$. This establishes the result in (33) for restricted additive locally adaptive splines.

For the statement in part (c), we can again consider $k\geq 2$ without a loss of generality. Then the same arguments as given for part (b) apply here, but now we use Lemma 13 for the additive falling factorial approximant $\check{f}$ to $f_{0}$, and we require $n/(\log n)^{2k+3}\geq n_{0}C_{n}^{4k+4}$ for the approximation error to be bounded by the estimation error.

A.10 Preliminaries for the proof of Theorem 2

Our first lemma is similar to Lemma 9, but concerns univariate functions. As in Lemma 9, this result relies on Lemma 3.5 in van de Geer (1990) (see also van de Geer (2000)).

Lemma 14**.**

Let $\epsilon^{i}$, $i=1,\ldots,n$ be uniformly sub-Gaussian random variables having variance proxy $\sigma^{2}>0$. Let $J$ be a seminorm satisfying Assumption Assumption C3 (Entropy bound), and let $\rho>0$ be arbitrary. Then there exist constants $c_{1},c_{2},n_{0}>0$, depending only on $\sigma,K,w,\rho$, such that for all $c\geq c_{1}$ and $n\geq n_{0}$,

[TABLE]

with probabilty at least $1-\exp(-c_{2}c^{2})$, where we write $\|\cdot\|_{Z_{n}}$ for the empirical norm defined over a set of univariate points $Z_{n}=\{z^{1},\ldots,z^{n}\}\subseteq[0,1]$.

Proof.

Assume without a loss of generality that $\rho\geq 1$. Note that for any $Z_{n}=\{z^{1},\ldots,z^{n}\}\subseteq[0,1]$,

[TABLE]

by Assumption Assumption C3 (Entropy bound). As the right-hand side in the above entropy bound does not depend on $Z_{n}$, we can apply Lemma 3.5 in van de Geer (1990) to get the desired uniform control over all subsets. ∎

We give a coupling between the empirical and $L_{2}$ norms over $B_{J}(1)\cap B_{\infty}(\rho)$, using Theorem 14.1 in Wainwright (2017) (see also van de Geer (2000); Bartlett et al. (2005); Raskutti et al. (2012)).

Lemma 15**.**

Let $z^{i}$, $i=1,\ldots,n$ denote an i.i.d. sample from a distribution $\Lambda$ on $[0,1]$. Write $\|\cdot\|_{2}$ for the $L_{2}(\Lambda)$ norm, and $\|\cdot\|_{n}$ for the $L_{2}(\Lambda_{n})$ norm. Let $J$ satisfy Assumption Assumption C3 (Entropy bound), and let $\rho>0$ be arbitrary. Then there are constants $c_{1},c_{2},c_{3},n_{0}>0$, that depend only on $K,w,\rho$, such that for any $t\geq c_{1}n^{-1/(2+w)}$ and $n\geq n_{0}$,

[TABLE]

with probability at least $1-c_{2}\exp(-c_{3}nt^{2})$.

Proof.

Abbreviate $\mathcal{F}=B_{J}(1)\cap B_{\infty}(\rho)$. We will analyze the local Rademacher complexity

[TABLE]

the expectation being taken over i.i.d. draws $z^{i}$, $i=1,\ldots,n$ from $\Lambda$ and i.i.d. Rademacher variables $\sigma^{i}$, $i=1,\ldots,n$, as usual. Define the critical radius $\tau_{n}>0$ to the smallest solution of the equation

[TABLE]

We will prove $\tau_{n}\leq c_{1}n^{-1/(2+w)}$ for a constant $c_{1}>0$. Applying Theorem 14.1 in Wainwright (2017) would then give the result.

In what follows, we will use $c>0$ to denote a constant whole value may change from line to line (but does not depend on $z^{i}$, $i=1,\ldots,n$). Consider the empirical local Rademacher complexity

[TABLE]

As we are considering $t\geq\tau_{n}$, Corollary 2.2 of Bartlett et al. (2005) gives

[TABLE]

with probability at least $1-1/n$. Denote by $\mathcal{E}$ the event that this occurs. Then on $\mathcal{E}$,

[TABLE]

where in second line we used Dudley’s entropy integral (Dudley, 1967), and in the third line we used Assumption Assumption C3 (Entropy bound). On $\mathcal{E}^{c}$, note that we have the trivial bound $\mathcal{R}_{n}(\mathcal{F}\cap B_{2}(t))\leq\rho$. Therefore we can upper bound the local Rademacher complexity, splitting the expectation over $\mathcal{E}$ and $\mathcal{E}^{c}$,

[TABLE]

where the second inequality holds when $n$ is large enough, as we may assume $t\geq n^{-1/2}$ without a loss of generality. An upper bound on the critical radius $\tau_{n}$ is thus given by the solution of

[TABLE]

which is $t=cn^{-1/(2+w)}$. This completes the proof. ∎

We extend Lemma 5 to give a uniform sup norm bound on the functions in $B_{J}(1)\cap B_{2}(\rho)$.

Lemma 16**.**

Assume the conditions of Lemma 5. Then there are constants $R_{1},R_{2}>0$ that depend only on $k,b_{1},b_{2},L$, such that

[TABLE]

Proof.

For $m\in B_{J}(1)\cap B_{2}(\rho)$, decompose $m=p+g$ where $p$ is a polynomial of degree $k$, and $g$ is orthogonal to all polynomials of degree $k$ with respect to the $L_{2}(\Lambda)$ inner product. By Lemma 5, we have $\|g\|_{\infty}\leq R_{0}$ for a constant $R_{0}>0$, and by the triangle inequality,

[TABLE]

Now write

[TABLE]

where $\phi_{\ell}$, $\ell=1,\ldots,k+1$ are orthonormal polynomials on $[0,1]$ with respect to the $L_{2}(\Lambda)$ inner product. Then $\|\alpha\|_{2}=\|p\|_{2}\leq\rho+R_{0}$, from the second to last display, and $\|\alpha\|_{2}\leq\sqrt{k+1}(\rho+R_{0})$, so for any $x\in[0,1]$,

[TABLE]

where $c_{k}=\max_{\ell=1,\ldots,k+1}\|\phi_{\ell}\|_{\infty}$ is a constant that depends only on $k,b_{1}$ from Aptekarev et al. (2016). Therefore

[TABLE]

and defining $R_{1},R_{2}>0$ appropriately, this is of the desired form, and completes the proof. ∎

Our last two lemmas pertain to the function space

[TABLE]

We derive a one-sided bound on the $L_{2}$ norm in terms of the empirical norm, over $\mathcal{M}_{2}(1)$. Our proof uses Theorem 14.2 in Wainwright (2017), which is a somewhat unique theorem, because it does not require a global sup norm bound on the function class in consideration (unlike many standard results of this type).

Lemma 17**.**

Let $X^{i}$, $i=1,\ldots,n$ denote an i.i.d. sample from a distribution $Q$ on $[0,1]^{d}$ satisfying Assumption Assumption A3 (Product input distribution), and let $J$ satisfy Assumption Assumption C3 (Entropy bound). Then there are constants $c_{1},c_{2},c_{3},n_{0}>0$, that depend only on $b_{1},b_{2},k,L,K,w$, such that for any $c_{1}\sqrt{d}n^{-1/(2+w)}\leq t\leq 1$ and $n\geq n_{0}$,

[TABLE]

with probability at least $1-c_{2}\exp(-c_{3}nt^{2})$, where $\mathcal{M}_{2}(1)$ is the space in (73).

Proof.

Let $m\in\mathcal{M}_{2}(1)$ with $\|m\|_{2}\leq 1$. Then as $\|m\|_{2}^{2}=\sum_{j=1}^{d}\|m_{j}\|_{2}^{2}$, it follows that $\|m_{j}\|_{2}\leq 1$, $j=1,\ldots,d$, and by Lemma 16, we have $\|m_{j}\|_{\infty}\leq R_{1}+R_{2}$, $j=1,\ldots,d$. From the calculation in Example 14.6 of Wainwright (2017), we have the property

[TABLE]

where $C^{2}=(R_{1}+R_{2})^{2}+6$. Abbreviating $\mathcal{F}=\mathcal{M}_{2}(1)$, we will study the local Rademacher complexity

[TABLE]

and the associated critical radius $\tau_{n}>0$, defined as usual to be the smallest solution of

[TABLE]

We will establish $\tau_{n}\leq c_{1}\sqrt{d}n^{-1/(2+w)}$ for a constant $c_{1}>0$. Applying Theorem 14.2 in Wainwright (2017) would then give the result.

For the rest of the proof, we will use $c>0$ for a constant whose value may change from line to line; also, many statements will hold for large enough $n$, but this will not always be made explicit. Fix some $0<t\leq 1$. By $L_{2}$ orthogonality of the components of functions in $\mathcal{F}$,

[TABLE]

We now bound the inner supremum above, for an arbitrary $j=1,\ldots,d$. Denote by $\tau_{nj}$ the critical radius of $B_{J}(1)\cap B_{2}(|\beta_{j}|)$, denote $r_{n}=n^{-1/(2+w)}$, and define the abbreviation $a\vee b=\max\{a,b\}$. Observe

[TABLE]

The first three inequalities above hold with probability at least $1-1/3n^{2}$ each. The first inequality is by Theorem 3.6 in Wainwright (2017) (see also Example 3.9 in Wainwright (2017)); the second and third are by Lemma A.4 and Lemma 3.6 in Bartlett et al. (2005), respectively. The fourth upper bounds the local Rademacher complexity of $B_{J}(1)\cap B_{2}(|\beta_{j}|)$, and the fifth upper bounds the critical radius $\tau_{nj}$ of this class, both following the proof of Lemma 15 (recall, the functions in $B_{J}(1)\cap B_{2}(|\beta_{j}|)$ have a uniform sup norm bound of $\rho=R_{1}+R_{2}$, by Lemma 16). The last step also uses $\log n/n\leq r_{n}\sqrt{\log n/n}$ for $n$ sufficiently large. The final result of the above display holds with probability at least $1-1/n^{2}$; by a union bound, it holds with probability at least $1-d/n^{2}$ simultaneously over $j=1,\ldots,d$. Call this event $\mathcal{E}$. Then on $\mathcal{E}$,

[TABLE]

In the second line, we use Holder’s inequality $a^{T}b\leq\|a\|_{p}\|b\|_{q}$ for the first term, with $p=4/(2+w)$ and $q=4/(2-w)$; we use $a\vee b\leq a+b$ for the second term, along the bound $\|\beta\|_{1}\leq\sqrt{d}t$, and the fact that $r_{n}\sqrt{\log n/n}\leq r_{n}^{2}$ for large enough $n$.

Meanwhile, on $\mathcal{E}^{c}$, we can apply the simple bound $\|m\|_{\infty}\leq\sum_{j=1}^{d}\|m_{j}\|_{\infty}\leq\rho d$ for functions in $\mathcal{F}\cap B_{2}(t)$, where $\rho=R_{1}+R_{2}$ (owing to Lemma 16), and thus

[TABLE]

Splitting the expectation defining the local Rademacher complexity over $\mathcal{E},\mathcal{E}^{c}$, and using (74), (75),

[TABLE]

It can be easily verified that for $t=c\sqrt{d}r_{n}^{2}$, the upper bound in (76) is at most $t^{2}/C$. Therefore this is an upper bound on the critical radius of $\mathcal{F}$, which completes the proof. ∎

Lastly, we bound the gap in the empirical and $L_{2}$ means of functions in $\mathcal{M}_{2}(1)\cap B_{2}(t)$, for small enough $t$. The proof uses Theorem 2.1 in Bartlett et al. (2005).

Lemma 18**.**

Let $X^{i}$, $i=1,\ldots,n$ denote an i.i.d. sample from a distribution $Q$ on $[0,1]^{d}$ satisfying Assumption Assumption A3 (Product input distribution), and let $J$ satisfy Assumption Assumption C3 (Entropy bound). There are constants $c_{0},n_{0}>0$, that depend only on $b_{1},b_{2},k,L,K,w$, such that for any $0<t\leq 1$ and $n\geq n_{0}$,

[TABLE]

with probability at least $1-1/n$, where $\mathcal{M}_{2}(1)$ is the space in (73).

Proof.

This follows by combining the local Rademacher bound in (76) from the proof of Lemma 17 with Theorem 2.1 in Bartlett et al. (2005), and simplifying by keeping the dominant terms for large enough $n$. ∎

A.11 Proof of Theorem 2

At a high-level, the difference between this proof and that of Theorem 1 is that here we do not try to directly control the sub-Gaussian complexity term (as this would lead to a poor dependence on the dimension $d$). Instead, we reduce the problem to controlling univariate sub-Gaussian complexities, and then assemble the result using ties between the empirical and $L_{2}$ norms, and the decomposition property (35). We will use the same general notation as in the proof of Theorem 1: $c>0$ denotes a constant that will multiply our final bound, and will control the probability with which the final result holds; we will use the empirical norms and inner products with vector arguments, to be interpreted appropriately; we use the abbreviations $r_{n},\hat{\Delta}$, and so on. Finally, in many lines that follow, we will redefine $c$ by absorbing constant factors into it, without explicit notice.

The same arguments that led us to (68) yield the basic inequality

[TABLE]

where we write $\hat{\Delta}=\sum_{j=1}^{d}\hat{\Delta}_{j}$.

Bounding the sub-Gaussian complexity terms. We now bound the univariate sub-Gaussian complexity terms, appearing in the sum on the right-hand side in (77). For $j=1,\ldots,d$, define $g_{j}=\hat{\Delta}_{j}/(2\delta+\|\hat{\Delta}_{j}\|_{n})$, and note that by construction $J(g_{j})\leq 1$ and $\|g_{j}\|_{n}\leq 1$. By Lemma 7, there are constants $c_{0},R>0$ such that $\|g_{j}\|_{\infty}\leq R$ on an event whose probability is at least $1-\exp(-c_{0}n)$. Thus by Lemma 14, there are constants $c_{1},c_{2}>0$ such that for all $c\geq c_{1}$,

[TABLE]

on an event $\Omega_{1}$ with probability at least $1-\exp(c_{0}n)-\exp(-c_{2}c^{2})$. Plugging this into (77) gives

[TABLE]

on $\Omega_{1}$, where we used Lemma 10 in the second inequality.

Converting empirical norms into $L_{2}$ norms. For each $j=1,\ldots,d$, let $\bar{\Delta}_{j}=\langle\hat{\Delta}_{j},1\rangle_{2}$ be the $L_{2}$ mean of $\hat{\Delta}_{j}$, and $\tilde{\Delta}_{j}=\hat{\Delta}_{j}-\bar{\Delta}_{j}$ the $L_{2}$ centered version of $\hat{\Delta}_{j}$. Note that, for each $j=1,\ldots,d$, we have by empirical orthogonality $\|\tilde{\Delta}_{j}\|_{n}^{2}=\|\hat{\Delta}_{j}\|_{n}^{2}+|\bar{\Delta}_{j}|^{2}$, which implies $\|\hat{\Delta}_{j}\|_{n}\leq\|\tilde{\Delta}_{j}\|_{n}$. Applying this to upper bound the right-hand side in (78) gives

[TABLE]

on $\Omega_{1}$. We bound each empirical norm in the sum on the right-hand side in (79) by its $L_{2}$ norm counterpart. Now, for each $j=1,\ldots,d$, define $g_{j}=\tilde{\Delta}_{j}/(2\delta+\|\tilde{\Delta}_{j}\|_{2})$. Since $J(g_{j})\leq 1$ and $\|g_{j}\|_{2}\leq 1$, by Lemma 16, there is a constant $R>0$ such that $\|g_{j}\|_{\infty}\leq R$. We can hence apply Lemma 15 to the measure $\Lambda=Q_{j}$, which gives constants $c_{3},c_{4},c_{5}>0$ such that

[TABLE]

on an event $\Omega_{2}$ with probability at least $1-c_{4}d\exp(-c_{5}nr_{n}^{2})$, where recall $r_{n}=n^{-1/(2+w)}$, i.e.,

[TABLE]

on $\Omega_{2}$, where we assume $n$ is large enough so that $c_{3}r_{n}\leq\sqrt{3/2}$. Returning to (79), and using the simple inequality in Lemma 10, we have

[TABLE]

on $\Omega_{1}\cap\Omega_{2}$.

Invoking $L_{2}$ decomposability. We recall the key $L_{2}$ decomposition property (35), of additive functions with $L_{2}$ mean zero components. Using Holder’s inequality $a^{T}b\leq\|a\|_{p}\|b\|_{q}$ to bound the first sum on the right-hand side in (80), with $p=4/(2+w)$ and $q=4/(2-w)$, and Cauchy-Schwartz to bound the second sum in (80), we get

[TABLE]

on $\Omega_{1}\cap\Omega_{2}$, where we denote $\tilde{\Delta}=\sum_{j=1}^{d}\tilde{\Delta}_{j}$.

Converting back to empirical norm. We bound the $L_{2}$ norm of the centered error vector on the right-hand side in (81) with its empirical norm counterpart. By Lemma 17 applied to $m=\tilde{\Delta}/(2\delta)$, provided $n$ is large enough so that $c_{6}\sqrt{d}r_{n}\leq 1$, there are constants $c_{6},c_{7},c_{8}>0$ such that

[TABLE]

on an event $\Omega_{3}$ with probability at least $1-c_{7}\exp(-c_{8}dnr_{n}^{2})$. Plugging this into the right-hand side in (81), and using Lemma 10, we have

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}$. Using Lemma 11 on the first term above, with $a=d^{(2+w)/4}\delta^{w/2}/\sqrt{n}$ and $b=\|\tilde{\Delta}\|_{n}$, and simplifying, gives

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}$.

Deriving an empirical norm error bound. Note that in (83), we have $\|\hat{\Delta}\|_{n}$ on the left-hand side and $\|\tilde{\Delta}\|_{n}$ on the right-hand side, where $\tilde{\Delta}=\hat{\Delta}-\bar{\Delta}$ is the centered error vector, and we are abbreviating $\bar{\Delta}=\sum_{j=1}^{d}\bar{\Delta}_{j}$. We seek to bound $|\bar{\Delta}|$. Define $t=c_{6}\sqrt{d}r_{n}$, where $c_{6}$ is the constant in (82), and define

[TABLE]

Note that $J(m_{j})\leq J(\tilde{\Delta}_{j})/(2\delta)\leq 1$, for $j=1,\ldots,d$, by construction, and also

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}$, recalling (82). By Lemma 18 applied to $m$, provided $n$ is large enough such that $t=c_{6}\sqrt{d}r_{n}\leq 1$, there is a constant $c_{9}>0$ such that $|\langle m,1\rangle_{n}|\leq c_{9}t^{2}$ on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$, where $\Omega_{4}$ is an event with probability at least $1-1/n$, i.e.,

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$, i.e.,

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$. Thus, by empirical orthogonality,

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$, and assuming $n$ is large enough so that $2(\sqrt{2}c_{9}t)^{2}\leq 1/2$ and $2(2c_{9})^{2}t^{2}\delta\leq 1$, this becomes

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$. Using this on the right-hand side in (83) gives

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$. Using the simple inequality $2ab\leq a^{2}+b^{2}$ on the first term on the right-hand side above, with $a=c\sqrt{d}r_{n}\sqrt{\delta}$ and $b=\|\hat{\Delta}\|_{n}$, gives

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$. The empirical norm result in (36) in the theorem follows by dropping the nonnegative term $\|\hat{\Delta}\|_{n}^{2}/2$ from the left-hand side, and adjusting the constants $c,c_{1},c_{2},c_{3}>0$ as needed.

Deriving an $L_{2}$ norm error bound. Note that (85) also implies

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$. Recalling (84), this gives

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$. By $L_{2}$ orthogonality,

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$, where in the second line we used (82) and $|\bar{\Delta}|\leq\|\tilde{\Delta}\|_{n}$, and in the third line we used (86). Finally,

[TABLE]

on $\Omega_{1}\cap\Omega_{2}\cap\Omega_{3}\cap\Omega_{4}$. The $L_{2}$ norm result in (37) in the theorem follows by simply adjusting the constants $c,c_{1},c_{2},c_{3}>0$ as needed.

A.12 Proof of Corollary 2

The proof of the statement in part (a) is exactly as in the proof of part (a) in Corollary 1.

For part (b), we can consider $k\geq 2$ without a loss of generality, and start with an additive spline approximant $\check{f}$ to $f_{0}$ from Lemma 12. Let $\tilde{f}$ denote the result of centering each component of $\check{f}$ to have zero empirical mean. Then $\mathrm{TV}(\tilde{f}_{j}^{(k)})\leq a_{k}c_{n}=\delta$, $j=1,\ldots,d$, and just as in the proof of part (b) in Corollary 1, letting $\|\cdot\|$ denote either the empirical or $L_{2}$ norm, we have

[TABLE]

But when $n\geq n_{0}(dc_{n})^{(2k+3)/(2k+2)}$, the right-hand side above is bounded by $a_{0}dn^{-(2k+2)/(2k+3)}c_{n}$ for a constant $a_{0}>0$, which shows the approximation error terms in (36), (37) are of the desired order. This proves the desired result for restricted locally adaptive splines.

For part (c), we follow the same arguments, the only difference being that we construct a falling factorial approximant $\check{f}$ to $f_{0}$ from Lemma 13.

A.13 Preliminaries for the proof of Theorem 3

The next two results in this subsection are helper lemmas for the last lemma.

Lemma 19**.**

Let $J$ be a functional that satisfies Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C4 (Matching packing and covering number bounds). Then there are constants $\tilde{K}_{1},\tilde{\delta}_{1}>0$, that depend only on $k,L,K_{1},w$, such that for all $0<\delta\leq\tilde{\delta}_{1}$,

[TABLE]

where $\|\cdot\|_{2}$ is the $L_{2}(U)$ norm, with $U$ the uniform distribution on $[0,1]$, and $\Pi_{k}^{\perp}$ is defined by

[TABLE]

with $\mathcal{P}_{k}$ denoting the space of polynomials of degree $k$. In other words, $\Pi_{k}^{\perp}$ is the projection operator onto the space orthogonal (in $L_{2}(U)$) to the polynomials of degree $k$.

Proof.

Let $R_{0}>0$ be the constant from Lemma 5, when we take $\Lambda=U$. Note that

[TABLE]

In general, for $S=S_{1}+S_{2}$ and a norm $\|\cdot\|$, observe that, from basic relationships between covering and packing numbers,

[TABLE]

so that

[TABLE]

Applying this to our decomposition in (87),

[TABLE]

where in the second inequality we used Assumption Assumption C4 (Matching packing and covering number bounds) (assuming without a loss of generality that $R_{0}\geq 1$), and a well-known entropy bound for a finite-dimensional ball (e.g., Mammen (1991)), with $A>0$ being a constant that depends only on $R_{0}$. For small enough $\delta$, the right-hand side above is of the desired order, and this completes the proof. ∎

Lemma 20**.**

Let $d,M>0$ be integers, and $I=\{1,\ldots,M\}$. Denote by $H(u,v)=\sum_{j=1}^{d}1\{u_{j}\neq v_{j}\}$ the Hamming distance between $u,v\in I^{d}$. Then there is a subset $S\subseteq I^{d}$ with $|S|\geq(M/4)^{d/2}$ such that $H(u,v)\geq d/2$ for any $u,v\in S$.

Proof.

Let $\Omega_{0}=I^{d}$, $u_{0}=(1,\ldots,1)\in\Omega_{0}$. For $j=0,1,\ldots$, recursively define

[TABLE]

where $u_{j+1}$ is arbitrarily chosen from $\Omega_{j+1}.$ The procedure is stopped when $\Omega_{j+1}$ is empty; denote the last set defined in this procedure by $\Omega_{E}$, and denote $S=\{u_{0},\ldots,u_{E}\}$. For $0\leq i,j\leq E$, by construction, $H(u_{i},u_{j})>a$. For $j=0,\ldots,E$,

[TABLE]

The last step is true because we can choose $d-a$ positions in which $u$ matches $u_{j}$ in ${d\choose d-a}$ ways, and the rest of the $a$ positions can be filled arbitrarily in $M$ ways. Also note $M^{d}=n_{0}+\cdots+n_{E}$. Therefore

[TABLE]

which implies

[TABLE]

∎

The lemma below gives a key technical result used in the proof of Theorem 3.

Lemma 21**.**

Let $J$ be a functional that satisfies Assumptions Assumption C1 (Seminorm regularizer, null space of polynomials), Assumption C2 (Relative boundedness of derivatives), Assumption C4 (Matching packing and covering number bounds). Then there are constants $\bar{K}_{1},\bar{\delta}_{1}>0$, that depend only on $w,\tilde{K}_{1},\tilde{\delta}_{1}$, where $\tilde{K}_{1},\tilde{\delta}_{1}>0$ are the constants from Lemma 19, such that for all $0<\delta\leq\bar{\delta}_{1}$,

[TABLE]

where $\|\cdot\|_{2}$ is the $L_{2}(U)$ norm, with $U$ the uniform distribution on $[0,1]^{d}$, and $\Pi_{k,d}^{\perp}$ is defined by

[TABLE]

and $\mathcal{P}_{k,d}$ contains all functions of the form $p(x)=\sum_{j=1}^{d}p_{j}(x_{j})$, for polynomials $p_{j}$, $j=1,\ldots,d$ of degree $k$. In other words, $\Pi_{k,d}^{\perp}$ is the projection operator onto the space orthogonal (in $L_{2}(U)$) to the space $\mathcal{P}_{k,d}$ of additive polynomials of degree $k$.

Proof.

It is easy to check that the decomposability property of the $L_{2}(U)$ norm, in (35), implies a certain decomposability of the $L_{2}(U)$ projection operators $\Pi_{k,d},\Pi_{k,d}^{\perp}$ over additive functions:

[TABLE]

where $\Pi_{k},\Pi_{k}^{\perp}$ are projection operators onto $\mathcal{P}_{k}$ and its orthocomplement, respectively, as defined in Lemma 19. The decomposability result for $\Pi_{k,d}^{\perp}$ in particular implies that

[TABLE]

Abbreviate $M=M(\delta/\sqrt{d/2},\|\cdot\|_{2},\Pi_{k}^{\perp}(B_{J}(1)))$. By Lemma 19, we have for small enough $\delta$,

[TABLE]

Now let $g_{1},\ldots,g_{M}$ denote a $(\delta/\sqrt{d/2})$-packing of $\Pi_{k}^{\perp}(B_{J}(1))$. Let $I=\{1,\ldots,M\}$, and for $u\in I^{d}$, define $f_{u}\in\Pi_{k}^{\perp}(B_{J}^{d}(1))$ by

[TABLE]

i.e., $f_{u}$ is an additive function with components $g_{u_{j}}$, $j=1,\ldots,d$. If the Hamming distance between indices $u,v$ satisfies $H(u,v)\geq d/2$, then

[TABLE]

where we have again used the $L_{2}(U)$ decomposability property in (35). Thus, it is sufficient to find a subset $S$ of $I^{d}$ such that $u,v\in S\Rightarrow H(u,v)\geq d/2$. By Lemma 20, we can choose such an $S$ with $|S|\geq(M/4)^{d/2}$. For small enough $\delta$, such that $M\geq 16$, this gives the desired result because

[TABLE]

∎

A.14 Proof of Theorem 3

Clearly, by orthogonality, for any functions $\hat{f},f_{0}$,

[TABLE]

where $\Pi_{k,d},\Pi_{k,d}^{\perp}$ are projection operators onto $\mathcal{P}_{k,d}$ and its orthocomplement, respectively, defined in Lemma 19. Thus it suffices to consider the minimax error over $\Pi_{k,d}^{\perp}(B_{J}^{d}(c_{n}))$.

First, we lower bound the packing number and upper bound the covering number of the class $\Pi_{k,d}^{\perp}(B_{J}^{d}(c_{n}))$. The upper bound is more straightforward:

[TABLE]

The second inequality follows from property (88) in the proof of Lemma 21 and similar arguments to those in the proof of Lemma 8—except that we leverage the decomposability of the $L_{2}$ norm, as in (35), instead of using the triangle inequality. The third inequality follows from Assumption Assumption C4 (Matching packing and covering number bounds).

The lower bound is less straightforward, and is given by Lemma 21:

[TABLE]

We note that (90) holds for $0<\delta\leq\bar{\delta}_{1}$, where $\bar{\delta}_{1}>0$ is the constant from Lemma 21.

Now, following the strategy in Yang and Barron (1999), we use these bounds on the packing and covering numbers, along with Fano’s inequality, to establish the desired result. Let $f_{1},f_{2},\ldots,f_{M}$ be a $\delta_{n}$-packing of $\Pi_{k,d}^{\perp}(B_{J}(c_{n}))$, for $\delta_{n}>0$ to be specified later. Fix an arbitrary estimator $\hat{f}$, and let

[TABLE]

We will use $P_{X,f}$ and $\mathbb{E}_{X,f}$ to denote the probability and expectation operators, respectively, over i.i.d. draws $X^{i}\sim U$, $i=1,\ldots,n$ (where $U$ is the uniform distribution on $[0,1]^{d}$), and i.i.d. draws $Y^{i}|X^{i}\sim N(f(X^{i}),\sigma^{2})$, $i=1,\ldots,n$. Then

[TABLE]

where in the last inequality we have used the fact that if $\hat{Z}\neq j$, then $\hat{f}$ must be at least $\delta_{n}/2$ away from $f_{j}$, for each $j=1,\ldots,M$.

Abbreviate $q_{j}$ for the distribution $P_{f_{j}}$, $j=1,\ldots,M$, and define the mixture $\bar{q}=\frac{1}{M}\sum_{j=1}^{M}q_{j}$. By Fano’s inequality,

[TABLE]

where $\mathrm{KL}(P_{1}\,\|\,P_{2})$ denotes the Kullback-Leibler (KL) divergence between distributions $P_{1},P_{2}$. Let $g_{1},g_{2},\ldots,g_{N}$ be an $\epsilon_{n}$-covering of $\Pi_{k,d}^{\perp}(B_{J}^{d}(c_{n}))$, for $\epsilon_{n}>0$ to be determined shortly. Abbreviate $s_{\ell}$ for the distribution $P_{g_{\ell}}$, $\ell=1,\ldots,N$, and $\bar{s}=\frac{1}{N}\sum_{\ell=1}^{N}s_{\ell}$. Also, write $p(N(f(X),\sigma^{2}I))$ for the density of a $N(f(X),\sigma^{2}I)$ random variable, where $f(X)=(f(X^{1}),\ldots,f(X^{n}))\in\mathbb{R}^{n}$. Then

[TABLE]

In the first line above, we used the fact that $\sum_{j=1}^{M}\mathrm{KL}(q_{j}\,\|\,\bar{q})\leq\sum_{j=1}^{M}\mathrm{KL}(q_{j}\,\|\,s)$ for any other distribution $s$; in the second and third, we explicitly expressed and manipulated the definition of KL divergence; in the fourth, we used $\mathrm{KL}(q_{j}\,\|\,s_{\ell})=\|f_{j}(X)-g_{\ell}(X)\|_{2}^{2}/(2\sigma^{2})$, and for each $j$, there is at least one $\ell$ such that $\mathbb{E}_{X}\|f_{j}(X)-g_{\ell}(X)\|_{2}^{2}=\|f_{j}-g_{\ell}\|_{2}^{2}\leq\epsilon_{n}^{2}$; in the fifth line, we used the entropy bound from (89). Minimizing (93) over $\epsilon_{n}>0$ gives

[TABLE]

for a constant $\bar{K}_{2}>0$. Returning to Fano’s inequality (91), (92), we see that a lower bound on the minimax error is

[TABLE]

Therefore, a lower bound on the minimax error is $\delta_{n}^{2}/8$, for any $\delta_{n}>0$ such that

[TABLE]

and for large enough $n$, the first term on the right-hand side above will be larger than $2\log{2}$, so it suffices to have

[TABLE]

Set $\delta_{n}=(\bar{K}_{1}/4\bar{K}_{2})^{1/w}\sqrt{d}n^{-1/(2+w)}c_{n}^{w/(2+w)}$. Provided that $\delta_{n}\leq\bar{\delta}_{1}$, our log packing bound (90) is applicable, and ensures that (94) will be satisfied. This completes the proof.

A.15 Proof of Corollary 3

We only need to check Assumption Assumption C4 (Matching packing and covering number bounds) for $J(g)=\mathrm{TV}(g^{(k)})$, $w=1/(k+1)$, and then we can apply Theorem 3. As before, the entropy bound upper bound is implied by results in Birman and Solomyak (1967) (see Mammen (1991) for an explanation and discussion). The packing number lower bound is verified as follows. For $f$ a $(k+1)$ times weakly differentiable function on $[0,1]$,

[TABLE]

Hence

[TABLE]

Results in Kolmogorov and Tikhomirov (1959) imply that the space on the right-hand side satisfies the desired log packing number lower bound. This proves the result.

A.16 Proof of the linear smoother lower bound in

(44)

We may assume without a loss of generality that each $f_{0j}$, $j=1,\ldots,d$ has $L_{2}$ mean zero (since $f_{0}$ does). By the decomposability property of the $L_{2}$ norm over additive functions with $L_{2}$ mean zero components, as in (35), we have for any additive linear smoother $\hat{f}=\sum_{j=1}^{d}\hat{f}_{j}$,

[TABLE]

where $\bar{f}_{j}$ denotes the $L_{2}$ mean of $\hat{f}_{j}$, $j=1,\ldots,d$. Note that the estimator $\hat{f}_{j}-\bar{f}_{j}$ is itself a linear smoother, for each $j=1,\ldots,d$, since if we write $\hat{f}_{j}(x_{j})=w_{j}(x_{j})^{T}Y$ for a weight function $w_{j}$ over $x_{j}\in[0,1]$, then $\hat{f}_{j}(x_{j})-\bar{f}_{j}=\tilde{w}_{j}(x_{j})^{T}Y$ for a weight function $\tilde{w}_{j}(x_{j})=w_{j}(x_{j})-\int_{0}^{1}w_{j}(t)\,dt$. This, and the last display, imply that

[TABLE]

Now fix an arbitrary $j=1,\ldots,d$, and consider the $j$th term in the sum on the right-hand side above. Here we are looking at a linear smoother $\hat{f}_{j}$ fit to data

[TABLE]

which depends on the components $f_{0\ell}$, for $\ell\not=j$. This is why the supremum in the $j$th term of the sum on the right-hand side in (95) must be taken over $f_{0}\in\mathcal{F}_{k}^{d}(c_{n})$, rather than $f_{0j}\in\mathcal{F}_{k}(c_{n})$. Our notation $\hat{f}_{j}(Y)$ is used as a reminder to emphasize the dependence on the full data vector in (96).

A simple reformulation, by appropriate averaging over the lattice, helps untangle this supremum. Write $\hat{f}_{j}(x_{j})=w_{j}(x_{j})^{T}Y$ for a weight function $w_{j}$ over $x_{j}\in[0,1]$, and for each $v=1,\ldots,N$, let $I_{j}^{v}$ be the set of indices $i$ such that $X^{i}_{j}=v/N$. Also let

[TABLE]

and $\bar{Y}_{j}=(\bar{Y}^{1}_{j},\ldots,\bar{Y}^{N}_{j})\in\mathbb{R}^{N}$. Then note that we can also write $\hat{f}_{j}(x_{j})=\bar{w}_{j}(x_{j})^{T}\bar{Y}_{j}$ for a suitably defined weight function $\bar{w}_{j}$, i.e., note that we can think of $\tilde{f}_{j}$ as a linear smoother fit to data $\bar{Y}_{j}$, whose components follow the distribution

[TABLE]

where we let $\mu_{j}=\mu+\frac{1}{N}\sum_{\ell\not=j}\sum_{u=1}^{n}f_{0\ell}(u/N)$, and $\bar{\epsilon}^{v}_{j}$, $v=1,\ldots,n$ are i.i.d. $N(0,\sigma^{2}/N^{d-1})$. Recalling that $f_{0j}\in\mathcal{F}_{k}(c_{n})$, we are in a position to invoke univariate minimax results from Donoho and Johnstone (1998). As shown in Section 5.1 of Tibshirani (2014), the space $\mathcal{F}_{k}(c_{n})$ contains the Besov space $B_{1,1}^{k+1}(c_{n}^{\prime})$, for a radius $c_{n}^{\prime}$ that differs from $c_{n}$ only by a constant factor. Therefore, by Theorem 1 of Donoho and Johnstone (1998) on the minimax risk of linear smoothers fit to data from the model (97), we see that for large enough $N$ and a constant $c_{0}>0$,

[TABLE]

As we have reduced the lower bound to the minimax risk of linear smoothers over a Besov ball, we can see that the same result (98) indeed holds simultaneously over all $j=1,\ldots,d$. Combining this with (95) gives the desired result (44).

A.17 Proof of Theorem 4 and derivation

details for Algorithm 2

We show that the dual of (48) is equivalent to the additive trend filtering problem (12), and further, the Lagrange multipliers corresponding to the constraints $u_{0}=u_{j}$, for $j=1,\ldots,d$, are equivalent to the primal variables $\theta_{j}$, $j=1,\ldots,d$. Let $M=I-\mathbbm{1}\mathbbm{1}^{T}/n$, and rewrite problem (48) as

[TABLE]

We can write the Lagrangian of this problem as

[TABLE]

and we want to minimize this over $u_{0},\ldots,u_{d}$ to form the dual of (48). This gives

[TABLE]

We use the fact that the support function of $U_{j}$ is just $\ell_{1}$ penalty composed with $S_{j}D_{j}$ (invoking the duality between $\ell_{\infty}$ and $\ell_{1}$ norms),

[TABLE]

where recall we abbreviate $D_{j}=D^{(X_{j},k+1)}$, for $j=1,\ldots,d$, and this allows us to rewrite the above problem (99) as

[TABLE]

which is precisely the same as the original additive trend filtering problem (12).

This realization has important consequences. In the ADMM iterations (49), the scaled parameters $\rho\gamma_{j}$, $j=1,\ldots,d$ correspond to dual variables $\theta_{j}$, $j=1,\ldots,d$ in problem (48), which from the above calculation, are precisely primal variables in (12). Under weak conditions, ADMM is known to produce convergent dual iterates, e.g., Section 3.2 of Boyd et al. (2011) shows that if (i) the criterion is a sum of closed, convex functions and (ii) strong duality holds, then the dual iterates from ADMM converge to optimal dual solutions. (Convergence of primal iterates requires stronger assumptions.) Our problem (48) satisfies these two conditions, and so for the ADMM algorithm outlined in (49), the scaled iterates $\rho\gamma_{j}^{(t)}$, $j=1,\ldots,d$ converge to optimal solutions in the dual of (48), i.e., optimal solutions in the additive trend filtering problem (12). This proves the first part of the theorem.

As for the second part of the theorem, it remains to show that Algorithm 2 is equivalent to the ADMM iterations (47). This follows by notationally swapping $\gamma_{j}$, $j=1,\ldots,d$ for $\theta_{j}/\rho$, $j=1,\ldots,d$, rewriting the updates

[TABLE]

as

[TABLE]

using (45), and lastly, eliminating $u_{j}$, $j=1,\ldots,d$ from the $u_{0}$ update by solving for these variables in terms of terms of $\theta_{j}$, $j=1,\ldots,d$, i.e., by using

[TABLE]

A.18 Cyclic versus parallel backfitting

We compare the performances of the usual cyclic backfitting method in Algorithm 1 to the parallel version in Algorithm 2, on a simulated data set generated as in Section 5.1, except with $n=2000$ and $d=24$. We computed the additive trend filtering estimate (12) (of quadratic order), at a fixed value of $\lambda$ lying somewhere near the middle of the regularization path, by running the cyclic and parallel backfitting algorithms until each obtained a suboptimality of $10^{-8}$ in terms of the achieved criterion value (the optimal criterion value here was determined by running Algorithm 1 for a very large number of iterations). We used simply $\rho=1$ in Algorithm 2.

Figure 6 shows the progress of the two algorithms, plotting the suboptimality of the criterion value across the iterations. The two panels, left and right, differ in how iterations are counted for the parallel method. On the left, one full cycle of $d$ component updates is counted as one iteration for the parallel method—this corresponds to running the parallel algorithm in “naive” serial mode, where each component update is actually performed in sequence. On the right, $d$ full cycles of $d$ component updates is counted as one iteration for the parallel method—this corresponds to running the parallel algorithm in an “ideal” parallel mode with $d$ parallel processors. In both panels, one full cycle of $d$ component updates is counted as one iteration for the cyclic method. We see that, if parallelization is fully utilized, the parallel method cuts down the iteration cost by about a factor of 2, compared to the cyclic method. We should expect these computational gains to be even larger as the number of components $d$ grows.

A.19 Simulated homogeneously-smooth data

Figure 7 shows the results of a homogeneous simulation, as in Section 5.1 and Figure 4, except that for the base component trends we used sinusoids of equal (and spatially-constant) frequency:

[TABLE]

and we defined the component functions as $f_{0j}=a_{j}g_{0j}-b_{j}$, $j=1,\ldots,d$, where $a_{j},b_{j}$ were chosen to standardize $f_{0j}$ (give it zero empirical mean and unit empirical norm), for $j=1,\ldots,d$.