# Multivariate Trend Filtering for Lattice Data

Veeranjaneyulu Sadhanalaa Yu-Xiang Wangb Addison J. Huc Ryan J. Tibshiranic

(aGoogle bUniversity of California Santa Barbara cCarnegie Mellon University)

###### Abstract

We study a multivariate version of trend filtering, called Kronecker trend
filtering or KTF, for the case in which the design points form a lattice in
$d$ dimensions. KTF is a natural extension of univariate trend filtering
(Steidl et al., 2006; Kim et al., 2009; Tibshirani, 2014), and is defined by
minimizing a penalized least squares problem whose penalty term sums the
absolute (higher-order) differences of the parameter to be estimated along
each of the coordinate directions. The corresponding penalty operator can be
written in terms of Kronecker products of univariate trend filtering penalty
operators, hence the name Kronecker trend filtering. Equivalently, one can
view KTF in terms of an $\ell_{1}$-penalized basis regression problem where
the basis functions are tensor products of falling factorial functions, a
piecewise polynomial (discrete spline) basis that underlies univariate trend
filtering.

This paper is a unification and extension of the results in Sadhanala et al.
(2016, 2017). We develop a complete set of theoretical results that describe
the behavior of $k^{\text{th}}$ order Kronecker trend filtering in $d$
dimensions, for every $k\geq 0$ and $d\geq 1$. This reveals a number of
interesting phenomena, including the dominance of KTF over linear smoothers in
estimating heterogeneously smooth functions, and a phase transition at
$d=2(k+1)$, a boundary past which (on the high dimension-to-smoothness side)
linear smoothers fail to be consistent entirely. We also leverage recent
results on discrete splines from Tibshirani (2020), in particular, discrete
spline interpolation results that enable us to extend the KTF estimate to any
off-lattice location in constant-time (independent of the size of the lattice
$n$).

##  1 Introduction

We consider a standard nonparametric regression model, relating real-valued
responses $y_{i}\in\mathbb{R}$, $i=1,\ldots,n$ to design points
$x_{i}\in\mathcal{X}\subseteq\mathbb{R}^{d}$, $i=1,\ldots,n$,

$y_{i}=f_{0}(x_{i})+\epsilon_{i},\quad i=1,\ldots,n,$ (1)

where $f_{0}:\mathcal{X}\to\mathbb{R}$ is the (unknown) regression function to
be estimated, and $\epsilon_{i}$, $i=1,\ldots,n$ are mean zero stochastic
errors. In this paper, we will focus on functions $f_{0}$ that display
heterogeneous smoothness across the domain $\mathcal{X}$, in a sense we will
make precise later. We will also focus on the case in which the design points
form a $d$-dimensional lattice: that is, we assume $n=N^{d}$, and

$\\{x_{1},\ldots,x_{n}\\}=\\{1/N,2/N,\ldots,1\\}^{d}:=Z_{n,d}.$ (2)

The assumption of a uniformly-spaced lattice, embedded in the unit cube
$[0,1]^{d}$, is used only for simplicity. Essentially all results (both
methods and theory) translate over to the case of a somewhat more general
lattice structure, a Cartesian product
$\\{z_{i1}\\}_{i=1}^{N_{1}}\times\\{z_{i2}\\}_{i=1}^{N_{2}}\times\cdots\times\\{z_{id}\\}_{i=1}^{N_{d}}$,
where $\prod_{j=1}^{d}N_{j}=n$, and the sets in this product are otherwise
arbitrary. We return to this point in Section 9.1.

This paper is a unification and extension of Sadhanala et al. (2016, 2017)
(more will be said about the relationship to these papers in Section 1.3). The
models of smoothness for $f_{0}$ that we will study are based on _total
variation_ (TV). For a univariate function $g:[a,b]\to\mathbb{R}$, recall that
its total variation is defined as

$\mathrm{TV}(g;[a,b])=\sup_{a<z_{1}<\cdots<z_{m+1}<b}\;\sum_{i=1}^{m}|g(z_{i})-g(z_{i+1})|.$

For a multivariate function $f:[0,1]^{d}\to\mathbb{R}$, we will consider
notions of smoothness that revolve around the following _discrete_ version of
multivariate total variation:

$\mathrm{TV}(f;Z_{n,d})=\sum_{j=1}^{d}\sum_{\begin{subarray}{c}x,z\in
Z_{n,d}\\\ z=x+e_{j}/N\end{subarray}}|f(x)-f(z)|.$

Here we use $e_{j}$ to denote the $j^{\text{th}}$ standard basis vector in
$\mathbb{R}^{d}$, and hence the inner sum is taken over all pairs of lattice
points $x,z\in Z_{n,d}$ that differ in the $j^{\text{th}}$ coordinate by $1/N$
(and match in all other coordinates). We will also consider higher-order
versions of discrete multivariate TV, which are based on higher-order
differences in the summands in the above display. We will connect our discrete
notions of TV smoothness to standard continuum notions of total variation in
Section 3.

Broadly speaking, there are many multivariate nonparametric regression methods
available. Many of the methods in common use are _linear smoothers_ :
estimators of the form $\hat{f}(x)=w(x)^{\mathsf{T}}y$, for a suitable weight
function $w:\mathcal{X}\to\mathbb{R}^{n}$ (which can depend on the design
$x_{1},\ldots,x_{n}$), where we use $y=(y_{1},\ldots,y_{n})\in\mathbb{R}^{n}$
for the response vector. Examples include kernel smoothing, thin-plate
splines, and reproducing kernel Hilbert space estimators. A critical
shortcoming of linear smoothers is that they cannot be _locally adaptive_
—they cannot adapt to different local levels of smoothness exhibited by
$f_{0}$ over $\mathcal{X}$. This is a phenomenon that has been well-documented
in various settings; see Section 1.3.

The limitations of linear smoothers—and the need for nonlinear adaptive
methods—is also a central theme in this paper. The central method that drives
this story is a multivariate extension of trend filtering, which is indeed
nonlinear and locally adaptive, a claim that will be supported by experiments
and theory in the coming sections. We must note at the outset that all of the
developments in this paper hinge on the assumption of lattice data (meaning, a
lattice structure for the design points). A multivariate extension of trend
filtering for scattered data would require a completely different approach
(unlike, say, kernel smoothing or reproducing kernel Hilbert space methods,
which apply regardless of the structure of the design points). The
intersection of multivariate nonparametric regression methods and locally
adaptive methods is actually quite small, especially when we further intersect
this with the set of simple methods that are easy to use in practice, are
well-understood theoretically. For this reason, we see the contributions of
the current paper, though limited to lattice data, as being worthwhile. The
development of new multivariate trend filtering methods for scattered data is
important, and an extension is discussed in Section 9.2, but a comprehensive
study is left to future work.

###  1.1 Review: trend filtering

Before describing the main proposal, we review _trend filtering_ , a
relatively recent method for univariate nonparametric regression,
independently proposed by Steidl et al. (2006); Kim et al. (2009). For a
univariate design, equally-spaced (say) on the unit interval, $x_{i}=i/n$,
$i=1,\ldots,n$, and an integer $k\geq 0$, the $k^{\text{th}}$ order trend
filtering estimate is defined by the solution of the optimization problem:

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;\frac{1}{2}\|y-\theta\|_{2}^{2}+\lambda\|D^{(k+1)}_{n}\theta\|_{1}.$
(3)

Here $\lambda\geq 0$ denotes a tuning parameter,
$y=(y_{1},\ldots,y_{n})\in\mathbb{R}^{n}$ is the response vector, and
$D^{(k+1)}_{n}\in\mathbb{R}^{(n-k-1)\times n}$ is the difference operator of
order $k+1$, which we will also loosely call discrete derivative operator of
order $k+1$. This can be defined recursively in the following manner:

$\begin{gathered}D_{n}=\left[\begin{array}[]{rrrrrr}-1&1&0&\ldots&0&0\\\
0&-1&1&\ldots&0&0\\\ \vdots&&&&&\\\ 0&0&0&\ldots&-1&1\end{array}\right],\\\
D^{(k+1)}_{n}=D_{n-k}D^{(k)}_{n},\quad k=1,2,3\ldots.\end{gathered}$ (4)

The intuition behind problem (3) is as follows. The penalty term, which
penalizes the discrete derivatives of $\theta$ of order $k+1$, can be
equivalently seen as penalizing the _differences_ in $k^{\text{th}}$ discrete
derivatives of $\theta$ at adjacent design points (due to (4)). By the
sparsity-inducing property of the $\ell_{1}$ norm, the $k^{\text{th}}$
discrete derivatives of the trend filtering solution $\hat{\theta}$ will be
exactly equal at a subset of adjacent design points, and $\hat{\theta}$ will
therefore exhibit the structure of a $k^{\text{th}}$ degree piecewise
polynomial, with adaptively-chosen knots.

Here is a summary of the properties of trend filtering, as a nonparametric
regression tool.111We have defined it in this subsection for an evenly-spaced
design, for simplicity, but trend filtering can still be defined for an
arbitrary design, and all of the following properties still hold; see Section
9.1.

  * •

The discrete trend filtering estimate, which is defined over the design
points, can be “naturally” extended to a $k^{\text{th}}$ degree piecewise
polynomial function, in fact, a $k^{\text{th}}$ degree discrete spline, on
$[0,1]$.

  * •

Trend filtering is computationally efficient (several fast algorithms exist
for the structured, convex problem (3)), and is not much slower to compute
than (say) the smoothing spline.

  * •

Trend filtering is more locally adaptive than the smoothing spline (or any
linear smoother). This not only carries theoretical backing (next point), but
is clearly noticeable in practice as well.

  * •

Trend filtering attains the minimax rate (in squared empirical norm) of
$n^{-(2k+2)/(2k+3)}$ for estimating a function $f_{0}$ whose $k^{\text{th}}$
weak derivative is of bounded variation. The minimax linear rate (best worst-
case risk attained by a linear smoother) over this class is
$n^{-(2k+1)/(2k+2)}$.

Support for the above facts can be found in Tibshirani (2014), and the
discrete spline (numerical analytic) perspective behind trend filtering is
further developed in in Tibshirani (2020). More will be said about all of
these properties in the coming sections, as analogous properties will be
developed for a multivariate extension of trend filtering.

To prepare for this multivariate extension, it helps to recast the discrete
problem (3) in just a slightly different form. First, some notation: for a
vector $\theta\in\mathbb{R}^{n}$, we will (when convenient) index it by the
underlying design points, and write its components as $\theta(x_{i})$,
$i=1,\ldots,n$ in place of $\theta_{i}$, $i=1,\ldots,n$. Next, we define a
difference operator, which we will again loosely refer to as a discrete
derivative operator, by

$(\Delta\theta)(x_{i})=\begin{cases}\theta(x_{i+1})-\theta(x_{i})&\text{if
$i\leq n-1$},\\\ 0&\text{else}.\end{cases}$

Naturally, we can view $\Delta\theta$ as a vector in $\mathbb{R}^{n}$ with
components $(\Delta\theta)(x_{i})$, $i=1,\ldots,n$. Higher-order discrete
derivatives are obtained by repeated application of the same formula; we
abbreviate $(\Delta^{2}\theta)(x_{i})=(\Delta(\Delta\theta))(x_{i})$, and so
on. In this new notation, we can now rewrite problem (3) as

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;\frac{1}{2}\sum_{i=1}^{n}\big{(}y_{i}-\theta(x_{i})\big{)}^{2}+\lambda\sum_{i=1}^{n}|(\Delta^{k+1}\theta)(x_{i})|.$
(5)

###  1.2 Kronecker trend filtering

Let us return to the setting of multivariate design points on a lattice, as in
(2) (where recall we use $N=n^{1/d}$, assumed to be integral). Building from
the univariate notation and definitions at the end of the last subsection, we
now introduce multivariate analogs. For a vector $\theta\in\mathbb{R}^{n}$, we
will (when convenient) index its components by their lattice positions,
denoted $\theta(x)$, $x\in Z_{n,d}$. For each $j=1,\ldots,d$, we define the
discrete derivative of $\theta$ in the $j^{\text{th}}$ coordinate direction at
a location $x$ by

$(\Delta_{x_{j}}\theta)(x)=\begin{cases}\theta(x+e_{j}/N)-\theta(x)&\text{if
$x,x+e_{j}/N\in Z_{n,d}$},\\\ 0&\text{else}.\end{cases}$

We write $\Delta_{x_{j}}\theta\in\mathbb{R}^{n}$ for the vector with
components $(\Delta_{x_{j}}\theta)(x)$, $x\in Z_{n,d}$. As before, higher-
order discrete derivatives are simply defined by repeated application of the
above definition; we use abbreviations
$(\Delta_{x_{j}^{2}}\theta)(x)=(\Delta_{x_{j}}(\Delta_{x_{j}}\theta))(x)$,
$(\Delta_{x_{j},x_{\ell}}\theta)(x)=(\Delta_{x_{j}}(\Delta_{x_{\ell}}\theta))(x)$,
and so on.

With this notation in place, we define a multivariate version of trend
filtering, that we call _Kronecker trend filtering_ (KTF). Given an integer
$k\geq 0$, the $k^{\text{th}}$ order KTF estimate is defined by the solution
of the optimization problem:

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;\frac{1}{2}\sum_{i=1}^{n}\big{(}y_{i}-\theta(x_{i})\big{)}^{2}+\lambda\sum_{j=1}^{d}\sum_{x\in
Z_{n,d}}|(\Delta_{x_{j}^{k+1}}\theta)(x)|.$ (6)

Note the close analogy between (5) and (6): the latter extends the former by
adding up absolute discrete derivatives of $\theta$ of order $k+1$ along each
one of the $d$ coordinate directions. A similar intuition carries over from
the univariate case, regarding the role of the penalty in (6), and the
structure of the solution. As we can see, the KTF problem penalizes the
differences in $k^{\text{th}}$ discrete derivatives of $\theta$ at lattice
positions $x$ and $z$, for all $x$ and $z$ that are adjacent along any one of
the $d$ coordinate directions. By the sparsifying nature of the $\ell_{1}$
norm, the KTF solution $\hat{\theta}$ will have equal $k^{\text{th}}$ discrete
derivatives between neighboring points on the lattice (and more so for larger
$\lambda$, generally speaking). Hence, along any line segment parallel to one
of the coordinate axes, the KTF solution $\hat{\theta}$ will have the
structure of a $k^{\text{th}}$ degree piecewise polynomial, with adaptively-
chosen knots. This intuition will be made rigorous in Section 2.3.

We can also rewrite the KTF problem (6) in a more compact form, so that it
resembles (3):

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;\frac{1}{2}\|y-\theta\|_{2}^{2}+\lambda\|D_{n,d}^{(k+1)}\theta\|_{1},$
(7)

where we define

$D_{n,d}^{(k+1)}=\left[\begin{array}[]{c}D^{(k+1)}_{N}\otimes
I_{N}\otimes\cdots\otimes I_{N}\\\ I_{N}\otimes
D^{(k+1)}_{N}\otimes\cdots\otimes I_{N}\\\ \vdots\\\ I_{N}\otimes
I_{N}\otimes\cdots\otimes D^{(k+1)}_{N}\end{array}\right].$ (8)

Here $D^{(k+1)}_{N}\in\mathbb{R}^{(N-k-1)\times N}$ is the discrete derivative
matrix from (4) (as would be used in $k^{\text{th}}$ order univariate trend
filtering on $N$ points); $I_{N}\in\mathbb{R}^{N\times N}$ denotes the
identity matrix; and $A\otimes B$ denotes the Kronecker product of matrices
$A,B$. Each block of rows in (8) is made up of a total of $d-1$ Kronecker
products (a total of $d$ matrices). The Kronecker product structure behind the
penalty matrix in (8) is what inspires the name Kronecker trend filtering. In
a similar vein, we will also refer to $\|D_{n,d}^{(k+1)}\theta\|_{1}$ as the
$k^{\text{th}}$ order _Kronecker total variation_ (KTV) of $\theta$. Note that
when $k=0$, the KTF problem (7) reduces to anisotropic TV denoising on the
$d$-dimensional lattice, and when $d=1$, it reduces to $k^{\text{th}}$ order
univariate trend filtering (3).

Before we delve any deeper into its properties, which will start in Section 2,
we give a simple example of KTF in Figure 1. The example portrays an
underlying function $f_{0}$ in $d=2$ dimensions that has two peaks in opposite
corners of the unit square $[0,1]^{2}$, one larger and one smaller, and is
otherwise very smooth. Based on noisy observations of $f_{0}$ over a 2d
lattice, the estimates from KTF of orders $k=0,1,2$ are able to capture this
behavior. Meanwhile, estimates from kernel smoothing are not, and they either
oversmooth the larger peak or undersmooth the valleys, depending on the choice
of the bandwidth. Advocates of kernel smoothing might say that this problem
can be solved by giving the kernel a locally-varying bandwidth (that is,
modeling the bandwidth itself as a function of $x\in[0,1]^{d}$). While this
can work in principle, for example, using Lepski’s method (see Section 1.3),
it can often be hard to implement in practice. More to the point: the
comparison in Figure 1 is not meant to portray the kernel smoothing framework
as being wholly untenable; rather, it is just meant to portray the differences
in adaptivity of KTF versus kernel smoothing when each method is allowed only
a single tuning parameter.

Figure 1: Top left: underlying regression function $f_{0}$ evaluated over a
square lattice with $n=40^{2}=1600$ points, and associated responses (formed
by adding noise) shown as black points. Top middle and top right: kernel
smoothing (with a spherical Gaussian kernel) fit to this data using large and
small bandwidth values, respectively. Bottom left, middle, and right:
Kronecker trend filtering estimates of orders $k=0,1,2$, respectively (recall,
KTF with $k=0$ reduces to anisotropic total variation denoising). We see that,
in order to capture the larger of the two peaks in $f_{0}$, kernel smoothing
must significantly undersmooth the other peak (and surrounding areas);
instead, with more regularization, it undersmooths throughout. The KTF
estimates are able to adapt to heterogeneity in the smoothness of $f_{0}$.
Also, each exhibits a distinct structure based on the polynomial order $k$.

###  1.3 Related work

Our paper builds off a line of work on trend filtering, which, recall,
enforces smoothness by penalizing the $\ell_{1}$ norm of discrete derivatives
of a given order (Steidl et al., 2006; Kim et al., 2009; Tibshirani, 2014;
Wang et al., 2014; Ramdas and Tibshirani, 2016). This can be seen as a
discrete analog of the _locally adaptive regression spline_ estimator, which
penalizes the TV of a given order derivative of a function (Koenker et al.,
1994; Mammen and van de Geer, 1997). For an overview of trend filtering, its
connection to classical theory on splines and divided differences, and related
results that bridge discrete and continuum representations (such as that
connecting trend filtering and locally adaptive regression splines), we refer
to Tibshirani (2020). Trend filtering was extended to general graphs by Wang
et al. (2016) (more on this below), and was studied in an additive model
setting by Sadhanala and Tibshirani (2019).

Kronecker trend filtering was proposed in Sadhanala et al. (2017). Minimax
theory for KTF has been developed in different special cases, distributed
among several papers. Minimax results for TV denoising on lattices can be
found in Hutter and Rigollet (2016) (upper bounds) and Sadhanala et al. (2016)
(lower bounds); with respect to KTF and KTV classes, this covers the case of
$k=0$ and all dimensions $d$. Meanwhile, Sadhanala et al. (2017) derived
minimax results (matching upper and lower bounds) for all $k$ and $d=2$. The
current paper completes the landscape, deriving minimax theory for all
smoothness orders $k$ and all dimensions $d$. The majority of this paper
(including the minimax theory) was completed in 2018 and can be found in the
Ph.D. thesis of the first author (Sadhanala, 2019).

#### Graph-based TV methods.

Our paper is complementary to the line of research on locally adaptive
nonparametric estimation over graphs, such as Gavish et al. (2010); Sharpnack
et al. (2013); Wang et al. (2016); Göbel et al. (2018); Padilla et al. (2018,
2020); Ye and Padilla (2021). The lattice structure that we consider in this
paper can be cast as a particular $d$-dimensional grid graph with $n$ nodes
(that is, with all side lengths equal to $N=n^{1/d}$). However, many of the
methods proposed in the aforementioned references seek to be far more general
and operate over arbitrary graph structures. This generality comes with
several challenges, from conceptual to theoretical.

For example, Wang et al. (2016) developed _graph trend filtering_ (GTF), which
estimates a signal that takes values over the nodes of an arbitrary graph by
penalizing the $\ell_{1}$ norm its graph derivatives, which are defined via an
iterated graph Laplacian operator. Applied to the $d$-dimensional grid graph,
and translated into the notation of our paper, the penalty term used in
$k^{\text{th}}$ order GTF for a signal $\theta$ at a point in the lattice
$x\in Z_{n,d}$ is:

$\begin{cases}\displaystyle\sum_{j_{1}=1}^{d}\Bigg{|}\sum_{j_{2},\ldots,j_{q}=1}^{d}\Big{(}\Delta_{x_{j_{1}},x_{j_{2}}^{2},\ldots,x_{j_{q}}^{2}}\theta\Big{)}(x)\Bigg{|}&\text{for
$k$ even, where $q=k/2$},\\\
\displaystyle\Bigg{|}\sum_{j_{1},\ldots,j_{q}=1}^{d}\Big{(}\Delta_{x_{j_{1}}^{2},x_{j_{2}}^{2},\ldots,x_{j_{q}}^{2}}\theta\Big{)}(x)\Bigg{|}&\text{for
$k$ odd, where $q=(k+1)/2$}.\end{cases}$

Compared to the analogous penalty term used in KTF (6), which is just
$\sum_{j=1}^{d}|(\Delta_{x_{j}^{k+1}}\theta)(x)|$, we see that the above is
much is harder to interpret. GTF clearly considers some form of mixed
derivatives (whereas KTF does not and is more anisotropic), but it is
generally conceptually unclear what kind of smoothness GTF is promoting.
Further, Sadhanala et al. (2017) show that using the GTF penalty operator to
define a smoothness class for the analysis of multivariate signals is
problematic, in the following sense: for any $k\geq 1$ and any dimension $d$,
there are $k^{\text{th}}$ order Holder smooth functions whose discretization
to the lattice is arbitrarily nonsmooth as measured by the $k^{\text{th}}$
order GTF penalty operator. This is due to issues in the way the GTF penalty
operator measures smoothness on the boundaries of the lattice.

#### Continuous-time multivariate TV methods.

There is a rich body of work in applied mathematics on the denoising of
signals or images by promoting total variation smoothness, beginning with the
seminal paper by Rudin et al. (1992), which gave rise to the so-called Rudin-
Osher-Fatemi (ROF) functional. This was then further developed by Rudin and
Osher (1994); Vogel and Oman (1996); Chambolle and Lions (1997); Chan et al.
(2000), and others. Penalized versions of the ROF functional have also been
combined with multiscale (Candès and Guo, 2002; Dong et al., 2011) and $L^{1}$
(Chan and Esedoglu, 2005) data-fidelity terms. Papers in this line of work
tend to be cast in continuous-time, which means that the estimand, estimator,
and typically even the data itself are each functions of one or more variables
on a continuous domain (such as $[0,1]^{d}$). A related line of work, inspired
by ROF, considers discretization in its numerical analysis, see, for example,
Chambolle (2004, 2005); Almansa et al. (2008).

More recently, del Álamo et al. (2021) studied estimation of a multivariate
function of bounded variation under the white noise model, in arbitrary
dimension $d$. This may be interpreted as a continuum analog of our setting,
albeit for $k=0$. They derive minimax rates under the $L^{p}$ norm, for any
$p\geq 1$. When $p=2$ (matching our study), their estimator obtains (up to log
factors), the minimax rate of $n^{-1/d}$ on the squared $L^{2}$ error scale,
for any $d\geq 2$, which agrees with the minimax rate in our discrete setting.
We note that our methods, while motivated from discrete principles, do in fact
bear rigorous connections to continuous-time formulations of multivariate
total variation; see Sections 2.3 and 3.

#### Other models for multivariate TV.

There has been a recent surge of work studying different generalizations of TV
and trend filtering penalties to multiple dimensions, for example, Bibaut and
van der Laan (2019); Fang et al. (2021); Ortelli and van de Geer (2021);
Chatterjee and Goswami (2021); Ki et al. (2021). Many of these works are based
on the notion of Hardy-Krause variation of a multivariate function, or the
related notion of Vitali variation. For a function to be smooth in the Hardy-
Krause sense, it must exihibit an order of smoothness $k$ that scales with the
dimension $d$. This leads to estimation error rates that are (nearly)
dimension-free. However, practically speaking, expecting the inherent
smoothness of the regression function to be on par with the ambient dimension
may not be a reasonable assumption in many applications. A distinctive feature
of our work is that the smoothness order $k$ (which translates into the max
degree of the fitted local polynomial) is a user-defined parameter, and is not
tied in any way to the dimension $d$.

#### Locally-adaptive kernel smoothing.

Lepski’s method, which originated in the seminal paper by Lepskii (1991), is a
procedure for selecting a local bandwidth in kernel smoothing; roughly
speaking, at each point $x$ in the domain, it chooses the largest bandwidth
(from a discrete set of possible values) such that the kernel estimate at $x$
is within a carefully-defined error tolerance to estimates at smaller
bandwidths. Since its introduction, many papers have studied and generalized
Lepski’s method, see, for example, Lepskii (1992, 1993); Lepski et al. (1997);
Lepski and Spokoiny (1997); Kerkyacharian et al. (2001, 2008); Goldenshluger
and Lepski (2008, 2009, 2011, 2013); Lepski (2015). Papers in this line of
work are focused on establishing novel theoretical guarantees—for example,
using kernel smoothing as a base estimator to achieve minimax rates over
heterogeneous function classes, and not on applied considerations—like the
computational difficulties or ambiguities (choice of constants) encountered in
practical implementation.

Most related to the current paper is Kerkyacharian et al. (2001, 2008), which
considered estimation in anisotropic Besov classes using Lepski’s method,
under the white noise model. These papers considered the “dense” and “sparse”
cases, which roughly correspond to the cases $s>1/2$ and $s\leq 1/2$ in our
theory, respectively (see Figure 2). This is revisited in Remark 12.

#### Tensor product and hyperbolic wavelets.

Wavelets have a rich history in signal processing, approximation theory, and
other disciplines; classic references include Daubechies (1992); Chui et al.
(1992a); Meyer and Roques (1993); Mallat (2009). Seminal work by Donoho and
Johnstone (1998), on minimax estimation over univariate Besov spaces using
wavelet-based estimators, contributed greatly to the popularity of wavelets in
statistics. The earliest work on multivariate wavelet approximation appears to
have been Meyer (1987, 1990) and Mallat (1989b, a), which considered
“separable” wavelet bases, formed from tensor products of univariate wavelet
basis functions within each resolution level. A second approach, well-studied
in approximation theory but less so in statistics, considers constructing
“truly multivariate” wavelet bases by generating multiresolution subspaces in
the ambient domain, from which a wavelet basis can be derived (Meyer, 1990;
Riemenschneider and Shen, 1992; Chui et al., 1992b; Lorentz and Madych, 1992;
DeVore and Lucier, 1992). This gives rise to nonseparable wavelet bases, but
in both cases, the support of the wavelet function has the same scale in each
coordinate direction. The desire to effectively represent functions with
different degrees of smoothness in different coordinate directions led to the
development of hyperbolic wavelets, whose basis functions are formed by taking
tensor products of univariate wavelet basis functions _across_ resolution
levels (Neumann and von Sachs, 1997; DeVore et al., 1998; Neumann, 2000). Of
particular relevance to our work is the latter paper, and connections will be
drawn in Remark 12.

###  1.4 Summary and outline

A summary of results in this paper and outline for this paper is given below.

  * •

In Section 2, we derive some basic properties of KTF, including an equivalent
continuous-time formulation for problem (7), which provides insights into the
local structure of KTF estimates.

  * •

In Section 3, we derive an expression for higher-order multivariate TV (in the
standard measure-theoretic sense) in terms of an integrating univariate TV on
line segments running parallel to the coordinate axes. We use this to motivate
the definition of KTV smoothness, the central notion of smoothness used in
this paper—it is simply a discretization of the aforementioned integral.

  * •

In Section 4, we introduce the smoothness classes of interest for our study of
minimax theory, and examine the relationships between them.

  * •

In Section 5, we derive a complete set of results on the minimax estimation
risk, as measured in the squared $\ell_{2}$ norm, over the set
$\mathcal{T}_{n,d}^{k}(C_{n})$ of vectors $\theta$ with $k^{\text{th}}$ order
KTV smoothness satisfying $\|D_{n,d}^{(k+1)}\theta\|_{1}\leq C_{n}$, for a
given sequence $C_{n}>0$. We prove that KTF is minimax rate optimal (up to log
factors) for any $k,d$, and derive lower bounds on the minimax linear risk
(that is, the best worst-case risk over all linear smoothers) which show that
linear estimators are suboptimal for any $k,d$. Interestingly, the minimax
rates reveal a phase transition at $2(k+1)=d$, and in the low smoothness-to-
dimension regime, linear smoothers fail to be consistent altogether. See
Figure 2 for a more detailed summary.

  * •

In Section 6, we describe and compare specialized convex optimization
algorithms that can be used to compute the KTF solution in (7).

  * •

In Section 7, we present an extremely efficient and simple algorithm for
interpolating the discrete KTF estimate $\hat{\theta}$ (the solution in (7)),
defined over the lattice, into a function $\hat{f}$ (the solution in (12)),
defined over the underlying continuum domain $[0,1]^{d}$. Remarkably, this
interpolation method, which builds off discrete spline interpolation results
from Tibshirani (2020), runs in _constant-time_ (independent of $n$).

  * •

In Section 8, we carry out empirical experiments that compare KTF and various
other estimators, and examine whether the empirical error rates match the
minimax theory derived in Section 5.

  * •

In Section 9, we conclude with a discussion, and cover some extensions and
directions for future work.

Figure 2: Summary of the minimax results developed in this paper. The central
object of our study is the set $\mathcal{T}_{n,d}^{k}(C_{n}^{*})$ of vectors
$\theta$ defined over the $d$-dimensional lattice $Z_{n,d}$, with
$k^{\text{th}}$ order KTV smoothness satisfying
$\|D_{n,d}^{(k+1)}\theta\|_{1}\leq C_{n}^{*}$, for a sequence $C_{n}^{*}>0$
obeying what we call the canonical scaling, to be made precise later. The
following two statements hold, generally (regardless of $k,d$):  1. KTF
achieves the minimax rate (up to log factors) over
$\mathcal{T}_{n,d}^{k}(C_{n}^{*})$; and 2. no linear smoother is able to
achieve the minimax rate over this class. However, the story is more
interesting, due to a phase transition occurring at $2(k+1)=d$. Defining a
notion of _effective smoothness_ by $s=(k+1)/d$, this can be explained as
follows. When $s>1/2$, the minimax rate has the more classical form
$n^{-2s/(2s+1)}$, matching the minimax rate for a $k^{\text{th}}$ Holder class
in dimension $d$ (or an $s^{\text{th}}$ order Holder class in the univariate
case). Indeed, the lower bound on the minimax rate that we derive is given by
embedding a Holder class into $\mathcal{T}_{n,d}^{k}(C_{n}^{*})$. Meanwhile,
the minimax linear risk (the best worst-case risk among linear smoothers)
scales as $n^{-(2s-1)/(2s)}$, which can be interpreted as the rate of KTF (or
any other minimax optimal method) for a problem with a half less degree of
effective smoothness. When $s\leq 1/2$, the minimax rate takes on the less
classical form $n^{-s}$, and the lower bound is obtained by embedding a
suitable $\ell_{1}$ ball into $\mathcal{T}_{n,d}^{k}(C_{n}^{*})$. Further, the
gap between the minimax linear and nonlinear rates is even more dramatic: the
minimax linear rate is constant, which means no linear smoother is even
consistent over $\mathcal{T}_{n,d}^{k}(C_{n}^{*})$ (in the sense of worst-case
risk). Finally, though not reflected in the figure, we note that when $s<1/2$
the KTV class and its embedded Holder class exhibit different minimax rates,
$n^{-s}$ versus $n^{-2s/(2s+1)}$, respectively. Whether KTF can adapt to the
latter (faster) Holder rate in the low smoothness-to-dimension regime,
$s<1/2$, is an open question.

##  2 Basic properties

In this section, we cover a number of basis properties that reflect the
structure and complexity of KTF estimates.

###  2.1 Unpenalized component

We start by examining the null space of the KTF penalty matrix in (8). A word
on notation here, and in general: when convenient, we will use $x_{j}$ to
denote the $j^{\text{th}}$ component of a vector $x$, which should not be
confused with our use of $x_{i}$ to denote the $i^{\text{th}}$ design point
(itself a $d$-dimensional vector). This is an unfortunate clash of notation,
however, the meaning should always be clear form the context. (Furthermore, we
will keep the use of indices $i,j$ for the two cases consistent throughout to
aid the interpretation—hence $x_{j}$ will always be univariate, the
$j^{\text{th}}$ component of a vector $x$, and $x_{i}$ will always be
$d$-dimensional, the $i^{\text{th}}$ design point.)

######  Proposition 1.

The null space of the KTF penalty matrix in (8) has dimension $(k+1)^{d}$.
Furthermore, it is spanned by a polynomial basis made up of elements

$p(x)=x_{1}^{a_{1}}x_{2}^{a_{2}}\cdots x_{d}^{a_{d}},\quad x\in Z_{n,d},$

for all $a_{1},\ldots,a_{d}\in\\{0,\ldots,k\\}$.

The proof is elementary and is deferred to Appendix A (which contains all
proofs in this paper). This proposition reveals that the KTF penalty matrix
has quite a rich null space, thus KTF lets a significant component of the
response vector $y$ “pass through” unpenalized. In contrast to univariate
trend filtering, which preserves univariate polynomials of degree $k$
(precisely, it preserves the projection of $y$ onto this subspace), KTF
preserves “much more” than multivariate polynomials of degree $k$: it
preserves multivariate polynomials of _max degree_ $k$.222By a multivariate
polynomial of degree $k$, we mean (adhering to the standard classification) a
sum of terms of the form $b\prod_{j=1}^{d}x_{j}^{a_{j}}$, where the sum of
degrees satisfies $\sum_{j=1}^{d}a_{j}\leq k$. By a multivariate polynomial of
max degree $k$, we mean the same, but where the degrees satisfy $a_{j}\leq k$,
$j=1,\ldots,d$. To give an example, when $k=1$ and $d=2$, the KTF
estimator—which we might be tempted to call “linear-order” KTF (to use an
analogous term as we do in univariate trend filtering)—preserves any
polynomial of the form $p(x)=ax_{1}+bx_{2}+cx_{1}x_{2}$. This is of course
_not_ a linear function, but a bilinear one (due to the cross-product term
$x_{1}x_{2}$).

###  2.2 Review: trend filtering in continuous-time

For univariate trend filtering (3), an equivalent continuous-time formulation
was derived in Tibshirani (2014):

$\mathop{\mathrm{minimize}}_{f\in\mathcal{H}^{k}_{n}}\;\frac{1}{2}\sum_{i=1}^{n}\big{(}y_{i}-f(x_{i})\big{)}^{2}+\lambda\,\mathrm{TV}(f^{(k)}).$
(9)

Problems (3), (9) are equivalent in the sense that their solutions
$\hat{\theta},\hat{f}$, respectively, satisfy
$\hat{\theta}_{i}=\hat{f}(x_{i})$, $i=1,\ldots,n$. In (9), we use $f^{(k)}$ to
denote the $k^{\text{th}}$ weak derivative of $f$, and
$\mathrm{TV}(\cdot)=\mathrm{TV}(\cdot;[a,b])$ to denote the total variation
operator defined with respect to any interval $[a,b]$ containing the design
points (henceforth, when convenient, we will drop the underlying domain from
our notation for univariate TV). The optimization in (9) is performed over all
functions $f$ that lie in a space $\mathcal{H}^{k}_{n}$, which is the span of
$h^{k}_{n,j}$, $j=1,\ldots,n$, the following $k^{\text{th}}$ degree piecewise
polynomials:

$\begin{gathered}h^{k}_{n,i}(x)=\frac{1}{(i-1)!}\prod_{j=1}^{i-1}(x-j/n),\quad
i=1,\ldots,k+1,\\\ h^{k}_{n,i}(x)=\frac{1}{k!}\prod_{j=i-k}^{i-1}(x-j/n)\cdot
1\\{x>(i-1)/n\\},\quad i=k+2,\ldots,n.\end{gathered}$ (10)

(Here, for convenience, we interpret the empty product to be equal to 1.) The
functions in (10) are called the $k^{\text{th}}$ degree _falling factorial
basis_. Observe that they depend on the $n$ underlying design points
$1/n,2/n,\ldots,1$, which we express through the first subscript $n$ in
$h^{k}_{n,i}$. Note also the similarity between the above basis and the
standard truncated power basis for splines; in fact, when $k=0$ or $k=1$, the
two bases are equal and $\mathcal{H}^{k}_{n}$ is just a space of splines with
knots at the design points. However, when $k\geq 2$, this is no longer
true—the falling factorial functions are $k^{\text{th}}$ degree piecewise
polynomials with (mildly) discontinuous derivatives of all orders
$1,\ldots,k-1$, and therefore they span a different space than that of
$k^{\text{th}}$ degree splines—they space the space of $k^{\text{th}}$ degree
_discrete splines_ , which are piecewise polynomials that have continuous
discrete derivatives (rather than derivatives) at their knot points. See
Tibshirani (2020).

To be clear, the original formulation (3) is more computationally convenient
(it is a structured convex problem for which several fast algorithms exist,
discussed in Section 1.3). But the variational formulation (9) is important
because it provides rigorous backing to the intuition that a $k^{\text{th}}$
order trend filtering estimate exhibits the structure of a $k^{\text{th}}$
degree piecewise polynomial, with adaptively-chosen knots (a feature of the
$\ell_{1}$ penalty in (3) or TV penalty in (9)); moreover, it shows how to
extend the trend filtering estimate from a discrete sequence, defined over the
design points, to a function on the continuum interval $[0,1]$.

###  2.3 Continuous-time

We now develop a similar continuous-time representation for KTF.

######  Proposition 2.

Let $h^{k}_{N,i}:[0,1]\to\mathbb{R}$, $i=1,\ldots,N$ denote the
$k^{\text{th}}$ degree falling factorial functions (10) with respect to design
points $1/N,2/N,\ldots,1$, and $\mathcal{H}^{k}_{n,d}$ denote the space
spanned by all $d$-wise tensor products of these functions. That is,
abbreviating $h_{i}=h^{k}_{N,i}$, $i=1,\ldots,N$, this is the space of all
functions $f:[0,1]^{d}\to\mathbb{R}$ of the form

$f(x)=\sum_{i_{1},\ldots,i_{d}=1}^{N}\alpha_{i_{1},\ldots,i_{d}}h_{i_{1}}(x_{1})h_{i_{2}}(x_{2})\cdots
h_{i_{d}}(x_{d}),\quad x\in[0,1]^{d},$ (11)

for coefficients $\alpha\in\mathbb{R}^{n}$ (whose components we denote as
$\alpha_{i_{1},\ldots,i_{d}}$, $i_{1},\ldots,i_{d}\in\\{1,\ldots,N\\}$). Then
the KTF estimator defined in (7) is equivalent to the optimization problem:

$\mathop{\mathrm{minimize}}_{f\in\mathcal{H}^{k}_{n,d}}\;\frac{1}{2}\sum_{i=1}^{n}\big{(}y_{i}-f(x_{i})\big{)}^{2}+\lambda\sum_{j=1}^{d}\sum_{x_{-j}}\mathrm{TV}\bigg{(}\frac{\partial^{k}f(\cdot,x_{-j})}{\partial
x_{j}^{k}}\bigg{)},$ (12)

where $f(\cdot,x_{-j})$ denotes $f$ as function of the $j^{\text{th}}$
coordinate with all other dimensions fixed at $x_{-j}$,
$(\partial^{k}/\partial x_{j}^{k})(\cdot)$ denotes the $k^{\text{th}}$ partial
weak derivative operator with respect to $x_{j}$, and inner sum in the second
term in (12) is interpreted as a sum over $Z_{m,d-1}$, where $m=N^{d-1}$ (that
is, a sum over the $(d-1)$-dimensional uniformly-spaced lattice with $N^{d-1}$
total points). The discrete (7) and continuous-time (12) problems are
equivalent in the sense that at their solutions $\hat{\theta},\hat{f}$,
respectively, we have: $\hat{\theta}_{i}=\hat{f}(x_{i})$, $i=1,\ldots,n$.

######  Remark 1.

Similar to the discrete-time terminology, we will refer to
$\sum_{j=1}^{d}\sum_{x_{-j}}\mathrm{TV}(\partial^{k}f(\cdot,x_{-j})/\partial
x_{j}^{k})$, the penalty functional in the continuous-time representation (12)
of the KTF problem, as the $k^{\text{th}}$ order _Kronecker total variation_
(KTV) of $f$. A key implication of Proposition 2 is that the solution
$\hat{f}$ in (12) not only interpolates the solution $\hat{\theta}$ in (7),
but it is exactly as smooth in continuous-time as $\hat{\theta}$ is in
discrete-time, as measured by KTV:

$\|D_{n,d}^{(k+1)}\hat{\theta}\|_{1}=\sum_{j=1}^{d}\sum_{x_{-j}}\mathrm{TV}\bigg{(}\frac{\partial^{k}\hat{f}(\cdot,x_{-j})}{\partial
x_{j}^{k}}\bigg{)}.$

How to form the interpolant $\hat{f}$ is discussed briefly in the remark after
the next, and covered in detail in Section 7.

######  Remark 2.

From (11), the basis underlying the continuous-time representation (12) of the
KTF optimization problem, we can see that a $k^{\text{th}}$ order KTF estimate
exhibits the structure of a tensor product of $k^{\text{th}}$ degree discrete
splines, with adaptively knots, chosen to promote higher-order TV smoothness
along the coordinate axes. In other words, locally, it exhibits the structure
of a multivariate polynomial of max degree $k$. When $k=1$ and $d=2$, for
example, the structure is locally of the form
$\hat{f}(x)=ax_{1}+bx_{2}+cx_{1}x_{2}$, which is a bilinear function and has
local curvature. Such curvature is somewhat visible in Figure 1 (bottom row,
middle panel), and it will be even more apparent later when we discuss
interpolation, see Figure 6.

######  Remark 3.

The proof of Proposition 2 reveals that (12) has an equivalent form,
transcribed here for convenience:

$\mathop{\mathrm{minimize}}_{\alpha\in\mathbb{R}^{n}}\;\frac{1}{2}\Big{\|}y-\Big{(}H^{(k+1)}_{N}\otimes\cdots\otimes
H^{(k+1)}_{N}\Big{)}\alpha\Big{\|}_{2}^{2}+\lambda
k!\left\|\left[\begin{array}[]{c}I^{0}_{N}\otimes
H^{(k+1)}_{N}\otimes\cdots\otimes H^{(k+1)}_{N}\\\ H^{(k+1)}_{N}\otimes
I^{0}_{N}\otimes\cdots\otimes H^{(k+1)}_{N}\\\ \vdots\\\ H^{(k+1)}_{N}\otimes
H^{(k+1)}_{N}\otimes\cdots\otimes
I^{0}_{N}\end{array}\right]\alpha\right\|_{1},$ (13)

where $H^{(k+1)}_{N}\in\mathbb{R}^{N\times N}$ is the falling factorial basis
matrix (with columns given by evaluations of the falling factorial functions
at the design points) and $I^{0}_{N}\in\mathbb{R}^{(N-k-1)\times N}$ denotes
the last $N-k-1$ rows of the identity $I_{N}$. Interestingly, the penalty in
(13) is not a pure $\ell_{1}$ penalty on the coefficients $\alpha$ (as it
would be in basis form in the univariate case) but an $\ell_{1}$ penalty on
aggregated (positive linear combinations of) coefficients.

The basis formulation in (13) gives us a natural recipe for how to extend a
KTF estimate from a discrete sequence, defined over the lattice points, to a
function on the hypercube $[0,1]^{d}$. This is simply:

$\hat{f}(x)=\sum_{i_{1},\ldots,i_{d}=1}^{N}\hat{\alpha}_{i_{1},\ldots,i_{d}}h_{i_{1}}(x_{1})h_{i_{2}}(x_{2})\cdots
h_{i_{d}}(x_{d}),\quad x\in[0,1]^{d},$ (14)

where $\hat{\alpha}$ is the solution in (13). Though it is not obvious from
the basis expansion in (14), it turns out that at any point $x\in[0,1]^{d}$,
we can form the prediction $\hat{f}(x)$ in constant-time, starting from the
fitted values $\hat{\theta}_{i}=\hat{f}(x_{i})$, $i=1,\ldots,n$. This
leverages recent advances in discrete spline interpolation from Tibshirani
(2020), and is covered in Section 7.

###  2.4 Degrees of freedom

Given data from a model (1), where the errors $\epsilon_{i}$, $i=1,\ldots,n$
are i.i.d. with mean zero and variance $\sigma^{2}$, recall that the _degrees
of freedom_ of an estimator $\hat{\theta}_{i}=\hat{f}(x_{i})$, $i=1,\ldots,n$
of the means $\theta_{0i}=f_{0}(x_{i})$, $i=1,\ldots,n$ is a quantitative
reflection of its complexity, defined as (Efron, 1986; Hastie and Tibshirani,
1990):

$\mathrm{df}(\hat{\theta})=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\mathrm{Cov}(y_{i},\hat{\theta}_{i}).$

When the errors are Gaussian, Tibshirani and Taylor (2011, 2012) derived an
expression for the degrees of freedom of any generalized lasso estimator,
based on Stein’s formula (Stein, 1981). This covers KTF in (7) as a special
case, and thus translates into the following result for our setting: if
$\epsilon_{i}\sim N(0,\sigma^{2}$), $i=1,\ldots,n$ in (1), and $\hat{\theta}$
denotes the solution in (7) with active set

$A=\mathrm{supp}\big{(}D_{n,d}^{(k+1)}\hat{\theta}\big{)}=\big{\\{}i:\big{[}D_{n,d}^{(k+1)}\theta\big{]}_{i}\not=0\big{\\}},$

then

$\mathrm{df}(\hat{\theta})=\mathbb{E}\Big{[}\mathrm{nullity}\Big{(}\big{[}D_{n,d}^{(k+1)}]_{-A}\Big{)}\Big{]},$
(15)

where $\mathrm{nullity}(M)$ denotes the nullity (dimension of the null space)
of a matrix $M$, and $M_{-S}$ denotes the submatrix of $M$ given by removing
all rows indexed by a set $S$. From the above, we of course have the natural
estimator of degrees of freedom

$\widehat{\mathrm{df}}(\hat{\theta})=\mathrm{nullity}\Big{(}\big{[}D_{n,d}^{(k+1)}]_{-A}\Big{)},$
(16)

which is unbiased for (15).

The expression in (16) is easy to interpret when $k=0$: in this case, it
reduces to the number of connected constant pieces in the KTF solution
$\hat{\theta}$, where connectivity is interpreted with respect to the
underlying $d$-dimensional grid graph. This follows from the fact that the
penalty matrix $D^{(1)}_{n,d}$ in this case is the edge incidence operator on
the grid graph, and any submatrix of this penalty matrix (defined by removing
a subset of rows) is itself the edge incidence operator with respect to a
subgraph of the original grid (induced by removing a subset of the edges).
This result and its interpretation was already given in Tibshirani and Taylor
(2011, 2012) in the context of TV denoising on a graph.

When $k\geq 1$, the unbiased estimator of degrees of freedom in (16) is not as
easy to interpret. This is because, at a high level, there is no longer a
clear link between the local structure exhibited by $\hat{\theta}$ and whether
or not a particular entry of $D_{n,d}^{(k+1)}\hat{\theta}$ is nonzero.
However, we show in Appendix D that it is possible to compute the right-hand
side in (16) with a simple, direct algorithm that runs in linear time (more
precisely, the algorithm requires $O(ndk)$ operations).

##  3 Interlude: total variation on lines

In this section, we take a continuum perspective, looking at total variation
defined over functions in $\mathbb{R}^{d}$, and connect it to the discrete
notion of total variation used in the previous sections used to define the KTF
estimator.

###  3.1 Measure-theoretic total variation

Let $U$ be an open, bounded subset of $\mathbb{R}^{d}$ and $L^{p}(U)$ denote
the space of real-valued functions on $U$ with finite $L^{p}$ norm,
$\int_{U}|f(x)|^{p}\,dx<\infty$. A function $f\in L^{1}(U)$ is said to be of
_bounded variation_ (BV) provided $\mathrm{TV}(f;U)<\infty$, where

$\mathrm{TV}(f;U)=\sup\bigg{\\{}\int_{U}f(x)\,\mathrm{div}\phi(x)\,dx:\phi\in
C^{\infty}_{c}(U;\mathbb{R}^{d}),\;\|\phi(x)\|_{\infty}\leq 1\;\text{for all
$x\in U$}\bigg{\\}}.$ (17)

Above, $C^{\infty}_{c}(U;\mathbb{R}^{d})$ denotes the space of infinitely
continuously differentiable functions from $U$ to $\mathbb{R}^{d}$ with
compact support, and $\mathrm{div}(\cdot)$ denotes the divergence operator,
$\mathrm{div}\phi=\sum_{j=1}^{d}\partial\phi_{j}/\partial x_{j}$. We call
$\mathrm{TV}(f;U)$ the _total variation_ of $f$; this is the standard measure-
theoretic definition used in modern analysis; see, for example, Chapter 5 of
Evans and Gariepy (2015). To be clear, it would be more precise to call our
definition in (17) the _aniostropic_ total variation of $f$ (due to the use of
the $\ell_{\infty}$ norm in the constraint in (17) on the test function
$\phi$), but we often drop the reference to the anisotropic qualifier for
simplicity.

To build intuition, we note that if $f\in W^{1,1}(U)$, that is, $f$ is in
$L^{1}(U)$ and it is weakly differentiable and its weak derivative $\nabla f$
is also in $L^{1}(U)$, then

$\mathrm{TV}(f;U)=\int_{U}\|\nabla f(x)\|_{1}\,dx.$ (18)

Writing $\mathrm{BV}(U)$ for the space of bounded variation functions, the
above shows that $W^{1,1}(U)\subseteq\mathrm{BV}(U)$. Importantly, this is a
strict inclusion, because, for example, the indicator function of a set that
has smooth boundary is of bounded variation, but it is not in $W^{1,1}(U)$ (it
is not weakly differentiable).

###  3.2 Univariate total variation revisited

In order to connect the discrete notions of TV that we use in this paper to
the standard measure-theoretic definition of TV defined in (17), we must first
refine our definition of univariate TV. For a function $g:[a,b]\to\mathbb{R}$,
we define its total variation as:

$\mathrm{TV}(g;[a,b])=\sup_{\begin{subarray}{c}a<z_{1}<\cdots<z_{m+1}<b\\\
z_{1},\ldots,z_{m+1}\in\mathrm{AC}(g)\end{subarray}}\;\sum_{i=1}^{m}|g(z_{i})-g(z_{i+1})|,$
(19)

where the supremum is only taken over the set $\mathrm{AC}(g)$ of points of
approximate continuity of $g$. Approximate continuity is a weak notion of
continuity that excludes, for example, point discontinuities; see Section
1.7.2 of Evans and Gariepy (2015). Observe that the definition in (19) differs
from that given in the introduction in that the latter does not require that
the supremum be taken over points of approximate continuity. Some authors,
including Evans and Gariepy (2015), differentiate these definitions by calling
the latter the _variation_ of $g$ and (19) the _essential variation_ of $g$.
An intuitive way of interpreting their connection is as follows: the essential
variation of $g$ is the infimum of the variation achievable by any function
$\tilde{g}$ that agrees with $g$ Lebesgue almost everywhere.

The reason the refinement in (19) is important, when using the measure-
theoretic definition in (17) as a basis for defining the BV space, is that BV
functions (as with $L^{p}$ functions and Sobolev functions) are only well-
defined up to a set of Lebesgue measure zero. That is, if $f$ and $\tilde{f}$
agree Lebesgue almost everywhere, then their TV as defined in (17) (as with
$L^{p}$ norms or Sobolev norms) must also agree. Therefore, it should be clear
that (19) is the proper univariate notion here, as otherwise redefining $g$ at
a point would change its univariate TV (without restricting the supremum to
points of approximate continuity). Lastly, and reassuringly, the multivariate
measure-theoretic definition in (17) reduces to the univariate definition in
(19) once we take $U=(a,b)$ (see Theorem 5.21 in Evans and Gariepy (2015)).

###  3.3 Total variation on lines

We now proceed in the opposite direction to the end of the last subsection:
instead of reducing the multivariate definition to the univariate case, we
will use the univariate definition of TV to approach the multivariate one.
Interestingly, as we will see next, it turns out that (19) can be used as a
building block for (17) for an open, bounded set $U\subseteq\mathbb{R}^{d}$.
In words, the next result says that the multivariate notion of TV on $U$ is
given by aggregating the univariate notion on all line segments parallel to
the coordinate axes, anchored at boundary points of $U$.

######  Theorem 1.

Let $U\subseteq\mathbb{R}^{d}$ be an open, bounded, convex set. Then for any
$f\in\mathrm{BV}(U)$,

$\mathrm{TV}(f;U)=\sum_{j=1}^{d}\int_{U_{-j}}\mathrm{TV}\big{(}f(\cdot,x_{-j});I_{x_{-j}}\big{)}\,dx_{-j},$
(20)

where for each $j=1,\ldots,d$, we define $U_{-j}=\\{x_{-j}:(x_{j},x_{-j})\in
U\;\,\text{for some $x_{j}$}\\}$, and $I_{x_{-j}}=[a_{x_{-j}},b_{x_{-j}}]$,
with

$\displaystyle a_{x_{-j}}$ $\displaystyle=\inf\\{x_{j}:(x_{j},x_{-j})\in
U\\},$ $\displaystyle b_{x_{-j}}$
$\displaystyle=\sup\\{x_{j}:(x_{j},x_{-j})\in U\\}.$

Recall $f(\cdot,x_{-j})$ denotes $f$ as function of the $j^{\text{th}}$
coordinate with all other dimensions fixed at $x_{-j}$. Lastly, the univariate
TV operator in the integrand in (20) is to be interpreted in the essential
variation sense, as in (19).

######  Remark 4.

The assumption of convexity of $U$ in Theorem 1 is used for simplicity, to
ensure that each coordinatewise slice of $U$—intersecting it with a line
segment parallel to the coordinate axis—is an interval. The proof trivially
extends to the case in which each slice is a finite union of intervals; more
complex structures could likely be handled via more complex arguments.

######  Remark 5.

The above result is inspired by Theorem 5.22 of Evans and Gariepy (2015). In
the proof, we mollify $f$ in order to invoke the representation in (18) for
the TV of a smooth function, and then we leverage the separability of the
$\ell_{1}$ norm (that is, we leverage the fact that the integrand in (18)
decomposes into a sum of absolute partial derivatives) in order to derive
(20). Curiously, an analogous result does not seem straightforward to derive
for the case of isotropic TV; this being defined by using an $\ell_{2}$ norm
constraint on the test function $\phi$ in (17) (and for functions in
$W^{1,2}(U)$, it would reduce to the integral of $\ell_{2}$ norm of the weak
derivative, instead of the $\ell_{1}$ norm as in (18)).

###  3.4 Connection to KTV smoothness

We connect the representation in (20) to the KTF penalty functional, which we
call KTV smoothness. First note that we can rewrite the definition of the
anisotropic TV of a function $f$ in (17) as

$\mathrm{TV}(f;U)=\sum_{j=1}^{d}\underbrace{\sup\bigg{\\{}\int_{U}f(x)\;\frac{\partial\phi(x)}{\partial
x_{j}}\,dx:\phi\in C^{\infty}_{c}(U),\;|\phi(x)|\leq 1\;\text{for all $x\in
U$}\bigg{\\}}}_{V_{j}(f;U)},$

where $C^{\infty}_{c}(U)$ is the space of infinitely continuously
differentiable real-valued functions on $U$. Note that $V_{j}(f;U)$, as
defined above, measures the variation of $f$ along the $j^{\text{th}}$
coordinate direction. Now consider the following definition of $k^{\text{th}}$
order multivariate TV, for an integer $k\geq 0$:

$\mathrm{TV}^{k}(f;U)=\sum_{j=1}^{d}V_{j}\bigg{(}\frac{\partial^{k}f}{\partial
x_{j}^{k}};U\bigg{)},$ (21)

where $\partial^{k}f/\partial x_{j}^{k}$ denotes the $k^{\text{th}}$ partial
weak derivative of $f$ with respect to $x_{j}$. The formulation in (21) is, in
a sense, among the many possible options for higher-order TV in the
multivariate setting, the “most” anisotropic. It only looks at the variation
in the partial derivatives along the coordinate directions with respect to
which they are defined (that is, the variation in the $j^{\text{th}}$ partial
derivative along the $j^{\text{th}}$ coordinate axis).

Thanks to the representation in Theorem 1, we can rewrite the definition of
$k^{\text{th}}$ order TV in (21) as:

$\mathrm{TV}^{k}(f;U)=\sum_{j=1}^{d}\int_{U_{-j}}\mathrm{TV}\bigg{(}\frac{\partial^{k}f(\cdot,x_{-j})}{\partial
x_{j}^{k}}\bigg{)}\,dx_{-j},$ (22)

where we have made the dependence on the domain $I_{x_{-j}}$ in the TV
operator in the integrand implicit. Finally, we are ready to draw the
connection to KTF. Observe that the penalty functional underlying KTF, the
second term in the criterion of its continuous-time formulation (12), is given
by taking the notion of $k^{\text{th}}$ order TV in (22) and approximating the
integral via discretization; that is, the integral over $U_{-j}$ is simply
replaced by a sum over an embedded lattice.

##  4 Smoothness classes

We present various discrete smoothness classes, then connect them to each
other and to traditional Holder smoothness classes defined in continuous-time,
to derive what we refer to as _canonical scalings_ for the radii of the
discrete classes. This paves the way for the minimax analysis in the next
section.

###  4.1 Discrete TV and Sobolev classes

First we define the $k^{\text{th}}$ order _Kronecker total variation_ (KTV)
class, for a radius $\rho>0$, by

$\mathcal{T}_{n,d}^{k}(\rho)=\big{\\{}\theta\in\mathbb{R}^{n}:\|D_{n,d}^{(k+1)}\theta\|_{1}\leq\rho\big{\\}}.$
(23)

It is a priori unclear what scaling for the radius $\rho$ in (23) makes for an
“interesting” smoothness class for theoretical analysis. This is discussed at
length in Sadhanala et al. (2016), and was one of the original motivations for
that paper. There, it is shown that taking $\rho$ to be a constant (as
$n\to\infty$) leads to seemingly very fast minimax rates, but in this regime
trivial estimators turn out to be rate optimal (such as the sample mean
estimator $\hat{\theta}_{i}=\bar{y}$, $i=1,\ldots,n$).

In order to begin reasoning about scalings for $\rho$ in (23), it helps to
define the order $k+1$ discrete $\ell_{2}$-Sobolev class:

$\mathcal{W}_{n,d}^{k+1}(\rho)=\big{\\{}\theta\in\mathbb{R}^{n}:\|D_{n,d}^{(k+1)}\theta\|_{2}\leq\rho\big{\\}}.$
(24)

Observe that $\mathcal{W}_{n,d}^{k+1}(\rho)$ only considers partial
derivatives of order $k+1$ aligned with one of the coordinate axes, rather
than considering all mixed derivatives of total order $k+1$, as we would in a
traditional Sobolev class. For simplicity, we drop reference to the $\ell_{2}$
prefix when referring to (24) henceforth.

By the inequality $\|v\|_{2}\leq\sqrt{p}\|\theta\|_{1}$ for vectors
$v\in\mathbb{R}^{p}$, and the fact that the number of rows of
$D_{n,d}^{(k+1)}$ can be upper bounded by $dn$, we have the following
embedding:

$\mathcal{W}_{n,d}^{k+1}(\rho)\subseteq\mathcal{T}_{n,d}^{k}\big{(}\sqrt{dn}\rho\big{)},\quad\text{for
any $\rho>0$}.$

This shows that any reasonable regime for analysis must have $\rho$ varying
with $n$ in (23), or in (24) (or both), because a constant radius in one class
would translate into a growing or diminishing radius in the other, by the
above display. However, it still leaves unspecified what precise scalings for
the radii in (23) and (24), would correspond to “interesting” classes,
comparable in some sense to choices of radii in analogous continuous-time TV
or Sobolev smoothness classes. We answer this question in the next subsection,
by introducing discrete and continuum Holder classes, and pursuing further
embeddings.

###  4.2 Discrete and continuum Holder classes

Now we recall the traditional definition for the $k^{\text{th}}$ order Holder
class of functions from $[0,1]^{d}$ to $\mathbb{R}$, of radius $L>0$:

$C^{k}(L;[0,1]^{d})=\Bigg{\\{}f:[0,1]^{d}\to\mathbb{R}:\text{$f$ is $k$ times
differentiable and for all integers $\alpha_{1},\ldots,\alpha_{d}\geq 0$},\\\
\text{with
$\alpha_{1}+\cdots+\alpha_{d}=k$},\;\,\bigg{|}\frac{\partial^{k}f(x)}{\partial
x_{1}^{\alpha_{1}}\cdots\partial
x_{d}^{\alpha_{d}}}-\frac{\partial^{k}f(z)}{\partial
x_{1}^{\alpha_{1}}\cdots\partial x_{d}^{\alpha_{d}}}\bigg{|}\leq
L\|x-z\|_{2},\;\,\text{for all $x,z\in[0,1]^{d}$}\Bigg{\\}}.$

We define a discretized version of this class by simply evaluating the
functions in $C^{k}(L;[0,1]^{d})$ on the lattice $Z_{n,d}$:

$\mathcal{C}_{n,d}^{k}(L)=\big{\\{}\theta\in\mathbb{R}^{n}:\text{there exists
some $f\in C^{k}(L;[0,1]^{d})$ such that $\theta(x)=f(x)$, $x\in
Z_{n,d}$}\big{\\}}.$ (25)

The next proposition derives embeddings for the discrete Holder class (25)
into the discrete Sobolev (24) and KTV (23) classes. It is a direct
consequence of Lemma A.6 in Sadhanala et al. (2017) (which is a classical
result of sorts that quantifies the error of the forward difference
approximation of the derivative of a Holder function).

######  Proposition 3 (Sadhanala et al. 2017).

The discrete classes in (23)–(25) satisfy, for any $L>0$,

$\mathcal{C}_{n,d}^{k}(L)\subseteq\mathcal{W}_{n,d}^{k+1}\big{(}c_{1}Ln^{\frac{1}{2}-\frac{k+1}{d}}\big{)}\subseteq\mathcal{T}_{n,d}^{k}\big{(}c_{2}Ln^{1-\frac{k+1}{d}}\big{)},$
(26)

where $c_{1},c_{2}>0$ are constants depending only on $k,d$.

Motivated by the last result, as in Sadhanala et al. (2017), we define the
canonical scalings for the discrete Sobolev and KTV classes as

$\displaystyle B_{n}^{*}$ $\displaystyle=n^{\frac{1}{2}-\frac{k+1}{d}},$ (27)
$\displaystyle C_{n}^{*}$ $\displaystyle=n^{1-\frac{k+1}{d}},$ (28)

so that
$\mathcal{C}_{n,d}^{k}(1)\subseteq\mathcal{W}_{n,d}^{k+1}(c_{1}B_{n}^{*})\subseteq\mathcal{T}_{n,d}^{k}(c_{2}C_{n}^{*})$,
for constants $c_{1},c_{2}>0$ that depend only on $k,d$. Thus, by analogy to
classical results on nonparametric estimation over Holder spaces, we should
expect the minimax rate over $\mathcal{W}_{n,d}^{k+1}(B_{n}^{*})$ (in the
squared $\ell_{2}$ norm) to be $n^{-2(k+1)/(2(k+1)+d)}$. This is indeed the
case, as we will show at the end of Section 5. The minimax rate over
$\mathcal{T}_{n,d}^{k}(C_{n}^{*})$, on the other hand, will turn out to be
more exotic, and is the focus of the majority of the next section.

##  5 Estimation theory

We derive a number of results on estimation theory over KTV classes. We begin
by deriving upper bounds on the error of the KTF estimator, and then study
lower bounds. Throughout, we assume the data model in (1) with
$\theta_{0i}=f_{0}(x_{i})$, $i=1,\ldots,n$ and i.i.d. normal errors, to be
precise:

$y_{i}\sim N(\theta_{0,i},\sigma^{2}),\quad\text{independently, for
$i=1,\ldots,n$}.$ (29)

To set some basic notation, based on estimators $\hat{\theta}$ of the mean
$\theta_{0}$ in (29), we define for a subset
$\mathcal{K}\subseteq\mathbb{R}^{n}$,

$R(\mathcal{K})=\inf_{\hat{\theta}}\;\sup_{\theta_{0}\in\mathcal{K}}\;\frac{1}{n}\mathbb{E}\|\hat{\theta}-\theta_{0}\|_{2}^{2},$

which is called the _minimax risk_ over $\mathcal{K}$. Also of interest will
be

$R_{L}(\mathcal{K})=\inf_{\hat{\theta}\;\text{linear}}\;\sup_{\theta_{0}\in\mathcal{K}}\;\frac{1}{n}\mathbb{E}\|\hat{\theta}-\theta_{0}\|_{2}^{2},$

called the _minimax linear risk_ over $\mathcal{K}$, the infimum being
restricted to linear estimators $\hat{\theta}$ (that is, of the form
$\hat{\theta}=Sy$ for a matrix $S\in\mathbb{R}^{n\times n}$). To finish our
discussion of notation, for deterministic sequences $a_{n},b_{n}$ we write
$a_{n}=O(b_{n})$ when $a_{n}/b_{n}$ is upper bounded for large enough $n$, we
write $a_{n}=\Omega(b_{n})$ when $a_{n}^{-1}=O(b_{n}^{-1})$, and $a_{n}\asymp
b_{n}$ when both $a_{n}=O(b_{n})$ and $a_{n}=\Omega(b_{n})$. For random
sequences $A_{n},B_{n}$, we write $A_{n}=O_{\mathbb{P}}(B_{n})$ when
$A_{n}/B_{n}$ is bounded in probability. In the theory that follows, all
asymptotics are for $n\to\infty$ with $k,d$ fixed.

As a side remark, although not the focus of the current paper, analogous
theory can be established for graph trend filtering on grids (see Section 1.3
for a discussion of its relation to KTF), which we defer to Appendix B.

###  5.1 Upper bounds on estimation risk

To derive upper bounds on the risk of KTF, we leverage the following simple
generalization of a key result from Wang et al. (2016). Here and henceforth,
for an integer $a\geq 1$, we abbreviate $[a]=\\{1,\ldots,a\\}$.

######  Theorem 2 (Wang et al. 2016).

Consider the generalized lasso estimator $\hat{\theta}$ with penalty matrix
$D\in\mathbb{R}^{r\times n}$, defined by the solution of

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;\frac{1}{2}\|y-\theta\|_{2}^{2}+\lambda\|D\theta\|_{1}$
(30)

Suppose that $D$ has rank $q$, and denote by $\xi_{1}\leq\cdots\leq\xi_{q}$
its nonzero singular values. Also let $u_{1},\ldots,u_{q}\in\mathbb{R}^{r}$ be
the corresponding left singular vectors. Assume that these vectors, except
possibly for those in a set $I\subseteq[q]$, are _incoherent_ , meaning that
for a constant $\mu\geq 1$,

$\|u_{i}\|_{\infty}\leq\mu/\sqrt{n},\quad i\in[q]\setminus I.$

Then under the data model (29), choosing

$\lambda\asymp\mu\sqrt{\frac{\log r}{n}\sum_{i\in[q]\setminus
I}\frac{1}{\xi_{i}^{2}}},$

the generalized lasso estimator satisfies

$\frac{1}{n}\|\hat{\theta}-\theta_{0}\|_{2}^{2}=O_{\mathbb{P}}\Bigg{(}\frac{\mathrm{nullity}(D)}{n}+\frac{|I|}{n}+\frac{\mu}{n}\sqrt{\frac{\log
r}{n}\sum_{i\in[q]\setminus
I}\frac{1}{\xi_{i}^{2}}}\cdot\|D\theta_{0}\|_{1}\,\Bigg{)}.$ (31)

We will now apply this result to KTF in (7), and choose the set $I$ in order
to balance the second and third terms on the right-hand side in (31).
Throughout, we will make reference to the _effective degree of smoothness_ (or
effective smoothness for short), defined by

$s=\frac{k+1}{d}.$

######  Theorem 3.

Let $\hat{\theta}$ denote the KTF estimator in (7). Under the data model (29),
denote $C_{n}=\|D_{n,d}^{(k+1)}\theta_{0}\|_{1}$, and assume $C_{n}>0$.
Choosing

$\displaystyle\lambda\asymp\begin{cases}\sqrt{\log n}&\text{if $s<1/2$},\\\
\log n\vphantom{n^{\frac{2s-1}{2s+1}}}&\text{if $s=1/2$},\\\ (\log
n)^{\frac{1}{2s+1}}(n/C_{n})^{\frac{2s-1}{2s+1}}&\text{if
$s>1/2$},\end{cases}$

the KTF estimator satisfies

$\frac{1}{n}\|\hat{\theta}-\theta_{0}\|_{2}^{2}=O_{\mathbb{P}}\bigg{(}\frac{1}{n}+\frac{\lambda}{n}C_{n}\bigg{)}.$

######  Remark 6.

The result in the above theorem for the case of $k=0$ (TV denoising) and any
$d$ was already established in Hutter and Rigollet (2016). The result for
$d=2$ and any $k$ was given in Sadhanala et al. (2017). For $d=1$ and any $k$,
the result was established in Mammen and van de Geer (1997); Tibshirani (2014)
(though the results in the latter papers are sharper by log factors).

Theorem 3 covers all $k$ and $d$. Its proof, an application of Theorem 2 to
KTF, requires checking the incoherence of the penalty matrix
$D=D_{n,d}^{(k+1)}$, and bounding the partial sum of squared reciprocal
singular values $\sum_{i\in[q]\setminus I}\,\xi_{i}^{-2}$, for an appropriate
set $I$ that corresponds to small eigenvalues. The former—checking the
incoherence of the left singular vectors of the KTF penalty matrix—turns out
to be the harder calculation. However, the hardest part of this calculation
was already done in Sadhanala et al. (2017), who established the incoherence
of $D^{(k+1)}_{n}$ (the trend filtering penalty matrix, or equivalently, the
KTF penalty matrix when $d=1$) by appealing to complex approximation arguments
for the eigenvectors of Toeplitz matrices. To handle KTF in arbitrary
dimension $d$, in the current paper, we use careful arguments that relate
singular vectors of Kronecker products to singular vectors of their
constituent matrices.

######  Remark 7.

It is interesting to note that the case $s\leq 1/2$ (that is, $2k+2\leq d$)
appears to be special, in that KTF ends up being adaptive to the underlying
level of smoothness $C_{n}$. In other words, the prescribed choice of tuning
parameter is $\lambda\asymp\sqrt{\log n}$ when $s<1/2$, and $\lambda\asymp\log
n$ when $s=1/2$, neither of which depend on $C_{n}$.

Under the canonical scaling $C_{n}\asymp C_{n}^{*}$ in (28), the error bound
in Theorem 3 reduces to the following.

######  Corollary 1.

Assume the conditions of Theorem 3. Under the canonical scaling in (28) for
the KTV smoothness level, $C_{n}\asymp C_{n}^{*}=n^{1-s}$, the KTF estimator
satisfies

$\displaystyle\frac{1}{n}\|\hat{\theta}-\theta_{0}\|_{2}^{2}=\begin{cases}O_{\mathbb{P}}\big{(}n^{-s}\sqrt{\log
n}\,\big{)}&\text{if $s<1/2$},\\\ O_{\mathbb{P}}(n^{-s}\log
n)\vphantom{n^{-\frac{2s}{2s+1}}}&\text{if $s=1/2$},\\\
O_{\mathbb{P}}\big{(}n^{-\frac{2s}{2s+1}}(\log
n)^{\frac{1}{2s+1}}\big{)}&\text{if $s>1/2$}.\end{cases}$ (32)

######  Remark 8.

Recall the continuous-time analog of the KTF penalty in (12). It turns out, as
we show in Appendix A.6, that the assumption that the mean vector $\theta_{0}$
is KTV smooth in Theorem 3 and Corollary 1 can be broadened into an assumption
that the KTV of the mean _function_ $f_{0}$ is KTV smooth, in the sense of the
penalty functional in (12).

###  5.2 Lower bounds on estimation risk

Next we give lower bounds on the minimax risk over KTV classes.

######  Theorem 4.

The minimax risk for KTV class defined in (23) satisfies, for any sequence
$C_{n}\leq n$,

$R\big{(}\mathcal{T}_{n,d}^{k}(C_{n})\big{)}=\Omega\Bigg{(}\frac{1}{n}+\frac{C_{n}}{n}+\bigg{(}\frac{C_{n}}{n}\bigg{)}^{\frac{2}{2s+1}}\Bigg{)}.$
(33)

######  Remark 9.

The result in Theorem 4 for $s\geq 1/2$ (that is, $2k+2\geq d$) was already
derived in Sadhanala et al. (2017). More precisely, these authors established
the third term on the right-hand side in (33), by using the Holder embedding
in (26) of Proposition 3, and suitably adapting classical results on minimax
bounds for Holder spaces (Korostelev and Tsybakov, 2003; Tsybakov, 2009),
which ends up being tight in the $s>1/2$ regime. Moreover, the lower bound in
the middle term in (33) was derived by Sadhanala et al. (2016), for $k=0$ and
all $d$, which was obtained by embedding an appropriate $\ell_{1}$ ball into
$\mathcal{T}_{n,d}^{k}(C_{n})$, and appealing to minimax theory over
$\ell_{1}$ balls from Birge and Massart (2001). This ends up being tight for
$k=0$ and all $d$.

In the current work, we establish tight lower bounds over all $k,d$, by
essentially combining these two strategies. As we will see comparing to upper
bounds in the next remark, the Holder embedding ends up being tight for
$s>1/2$, and the $\ell_{1}$ embedding for $s\leq 1/2$.

######  Remark 10.

Plugging in the canonical scaling $C_{n}\asymp C_{n}^{*}=n^{1-s}$ from (28),
and simplifying to dominant terms, we see that (33) becomes

$R\big{(}\mathcal{T}_{n,d}^{k}(C_{n}^{*})\big{)}=\begin{cases}\Omega(n^{-s})&\text{if
$s\leq 1/2$},\\\ \Omega(n^{-\frac{2s}{2s+1}})&\text{if $s>1/2$}.\end{cases}$

By comparing this to (32), we can see that KTF is minimax rate optimal for
estimation over KTV classes, under the canonical scaling, up to log factors.

######  Remark 11.

From the upper bound in (32) and the Holder embedding in (26), we see that for
$s\geq 1/2$, KTF achieves the rate of $n^{-2s/(2s+1)}$ (up to log factors)
over $\mathcal{H}_{d}^{k}(1)$. This matches the optimal rate for estimation
over Holder classes (see Sadhanala et al. (2017) for a formal statement and
proof for the discretized class $\mathcal{H}_{d}^{k}(1)$), which means that
KTF automatically adapts Holder smooth signals.

However, when $s<1/2$, the same results show that KTF achieves a rate of
$n^{-s}$ (ignoring log factors) over $\mathcal{H}_{d}^{k}(1)$, which is slower
than the optimal rate of $n^{-2s/(2s+1)}$. Whether this upper bound is
pessimistic and KTF can actually adapt to $\mathcal{H}_{d}^{k}(1)$, or whether
this upper bound is tight and KTF fails to adapt to Holder smooth signals,
remains to be formally resolved. Later, we investigate this empirically in
Section 8.3.

######  Remark 12.

It is interesting to compare minimax results for anisotropic Besov spaces,
under the white noise model. Past work on this topic includes Neumann (2000)
with a focus on hyperbolic wavelets, as well as Kerkyacharian et al. (2001,
2008) with a focus on Lepski’s method applied to kernel smoothing. A nice
summary of past work, and what appears to be the most comprehensive results,
can be found in Lepski (2015). It should be noted that this line of work
considers a more general setup than ours (albeit in the white noise model), in
a few ways: anisotropic Besov classes with an arbitrary smoothness index in
each coordinate direction; error measured in $L^{p}$ norm, for arbitrary
$p\geq 1$; and so on. That said, translating their results to match our
setting as best as we can (recalling general embeddings of BV spaces into
Besov spaces; for example, DeVore and Lorentz (1993)), we find that the
minimax rate under the squared $L^{2}$ loss, for the anisotropic Besov class
with integrability index 1, smoothness index $k+1$ in each coordinate
direction, and any third index $q\geq 1$, is indeed $n^{-2s/(2s+1)}$ for
$s>1/2$ (or $2k+2>d$). This regime is what Lepski and others refer to as the
“dense zone”.

In what they refer to as the “sparse zone”, $s\leq 1/2$, the minimax risk
under the $L^{2}$ loss for the same Besov class is a constant. This is due to
the fact that this Besov space fails to embed compactly into $L^{2}$ (see
Section 5.5 of Johnstone (2015) for a general statement of the implications of
this phenomenon). If the underlying regression function is itself additionally
assumed to be bounded in $L^{\infty}$ norm, then we believe the minimax rate
in their white noise setting will be $n^{-s}$, matching that in our discrete
setting. Evidence of this claim includes the result in del Álamo et al. (2021)
for the BV space (under the white noise model), who derive a rate of
$n^{-1/d}$, assuming such $L^{\infty}$ boundedness; as well as the Besov
results in Goldenshluger and Lepski (2014), who also assume $L^{\infty}$
boundedness, but study density estimation.

###  5.3 Minimax rates for linear smoothers

We now study whether linear smoothers can achieve the minimax rate over the
KTV class in (23). Before stating our result, we define a truncated eigenmaps
estimator based on $D_{n,d}^{(k+1)}$ as follows. Denote by $\xi_{i}\geq 0$,
$i\in[N]^{d}$ its singular values (noting that $(k+1)^{d}$ of these are zero),
where along each dimension in the multi-index, the singular values are sorted
in increasing order. Denote also by $v_{i}\in\mathbb{R}^{n}$, $i\in[N]^{d}$
its corresponding right singular vectors. Then, for a subset
$Q\subseteq[N]^{d}$, we denote by $V_{Q}\in\mathbb{R}^{n\times|Q|}$ the matrix
with columns given by $v_{i}$, $i\in Q$, and define the projection estimator

$\hat{\theta}=V_{Q}V_{Q}^{\mathsf{T}}y.$ (34)

Note that this reduces to the Laplacian eigenmaps estimator (here the
Laplacian is that of the grid graph) when $k=0$.

######  Theorem 5.

The minimax linear risk over the KTV class in (23) satisfies, for any sequence
$C_{n}\leq\sqrt{n}$,

$R_{L}\big{(}\mathcal{T}_{n,d}^{k}(C_{n})\big{)}=\begin{cases}\Omega(1/n+C_{n}^{2}/n)&\text{if
$s<1/2$},\\\ \Omega\big{(}1/n+C_{n}^{2}/n\log(1+n/C_{n}^{2})\big{)}&\text{if
$s=1/2$},\\\ \Omega\big{(}1/n+(C_{n}^{2}/n)^{\frac{1}{2s}}\big{)}&\text{if
$s>1/2$}.\end{cases}$ (35)

This is achieved in rate by the projection estimator in (34), by setting
$Q=[\tau]^{d}$ for $\tau^{d}\asymp(C_{n}n^{s-1/2})^{1/s}$, in the case
$s>1/2$. When $s<1/2$, the simple polynomial projection estimator, which
projects onto all multivariate polynomials of max degree $k$ (equivalently,
the estimator in (34) with $Q=[k+1]^{d}$), achieves the rate in (35). When
$s=1/2$, either estimator achieves the rate in (35) up to a log factor.
Lastly, if $C_{n}^{2}=O(n^{\alpha})$ for $\alpha<1$, and still $s=1/2$, then
either estimator achieves the rate in (35) without the additional log factor.

######  Remark 13.

Plugging in the canonical scaling $C_{n}\asymp C_{n}^{*}=n^{1-s}$ from (28),
and simplifying to dominant terms, we see that (35) becomes

$R\big{(}\mathcal{T}_{n,d}^{k}(C_{n}^{*})\big{)}\asymp\begin{cases}1&\text{if
$s\leq 1/2$},\\\ n^{-\frac{2s-1}{2s}}&\text{if $s>1/2$}.\end{cases}$

These minimax linear rates display a stark difference to the minimax rates for
the KTV class in Remark 10 (achieved by KTF up to log factors, in (32)). When
$s>1/2$, the minimax linear rate of $n^{-(2s-1)/(2s)}$ can be interpreted as
the rate of the optimal nonlinear estimator when the problem has a half degree
less of effective smoothness ($s-1/2$ in place of $s$). When $s\leq 1/2$, even
more dramatically, we see that _no linear smoother can even be consistent_
over the KTV class, in the sense of worst-case risk. This contributes an
interesting addition to the line of work on the suboptimality of linear
smoothers for nonparametric regression over heterogeneous spaces, which begun
with the seminal work of Donoho and Johnstone (1998).

###  5.4 Summary of rates

Table 1 presents a summary of the minimax rates derived in the previous three
subsections. (It offers a more detailed summary than Figure 2.) The upper
bound on minimax risk is from Corollary 1, the lower bound on the minimax risk
from Theorem 4 and Remark 10, and the minimax linear risk is from Theorem 5
and Remark 13.

Regime |  $R$, upper bound |  $R$, lower bound |  $R_{L}$, linear risk
---|---|---|---
$s<1/2$ | $n^{-s}\sqrt{\log n}$ | $n^{-s}$ | 1
$s=1/2$ | $n^{-\frac{1}{2}}\log n$ | $n^{-\frac{1}{2}}$ | 1
$s>1/2$ | $n^{-\frac{2s}{2s+1}}(\log n)^{\frac{1}{2s+1}}$ | $n^{-\frac{2s}{2s+1}}$ | $n^{-\frac{2s-1}{2s}}$
Table 1: Minimax rates over the KTV class $\mathcal{T}_{n,d}^{k}(C_{n}^{*})$,
where the canonical scaling is $C_{n}^{*}=n^{1-s}$, and recall $s=(k+1)/d$. We
use the abbreviations $R=R(\mathcal{T}_{n,d}^{k}(C_{n}^{*}))$ and
$R_{L}=R_{L}(\mathcal{T}_{n,d}^{k}(C_{n}^{*}))$.

###  5.5 Minimax rates over Sobolev classes

For completeness, we establish the minimax rate for the (discrete) Sobolev
class defined in (24). The lower bounds are simply those from the Holder class
(due to the embedding in Proposition 3), and as we show next, this is achieved
in rate by the eigenmaps estimator in (34).

######  Theorem 6.

The minimax risk over the (discrete) Sobolev class in (24) satisfies, for any
sequence $B_{n}\leq\sqrt{n}$,

$R\big{(}\mathcal{W}_{n,d}^{k+1}(B_{n})\big{)}\asymp\Omega\Bigg{(}\frac{1}{n}+\bigg{(}\frac{B_{n}^{2}}{n}\bigg{)}^{\frac{1}{2s+1}}\Bigg{)}.$

The lower bound is due to the Holder embedding in (26) (and the lower bound on
the discretized Holder class derived in Sadhanala et al. (2017)), and the
upper bound is from the estimator in (34), with $Q=[\tau]^{d}$ for
$\tau^{d}\asymp(B_{n}^{2}n^{2s})^{1/(2s+1)}$. Finally, under the canonical
scaling $B_{n}\asymp B_{n}^{*}=n^{1/2-s}$ from (27), the minimax rate
simplifies to $n^{-2s/(2s+1)}$.

##  6 Optimization algorithms

In this section, we describe a number of numerical algorithms for solving the
convex KTF problem (7), analyze their asymptotic time complexity, and
benchmark their performance. In particular, we will describe a family of
specialized ADMM algorithms that adapt to the structure of the KTF problem,
and as we will show, can find moderately accurate solutions much faster than a
general purpose “off-the-shelf” solver. This general purpose solver can be
applied to the dual of (7), which is a simple box-constrained quadratic
program, and is amenable to standard interior point methods. Details are
deferred to the appendix.

###  6.1 Specialized ADMM algorithms

Motivated by the popularity of operating splitting methods in machine learning
over the last decade, and specifically by their success in application to
trend filtering and TV denoising problems (for example, see Ramdas and
Tibshirani (2016); Wang et al. (2016); Barbero and Sra (2014); Tansey and
Scott (2015)), we investigate application of similar methods to Kronecker
trend filtering. We study a proximal Dykstra algorithm, Douglas-Rachford
splitting, and a family of specialized ADMM algorithms. For brevity, the
details on the former two are deferred to the appendix.

Our specialized ADMM approach is inspired by that of Ramdas and Tibshirani
(2016), for univariate trend filtering. To reformulate (7) into “ADMM form”,
where the criterion decomposes as a sum of functions of separate optimization
variables, we must introduce auxiliary variables. To do so, we rely on the
following observation that decomposes the KTF penalty operator into the
product of a block-diagonal matrix and a lower-order KTF penalty operator.

######  Proposition 4.

For each $j=1,2,\ldots,k+1$, the KTF penalty operator in (8) obeys (with
$D^{(0)}_{N}=I_{N}$ for convenience):

$D_{n,d}^{(k+1)}=\underbrace{\begin{bmatrix}D^{(k+1-j)}_{N}\otimes
I_{N}\otimes\cdots\otimes I_{N}&&\\\ &I_{N}\otimes
D^{(k+1-j)}_{N}\otimes\cdots\otimes I_{N}&\\\ &&\ddots&\\\ &&&I_{N}\otimes
I_{N}\otimes\cdots\otimes D^{(k+1-j)}_{N}\\\
\end{bmatrix}}_{M^{(k+1-j)}_{n,d}}D_{n,d}^{(j)}.$

This follows directly from the the univariate recursion in (4), and the
Kronecker structure in (8) (particularly, the mixed-product property of
Kronecker products), and therefore we omit its proof. Observe that each
diagonal block of the block-diagonal matrix $M^{(k+1-j)}_{n,d}$ is
itself—possibly after appropriate permutation of the row and column order—a
block-diagonal matrix with all diagonal blocks equal to $D^{(k+1-j)}_{N}$.
This is a key fact that we will leverage shortly.

Now fix any $j\in\\{1,\ldots,k+1\\}$. We can reformulate (7) as:

$\displaystyle\mathop{\mathrm{minimize}}_{\theta,z}\quad$
$\displaystyle\frac{1}{2}\|y-\theta\|_{2}^{2}+\lambda\|M^{(k+1-j)}_{n,d}z\|_{1}$
(36) $\displaystyle\mathop{\mathrm{subject\,\,to}}$ $\displaystyle
z=D_{n,d}^{(j)}\theta.$

The augmented Lagrangian associated with (36), for an augmented Lagrangian
parameter $\rho\geq 0$, is:

$L_{\rho}(\theta,z,u)=\frac{1}{2}\|y-\theta\|_{2}^{2}+\lambda\|M^{(k+1-j)}_{n,d}z\|_{1}+\frac{\rho}{2}\|z-D_{n,d}^{(j)}\theta+u\|_{2}^{2}-\frac{\rho}{2}\|u\|_{2}^{2}.$

ADMM iteratively performs a separate minimization over each of primal
variables $\theta,z$, then updates the dual variable $u$ via gradient ascent.
Namely, given some initialization $\theta^{(0)},z^{(0)},u^{(0)}$, it repeats
the following, for $t=1,2,3,\ldots$:

$\displaystyle\theta^{(t+1)}$
$\displaystyle=\Big{(}I_{n}+\rho[D_{n,d}^{(j)}]^{\mathsf{T}}D_{n,d}^{(j)}\Big{)}^{-1}\Big{(}y+\rho[D_{n,d}^{(j)}]^{\mathsf{T}}(z^{(t)}+u^{(t)})\Big{)},$
(37) $\displaystyle z^{(t+1)}$
$\displaystyle=\mathsf{prox}_{\frac{\lambda}{\rho}\|M^{(k+1-j)}_{n,d}(\cdot)\|_{1}}\Big{(}D_{n,d}^{(j)}\theta^{(t+1)}-u^{(t)}\Big{)},$
(38) $\displaystyle u^{(t+1)}$
$\displaystyle=u^{(t)}+z^{(t+1)}-D_{n,d}^{(j)}\theta^{(t+1)}.$ (39)

Here we use the notation $\mathsf{prox}_{h}(\cdot)$ for the proximal operator
associated with a function $h$. Below we make a few remarks on the
computation.

  * •

When $j=1$, the matrix $[D_{n,d}^{(j)}]^{\mathsf{T}}D_{n,d}^{(j)}$ in (37) is
the graph Laplacian of the $d$-dimensional grid, which decomposes into the
Kronecker sum of Laplacians of (univariate) chain graphs. This can be
diagonalized by a $d$-dimensional discrete cosine transform (DCT) (see, for
example, the proof of Corollary 8 in Wang et al. (2016)). Computationally,
this is essentially sequentially applying univariate DCTs to every dimension.
This implies that the $\theta$-update in (37) can be done in $O(n\log n)$
time. Further improvements (to linear-time) should be possible with multi-grid
methods.

  * •

By the key fact mentioned after the proposition, the $z$-update in (38) can be
decomposed into $dN^{d-1}$ univariate trend filtering problems, each of order
$k+1-j$ and each with sample size $N$. These can be solved in parallel, in a
total time that is nearly-linear in $n$; for example, using either the
univariate ADMM approach of Ramdas and Tibshirani (2016) or the primal-dual
interior point method (PDIP) of Kim et al. (2009). PDIP is likely the best
option for reasonably small $N$, and using it, the update (38) can be done in
$O(dN^{d-1}N^{1.5})=O(n^{1+1/(2d)})$ time.

  * •

When $j=k$ or $k+1$, the $z$-update (38) can be performed even more
efficiently, in $O(n)$ time. This is because it reduces to soft-thresholding
for $j=k+1$, and to separate univariate TV denoising problems for $j=k$. In
the former case, the linear time complexity is obvious; in the latter, it is
due to the dynamic programming (DP) of Johnson (2013).

  * •

The case $k=0$ is quite favorable, as we can use DCT in (37) and soft-
thresholding in (38), by choosing $j=1$, or simple coordinatewise shrinkage in
(37) and DP in (38), by choosing $j=0$. Both are efficient, but the latter
ends up being generally the better approach, and can be seen as the ADMM-
analog of Barbero and Sra (2014).

  * •

Among the higher-order cases $k\geq 1$, the case $k=1$ ends up being quite
special, because $j=1$ simultaneously supports the DCT solver in (37) and DP
in (38). When $k\geq 2$, we essentially need to decide in between these highly
efficient subroutines (choosing either $j=1$ or $j=k$).

In summary, each iteration (cycle of updates over $\theta,z,u$) of the
proposed ADMM algorithm is $O(n)$ for $k=0,1$. For $k\geq 2$, the time
complexity is $O(n^{1+1/(2d)})$ when we choose $j=1$, which we refer to as
ADMM _Type I_. When we choose $j=k$, which we refer to as ADMM _Type II_ , the
time complexity is dominated by the sparse linear system solve in (37). This
linear system should be well-conditioned for reasonable ranges of $\rho$, thus
the standard conjugate gradient method will be able to solve it in
approximately linear-time (proportional to the number of nonzero elements).

###  6.2 Empirical comparisons

We now compare the ADMM algorithms developed in the last subsection to
Douglas-Rachford and proximal Dykstra methods applied to (7), as well as as
general purpose solver—Gurobi, which is free for academic use—applied to the
dual of (7).333We add tiny amount of regularization to the dual problem to
avoid numerical issues that cause Gurobi to fail. The solution of this
regularized problem is first used to reconstruct the primal solution, but then
evaluated on the objective function of the original problem, when computing
the suboptimality gaps. From the family of specialized ADMM algorithms, we pay
particular attention to ADMM Types I and II, which correspond to $j=1$ and
$j=k$, respectively. We also consider ADMM Type 0, which sets $j=0$, because
it is closely related and should perform similarly to the Douglas-Rachford and
proximal Dykstra methods. In all ADMM algorithms, we adopt an adaptive choice
of $\rho$ that balances the primal and dual suboptimality (Boyd et al., 2011).

To ensure a fair comparison between ADMM Types I and II, which we will see are
generally the best performing methods (and thus the comparison between them is
of particular interest), we use optimized C++ implementations for each of
their prox subroutines; for Type I, this is the DP algorithm for univariate TV
denoising, and for Type II, this is the PDIP algorithm for univariate trend
filtering; and in both cases, we use C++ implementations from Ramdas and
Tibshirani (2016). Aside from specialized subroutines, implementation of all
iterative algorithms is in MATLAB.

The results are presented in Figures 4 and 4. Figure 4 compares the operator
splitting algorithms for denoising the standard “Lena” method at a resolution
of $256\times 256$. The KTF orders are taken to be $k=1,2,3$, corresponding to
the columns in the figure. For each $k$, the solution returned by Gurobi is
used to define the optimal criterion value, which is is then used to measure
the suboptimality gap of solutions returned the iterative methods. ADMM Type I
is generally the winner in all cases, whether measured by iteration or
especially by wall-clock time.

Figure 4 compares ADMM Type I to Gurobi for varying $k$, and also for varying
resolutions of the underlying Lena image. In all cases, ADMM Type I obtains a
moderate-quality solution in less one second—which is sometimes two orders of
magnitude faster than the off-the-shelf solver provided by Gurobi. It should
be further noted that Gurobi is highly-optimized, whereas our ADMM Type I
implementation is not—recall, only the prox subroutine is optimized, and the
outer looping is performed in MATLAB. Transporting the entire algorithm to C++
would clearly yield further improvements in efficiency. Of course, if a truly
high-accuracy solution is required, then Gurobi may be the best option.
However, its strong performance in this subsection suggests that ADMM Type I
is an efficient and useful approach for many applications in statistics and
machine learning, where moderate-accuracy solutions suffice.

Figure 3: Comparison of iterative algorithms for KTF on the standard Lena
image of resolution $256\times 256$ (that is, $n=65536$), when $k=1,2,3$,
corresponding to the three columns, from left to right. The top row compares
the convergence of the suboptimality gap as a function of the number of
iterations. The bottom row shows the same but parametrized by wall-clock time
in seconds. While these methods have similar sublinear convergence rates (top
row), ADMM Types I and II are clearly the fastest (bottom row) to reach a
small suboptimality gap, due to their low per-iteration cost. Type I is the
overall winner. (Recall that when $k=1$, Types I and II coincide, so the blue
curves is hidden behind the red curve).

Figure 4: Comparison of ADMM Type I to Gurobi for the same Lena problem. The
left panel compares the two methods for varying $k$ at a fixed resolution of
$256\times 256$. The middle panel fixes $k=2$ and compares them across varying
resolutions. The right panel plots the time needed to achieve a certain
relative error, defined by the ratio of the suboptimality gap to the objective
value being less than $10^{-2}$. Altogether, we see that ADMM Type I achieves
a moderate-quality solution several orders of magnitude faster than Gurobi.
Furthermore, the right panel shows that Gurobi appears to scale as
$O(n^{1.5}$) (to be expected, if it is based on interior point methods
internally) whereas ADMM Type I appears to scale closer to $O(n)$.

##  7 Interpolation algorithm

In this section, we derive an algorithm to extend the KTF solution in (7),
defined only at points in the lattice $Z_{n,d}$, to a function defined on all
of $[0,1]^{d}$. As explained in Remark 3, this is made possible by the
continuous-time formulation for KTF in (12) of Proposition 2, which, recall,
relates to the original discrete-time problem (7) in that at their solutions,
we have $\hat{\theta}_{i}=\hat{f}(x_{i})$, for $i=1,\ldots,n$. Given the
coefficients that define the function $\hat{f}\in\mathcal{H}_{n,d}^{k}$ in its
expansion in terms of the tensor product of univariate falling factorial basis
functions, that is, given the solution $\hat{\alpha}$ in (13), we can form the
interpolated prediction $\hat{f}(x)$ at an arbitrary point $x\in[0,1]^{d}$ by
evaluating the basis expansion at $x$, as in (14).

While this is conceptually the easiest way to interpolate the fitted values
$\hat{f}(x_{i})$, $i=1,\ldots,n$ to form the prediction at an arbitrary
$x\in[0,1]^{d}$, it is not the most efficient. The expression in (14) takes
$O(\|\hat{\alpha}\|_{0}(k+1)^{d})$ operations, where $\|\hat{\alpha}\|_{0}$
denotes the number of nonzero elements in $\hat{\alpha}$ (the number of active
basis functions), because evaluating each univariate basis function
$h_{i_{j}}(x_{j})$ takes $O(k+1)$ operations (it being a product of $k+1$
terms, recall (10)). In what follows, we present an interpolation algorithm
that takes only $O((k+1)^{d+1})$ operations, a big savings over (14) when
$\|\hat{\alpha}\|_{0}$ is large. Moreover, our algorithm acts directly on the
solution $\hat{\theta}$ in (7), meaning that we never have to solve (13) in
the first place.

###  7.1 Review: univariate interpolation

Our interpolation algorithm for KTF estimates in the multivariate case builds
off the univariate discrete spline interpolation algorithm derived in
Corollary 2 of Tibshirani (2020). For convenience, we transcribe this in
Algorithm 1. Here and henceforth, we use the abbreviation
$x_{a:b}=(x_{a},\ldots,x_{b})$ for integers $a\leq b$. Also, we use
$f[z_{1},\ldots,z_{r}]$ for the divided difference of a function $f$ at
distinct points $z_{1},\ldots,z_{r}$. Recall, this is defined for $r=2$ by

$f[z_{1},z_{2}]=\frac{f(z_{2})-f(z_{1})}{z_{2}-z_{1}},$

and for any $r\geq 3$ by the recursion

$f[z_{1},\ldots,z_{r}]=\frac{f[z_{2},\ldots,z_{r}]-f[z_{1},\ldots,z_{r-1}]}{z_{r}-z_{1}}.$

Note that for evenly-spaced points, this simply coincides with a scaled
forward difference; in particular, recalling the notation introduced in
Section 1.1, we have

$f[z,\ldots,z+(r-1)/n]=\frac{(r-1)!}{n^{r}}(\Delta^{r}f)(z).$

For more background on divided differences, discrete splines, their connection
to the falling factorial basis and to trend filtering, we refer to Tibshirani
(2020).

Algorithm 1 takes $O((k+1)^{2})$ operations, as (40), (41) are each linear
systems in just one unknown, and forming the coefficients in either linear
system can be done in $O((k+1)^{2})$ operations (this uses a representation of
a divided difference as an explicit linear combination of the underlying
function evaluations; see, for example, Section 2.1 of Tibshirani (2020)).
Note that this assumes $x_{1:n}$ are evenly-spaced design points, because in
this case identifying the smallest index $i$ such that $x_{i}>x$ can be done
with integer division. For general design points, the identification step
takes $O(\log n)$ operations via binary search, so the total cost would be
$O(\log n+(k+1)^{2}$).

Moreover, Corollary 2 in Tibshirani (2020) establishes that the value $f(x)$
returned by Algorithm 1 is equal to that produced by the falling factorial
basis representation in (11) (for $d=1$), where $\alpha$ is the unique
coefficient vector such that $f(x_{i})=\theta_{i}$, $i=1,\ldots,n$.

Input: design points $x_{1:n}$ with entries in increasing order; values
$\theta_{1:n}$ to interpolate; query point $x$; integer $k\geq 0$.
Output: interpolated value $f(x)$, where $f$ is the unique $k^{\text{th}}$
order discrete spline with knots in $x_{(k+1):(n-1)}$, such that
$f(x_{i})=\theta_{i}$, $i=1,\ldots,n$.

  1. 1.

If $x=x_{i}$ for some $i=1,\ldots,n$, then return $\theta_{i}$.

  2. 2.

Else, if $x>x_{k+1}$ and $i$ is the smallest index such that $x_{i}>x$ (with
$i=n$ when $x>x_{n}$), then return $f(x)$ as the unique solution of the linear
system:

$f[x_{i-k},\ldots,x_{i},x]=0.$ (40)

  3. 3.

Else, if $x<x_{k+1}$, then return $f(x)$ as the unique solution of the linear
system:

$f[x_{1},\ldots,x_{k+1},x]=0.$ (41)

(Note that both (40), (41) are linear systems in just one unknown, $f(x)$,
since we interpret $f(x_{i})=\theta_{i}$, $i=1,\ldots,n$.)

Algorithm 1 Interpolate-1d$(x_{1:n},\theta_{1:n},x,k)$

Input: lattice
$\\{z_{i1}\\}_{i=1}^{N_{1}}\times\cdots\times\\{z_{id}\\}_{i=1}^{N_{d}}$ where
each set $\\{z_{ij}\\}_{i=1}^{N_{j}}$ in the Cartesian product is sorted in
increasing order; values
$\\{\theta_{i}\\}_{i\in[N_{1}]\times\cdots\times[N_{d}]}$ over the lattice to
interpolate; query point $x$; integer $k\geq 0$.
Output: interpolated value $f(x)$, for the unique function $f$ in the tensor
product space of $k^{\text{th}}$ degree discrete splines with knots in
$\\{z_{i1}\\}_{i=k+1}^{N_{1}-1}\times\cdots\times\\{z_{id}\\}_{i=k+1}^{N_{d}-1}$,
such that $f(z_{i_{1},1},\ldots,z_{i_{d},d})=\theta_{i_{1},\ldots,i_{d}}$,
$(i_{1},\ldots,i_{d})\in[N_{1}]\times\cdots\times[N_{d}]$.

  1. 1.

If $d=1$, then return Interpolate-1d$(z_{1:N_{1}},\theta_{1:N_{1}},x,k)$.

  2. 2.

Else, let $i_{1}$ denote the smallest index such that $x_{i_{1},1}\geq
z_{i_{1},1}$.

  3. 3.

Let $\ell_{1}=\min\\{\max\\{i_{1}-k,1\\},N_{1}-k\\}$.

  4. 4.

Let
$\vartheta_{p}=\textsc{Interpolate}\big{(}\\{z_{i2}\\}_{i=1}^{N_{2}}\times\cdots\times\\{z_{id}\\}_{i=1}^{N_{d}},\,\\{\theta_{i}\\}_{i\in\\{i_{1}+p-1\\}\times[N_{2}]\times\cdots\times[N_{d}]},\,x_{2:d},\,k\big{)}$,
for $p\in[k+1]$.

  5. 5.

Return
$\textsc{Interpolate-1d}(z_{\ell_{1}:(\ell_{1}+k),1},\vartheta_{1:(k+1)},x_{1},k)$.

Algorithm 2
Interpolate$\big{(}\\{z_{i1}\\}_{i=1}^{N_{1}}\times\cdots\times\\{z_{id}\\}_{i=1}^{N_{d}},\,\\{\theta_{i}\\}_{i\in[N_{1}]\times\cdots\times[N_{d}]},\,x,\,k\big{)}$

###  7.2 Multivariate interpolation

In the multivariate case, it turns out that we can interpolate within the
space of tensor products of $k^{\text{th}}$ degree discrete splines in
$O((k+1)^{d+1})$ time, assuming a uniformly-spaced lattice. The idea is to
reduce the $d$-dimensional problem calculation to $k+1$ interpolation
problems, each in dimension $d-1$. Figure 6 provides the intuition for $d=2$.
The algorithm in described in Algorithm 2. We recall the notation from Section
5, where we abbreviate $[a]=\\{1,\ldots,a\\}$ for an integer $a\geq 1$.

Algorithm 2 assumes each side length of the lattice is at least $k+1$. Its
proof of correctness, as well as the running time of $O((k+1)^{d+1})$ for a
uniformly-spaced lattice, follows from a straightforward inductive argument
over $d$, whose proof we omit. We note that for an arbitrary lattice, the
running time is $O(\sum_{j=1}^{d}\log N_{j}+(k+1)^{d+1})$. Figure 6 gives some
examples of interpolation for $d=2$.

Figure 5: Illustration of multivariate interpolation from Algorithm 2, when
$d=2$ and $k=2$. The value to be interpolated is marked by a red star. The
algorithm first interpolates the $k+1=3$ values indicated by red dots, each
time using univariate interpolation along the y-axis, from Algorithm 1. As the
final step, these red dots are used to interpolate the value at the red star,
using univariate interpolation along the x-axis, again from Algorithm 1.

Figure 6: Top left: KTF solution, with $k=1$, for a problem over a square
lattice with side length $N=25$ (that is, $n=25^{2}=625$ points). Top right:
the interpolated surface from Algorithm 2 (itself evaluated over a grid with
4x the resolution along each dimension; that is, $N=100$). Bottom row:
analogous but for $k=2$. In either case the interpolated function is in the
tensor space of $k^{\text{th}}$ order discrete splines.

##  8 Experiments

We explore the empirical properties of KTF through a series of experiments
that compare its performance against other nonparametric methods.

###  8.1 Comparison of methods

We study the performance of KTF against other nonparametric estimators by
analyzing their empirical mean squared error (MSE) in estimating the function
$f_{0}$ defined in Figure 1, over a square lattice with $n=256^{2}=65536$
points. In the experiments that follow, we generate data by adding Gaussian
noise to evaluations of $f_{0}$, of the same magnitude illustrated in Figure
1, yielding a signal-to-noise ratio (SNR) of 0.5. Aside from KTF and graph
trend filtering (GTF) of orders $k=0,1,2$ (for $k=0$, they both reduce to TV
denoising), we also consider second-order Laplacian smoothing (using the
squared Laplacian $L^{2}$, where $L$ is the Laplacian matrix of the 2d grid
graph), the eigenmaps estimator in (34) with $k=1$, kernel smoothing (using a
Gaussian kernel), and wavelet smoothing (using Daubechies’ least asymmetric
wavelets, ten vanishing moments, and hard thresholding). For the latter two
estimators, we use the implementations from the R packages np and wavethresh,
respectively.

Figure 7: MSE of various methods for estimating the signal function $f_{0}$ in
Figure 1, over a square lattice with $n=256^{2}=65536$ points, and subject to
Gaussian noise to yield an SNR of 0.5. Left panel: the average performance
(over 20 repetitions) of trend filtering estimators over a common range of
tuning parameters $\lambda$. Right panel: the average performance of other
nonparametric methods. Since their tuning parameters lie on different scales,
the tuning parameter values are scaled to fit onto a single x-axis, and we
notate the MSE-optimal tuning parameter value for each method. The takeaway is
that KTF and GTF, in particular for $k=2$, achieve clearly the best
performance.

Each estimator is fit over a range of tuning parameters, and its MSE as a
function of the tuning parameter is shown in Figure 7. The MSE curves here are
averaged over 20 repetitions (each repetition forms a data set by adding noise
to $f_{0}$). Error bars are shown as well, denoting the standard deviations of
the MSE over these 20 repetitions. We can see that the trend filtering
estimators (with the exception of TV denoising, which returns a piecewise
constant fit that is not well-suited to an underlying signal with this level
of smoothness) perform quite a bit better than all competitors, and factoring
in the variability over repetitions, KTF and GTF achieve essentially
equivalent performance (for common values of $k$). The superiority of KTF over
the linear estimators (Laplacian smoothing, eigenmaps projection, and kernel
smoothing), should not come as a surprise—our theory prescribes that KTF
should perform better in a minimax sense than any linear smoother when the
underlying signal displays heterogeneous smoothness. These results thus serve
as a quantitative complement to the qualitative findings in Figure 1, and the
theoretical findings in Section 5.

###  8.2 Rates for heterogeneous signals

We now examine the empirical error rates of KTF and a linear smoother in
estimating signals belonging to KTV classes, for $k=0,1$. In both cases, we
choose $k,d$ so that the effective degree of smoothness remains $s=1/2$, and
we use the canonical scaling in (28), so that the two classes under
consideration are:

$\mathcal{T}_{n,2}^{0}(\sqrt{n})\quad\text{and}\quad\mathcal{T}_{n,4}^{1}(\sqrt{n}).$
(42)

As representatives for “hard” signals within these two classes, we take the
true mean $\theta_{0}$ to be a “one-hot” signal when $k=0$, and a linear
“spike” signal when $k=1$ (each scaled to the appropriate magnitude). For each
sample size $n$, we compute the MSE of KTF, and the eigenmaps projection
estimator in (34) with $k=0,1$, averaged over 20 repetitions, where each
method is tuned to have the optimal average MSE over a range of tuning
parameter values.

Figure 8 reports these averaged MSE curves as functions of $n$, where the
error bars again denote standard deviations. We can see that in both cases,
the KTF error decays faster than the minimax rate (suggesting that the
particular signals under consideration do not achieve the worst-case rate for
KTF), while the linear method exhibits a perfectly flat error curve,
suggesting that it fails to be consistent entirely.

Figure 8: Empirical error rates of KTF versus the eigenmaps projection
estimator, for two cases: $k=0,\,d=2$, and $k=1,\,d=4$. The true mean in each
case was chosen to be representative of a “hard” signal in the appropriate KTV
class in (42). KTF converges faster than the minimax rate, whereas the
eigenmaps estimator fails to be consistent entirely.

###  8.3 Rates for homogeneous signals

Lastly, we examine the empirical error rates of KTF and the eigenmaps
projector for estimating homogeneously smooth signals. Recall that we showed
the latter is minimax rate optimal over Sobolev (and Holder) classes in
Theorem 6, whereas for $s>1/2$, it is not clear (from its minimax rate over
the inscribing KTV class) that KTF achieves the faster rate for Sobolev (or
Holder) signals. We fix $k=0$ and consider two values for the dimension, $d=2$
and $d=3$, which correspond to $s=1/2$ and $s>1/2$, respectively. The true
mean $\theta_{0}$ was taken to be the evaluations of a linear function over
the lattice, scaled appropriately (according to the canonical scaling). Figure
9 reports the MSE curves from each method, under the same general setup as in
the last subsection (averaged over 20 repetitions, optimal tuning per $n$).
Moving from $s=1/2$ to $s>1/2$, we do not find that the empirical performance
of TV denoising becomes markedly worse than its linear competitor, and the
adaptivity of KTF to the smooth signals remains an open question.

Figure 9: Empirical error rates of KTF versus the eigenmaps projection
estimator, for two cases: $k=0,\,d=2$, and $k=0,\,d=3$. The true mean was
chosen to be a linear function in either case. KTF does not appear to be
outperformed (in rate) by the eigenmaps estimator, leaving its adaptivity to
smooth signals when $s>1/2$ an open question.

##  9 Discussion

###  9.1 General lattice structures

The results and analysis in this paper have heretofore assumed an equally
spaced lattice with the same number of knots in each coordinate direction, but
as alluded to in the introduction, many of our results also apply to a more
general lattice structure, namely a Cartesian product
$Z_{N_{1},\dotsc,N_{d}}=\\{z_{i1}\\}_{i=1}^{N_{1}}\times\\{z_{i2}\\}_{i=1}^{N_{2}}\times\cdots\times\\{z_{id}\\}_{i=1}^{N_{d}}$,
with $n=\prod_{j=1}^{d}N_{j}$. With this general lattice, we generalize the
$(k+1)^{\text{th}}$-order discrete derivative matrix (8),

$D^{(k+1)}_{z_{1},\dotsc,z_{d}}=\left[\begin{array}[]{c}D^{(k+1)}_{z_{1}}\otimes
I_{N_{2}}\otimes\cdots\otimes I_{N_{d}}\\\ I_{N_{1}}\otimes
D^{(k+1)}_{z_{2}}\otimes\cdots\otimes I_{N_{d}}\\\ \vdots\\\ I_{N_{1}}\otimes
I_{N_{2}}\otimes\cdots\otimes D^{(k+1)}_{z_{d}}\end{array}\right],$ (43)

where $z_{j}:=\\{z_{ij}\\}_{i=1}^{N_{j}}$ and $D_{z_{j}}$ is a univariate
discrete derivative matrix posed over that set of points. When the lattice
points may be embedded into Euclidean space, i.e., the data have been rotated
such that the elements of $z_{j}$ are increasing and vary only in the
$j^{\text{th}}$ coordinate, then $D^{(k+1)}_{z_{j}}$ is given (Tibshirani,
2014, 2020) by

$D^{(k+1)}_{z_{j}}:=D^{(1)}\cdot\mathrm{diag}\left(\frac{k}{z_{(k+1)j}-z_{1j}},\frac{k}{z_{(k+2)j}-z_{2j}},\cdots,\frac{k}{z_{N_{j}j}-z_{(N_{j}-k)j}}\right)\cdot
D^{(k)},$ (44)

where we have abused notation to only consider the $j^{\text{th}}$ coordinate
in the $z_{ij}$ differences (as they are the only nonzero differences).

The KTF estimator with (43) satisfies Proposition 1, which states that the
nullspace of (43) consists of polynomials with maximum degree $k$. Similarly,
we obtain a version of Proposition 2, with appropriate adjustments for the
number of basis elements oriented in each coordinate direction, and since KTF
on a general lattice is still a generalized lasso problem, we inherit results
on the degrees of freedom. Finally, we expect that under mild assumptions on
the spacing of the points $z_{j}$ in each direction $j=1,\dotsc,d$, one also
obtains the discrete Holder class embedding of Proposition 3 and the
estimation error bounds of Theorem 3.

###  9.2 Scattered data

The continuous-time representation of KTF, elucidated in Section 2.3, suggests
an approach for fitting a KTF estimate to scattered data. Recall that
Proposition 2 provides that the KTF problem (7) is equivalent to a variational
problem over a function space spanned by tensor product of falling factorial
basis functions, penalized by a kind of anisotropic total variation norm.
Moreover, this “synthesis” viewpoint of KTF has a basis regression form made
explicit in (13). Given scattered data
$(x_{1},y_{1}),\dotsc,(x_{n},y_{n})\in\mathbb{R}^{d}\times\mathbb{R}$ and
assuming a knot set
$Z_{N_{1},\dotsc,N_{d}}=\\{z_{i1}\\}_{i=1}^{N_{1}}\times\\{z_{i2}\\}_{i=1}^{N_{2}}\times\cdots\times\\{z_{id}\\}_{i=1}^{N_{d}}$,
we let

$h_{z_{j},x_{i}}:=(h_{N_{j},1}(x_{ij}),h_{N_{j},2}(x_{ij}),\dotsc,h_{N_{j},N_{j}}(x_{ij}))\in\mathbb{R}^{N_{j}}$
(45)

be the evaluations of the marginal falling factorial bases (with knots given
by $z_{j}=\\{z_{ij}\\}_{i=1}^{N_{j}}$; see Tibshirani (2020) for the unevenly
spaced case) in coordinate direction $j$ for the $i^{\text{th}}$ observation.
Then we solve

$\mathop{\mathrm{minimize}}_{\alpha\in\mathbb{R}^{\prod_{j=1}^{d}N_{j}}}\frac{1}{2}\Bigg{\lVert}y-\underbrace{\begin{bmatrix}h_{z_{1},x_{1}}\otimes
h_{z_{2},x_{1}}\otimes\dots\otimes h_{z_{d},x_{1}}\\\ h_{z_{1},x_{2}}\otimes
h_{z_{2},x_{2}}\otimes\dots\otimes h_{z_{d},x_{2}}\\\ \vdots\\\
h_{z_{1},x_{n}}\otimes h_{z_{2},x_{n}}\otimes\dots\otimes
h_{z_{d},x_{n}}\end{bmatrix}}_{=:\tilde{H}^{(k+1)}}\alpha\Bigg{\rVert}_{2}^{2}+\lambda\left\lVert
D^{(k+1)}_{z_{1},\dotsc,z_{d}}\begin{bmatrix}H_{z_{1}}^{(k+1)}\otimes
H_{z_{2}}^{(k+1)}\otimes\dots
H_{z_{d}}^{(k+1)}\end{bmatrix}\alpha\right\rVert_{1},$ (46)

and take $\tilde{H}^{(k+1)}\hat{\alpha}$ to be fitted values for a given
$\lambda$, where $D^{(k+1)}_{z_{1},\dotsc,z_{d}}$ is defined as in (43) and
$H_{z_{j}}^{(k+1)}$ is the $(k+1)^{\text{th}}$ order falling factorial basis
with knots given by $z_{j}$. Since the scattered data do not immediately imply
an appropriate lattice of knots with which to form the estimator (46), we
propose a simple scheme, absent countervailing considerations. Take
$N_{1}=\dots=N_{d}=N:=\lfloor n^{1/d}\rfloor$, and project the data
$\\{x_{i}\\}_{i=1}^{n}$ onto each of the $d$ coordinate axes to obtain
$\\{p_{ij}\\}_{i=1}^{n}\subset\mathbb{R}$, $j=1,\dotsc,d$. Then, for each $j$,
let the knots $\\{z_{ij}\\}_{i=1}^{N}$ be the $1/(N+1),2/(N+1),\dotsc,N/(N+1)$
quantiles of $\\{p_{ij}\\}_{i=1}^{n}$.

###  9.3 Generalized linear models

The KTF problem may be further extended through the replacement of the squared
loss in (7) with a generalized linear model (GLM) loss. Specifically, in lieu
of (7), we solve

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;-\theta^{\top}y+A(\theta)+\lambda\|D_{n,d}^{(k+1)}\theta\|_{1},$
(47)

where $y_{i}$ is assumed to be distributed according to an exponential family
distribution with canonical parameter $\theta$ and log-partition function
$A(\theta)$. This reduces to the squared loss when $y_{i}$ follows a Gaussian
distribution with mean $\mu_{i}$ and variance $\sigma^{2}$. A common
application of the GLM loss framework is density estimation via Poisson
regression. Supposing observations from a distribution supported on a
hyperrectangle, one discretizes the domain using a grid and bins the counts.
The centers of the bins may then be treated as lattice points and their
corresponding counts modeled using a KTF GLM where the response is a Poisson
mean. Specifically, letting $y$ be a vector of binned counts and $\theta$ be
the canonical parameter, we solve

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;-\theta^{\top}y+\mathds{1}^{\top}\exp{\theta}+\lambda\|D_{n,d}^{(k+1)}\theta\|_{1}.$
(48)

This learns a spatially dependent Poisson mean that varies smoothly over the
domain, with the notion of smoothness dictated by $D_{n,d}^{(k+1)}$. We note
that this has previously been studied in univariate and a graph fused lasso
setting in Padilla and Scott (2016); Bassett and Sharpnack (2019); this
proposal would extend the approach to multiple dimensions and smoother
signals. The theory for GLMs is known to be more complicated (see, e.g., Haris
et al., 2019), but some simple truncation arguments may be a viable approach
(e.g., Lin et al., 2017).

###  9.4 Mixed discrete derivatives

The KTF penalty defined via (8) applies an $\ell_{1}$ penalty to the
$(k+1)^{\text{th}}$ discrete derivatives, marginally along each coordinate
direction. One could instead penalize mixed discrete derivatives of maximum
order $k+1$, in which case the nullspace of the penalty operator is
polynomials of _total_ degree $k$ rather than those of _maximum_ degree $k$.
The tradeoff is that this is computationally more difficult since the penalty
matrix becomes more complex. In particular, the penalty for a
$\theta\in\mathbb{R}^{n}$ is given by

$\sum_{\alpha:\alpha^{\mathsf{T}}1=k+1,\alpha\geq 0}\;\;\sum_{x\in
Z_{n,d}}\big{(}\Delta_{x_{1}^{\alpha_{1}}\dots
x_{d}^{\alpha_{d}}}\theta\big{)}(x)$

in the notation of Section 1.2. It would be interesting to verify that this
penalty formulation matches with the obvious continuous equivalent, in the
sense of Proposition 2 for KTF penalty. Error analysis for this method is also
left to future work.

#### Different smoothness levels across dimensions.

We can naturally extend KTF in (6) to handle differing orders of smoothness
across dimensions:

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;\frac{1}{2}\sum_{i=1}^{n}\big{(}y_{i}-\theta(x_{i})\big{)}^{2}+\sum_{j=1}^{d}\lambda_{j}\sum_{x\in
Z_{n,d}}|(\Delta_{x_{j}^{k_{j}+1}}\theta)(x)|.$ (49)

with $d$ tuning parameters $\lambda_{1},\dots,\lambda_{d}$. This estimator
satisfies key properties in Section 2 and the optimization and interpolation
algorithms in Sections 6, 7 apply to this method with appropriate change in
notation. We believe estimation error results similar in spirit to those in
Section 5 will hold, but we haven’t pursued the analysis. Two key pre-
conditions to apply Theorem 2 are incoherence of eigenvectors of
$DD^{\mathsf{T}}$ and tight bounds on the summation $\sum_{i}\xi_{i}^{-2}$.
Incoherence is easy to verify for the difference operator in (49) but we leave
bounding the summation to future work.

##  Appendix A Proofs

###  A.1 Proof of Proposition 1

Abbreviating $D=D_{n,d}^{(k+1)}$, the null space of $D$ is the number of its
nonzero singular values or equivalently, the number of nonzero eigenvalues of
$D^{\mathsf{T}}D$. Following from (8), and abbreviating $Q=D^{(k+1)}_{N}$ and
$I=I_{N}$,

$D^{\mathsf{T}}D=Q^{\mathsf{T}}Q\otimes I\otimes\cdots\otimes I\;+\;I\otimes
Q^{\mathsf{T}}Q\otimes\cdots\otimes I\;+\;\ldots\;+\;I\otimes
I\otimes\cdots\otimes Q^{\mathsf{T}}Q,$

the Kronecker sum of $Q^{\mathsf{T}}Q$ with itself, a total of $d$ times.
Using a standard fact about Kronecker sums, if we denote by $\rho_{i}$,
$i=1,\ldots,N$ the eigenvalues of $Q^{\mathsf{T}}Q$ then

$\rho_{i_{1}}+\rho_{i_{2}}+\cdots+\rho_{i_{d}},\quad
i_{1},\ldots,i_{d}\in\\{1,\ldots,N\\}$

are the eigenvalues of $D^{\mathsf{T}}D$. By counting the multiplicity of the
zero eigenvalue, we arrive at a nullity for $D$ of $(k+1)^{d}$. It is
straightforward to check that the vectors specified in the proposition, given
by evaluations of polynomials of max degree $k$, are in the null space, and
that these are linearly independent, which completes the proof. ∎

###  A.2 Proof of Proposition 2

Let us define

$B^{(k+1)}_{N}=\begin{bmatrix}C^{(k+1)}_{N}\\\
D^{(k+1)}_{N}\end{bmatrix}\in\mathbb{R}^{N\times N},$

where the first $k+1$ rows are given by a matrix
$C^{(k+1)}_{N}\in\mathbb{R}^{(k+1)\times N}$ that completes the row space, as
in Lemma 2 of Wang et al. (2014), or Section 6.2 of Tibshirani (2020). Now,
again by Lemma 2 of Wang et al. (2014), or Section 6.3 of Tibshirani (2020),

$\big{(}H^{(k+1)}_{N}\big{)}^{-1}=\frac{1}{k!}B^{(k+1)}_{N}$ (50)

where $H^{(k+1)}_{N}\in\mathbb{R}^{N\times N}$ is the falling factorial basis
matrix of order $k$, which has elements

$\big{[}H^{(k+1)}_{N}\big{]}_{ij}=h^{k}_{N,j}(i/N),\quad i,j=1,\ldots,N,$

with $h^{k}_{N,i}$, $i=1,\ldots,N$ denoting the falling factorial functions in
(10) with respect to design points $1/N,2/N,\ldots,1$.

We now transform variables in (7) by defining

$\theta=\Big{(}H^{(k+1)}_{N}\otimes\cdots\otimes H^{(k+1)}_{N}\Big{)}\alpha,$

and using (50), this turns (7) into an equivalent basis form,

$\mathop{\mathrm{minimize}}_{\alpha\in\mathbb{R}^{n}}\;\frac{1}{2}\Big{\|}y-\Big{(}H^{(k+1)}_{N}\otimes\cdots\otimes
H^{(k+1)}_{N}\Big{)}\alpha\Big{\|}_{2}^{2}+\lambda
k!\left\|\left[\begin{array}[]{c}I^{0}_{N}\otimes
H^{(k+1)}_{N}\otimes\cdots\otimes H^{(k+1)}_{N}\\\ H^{(k+1)}_{N}\otimes
I^{0}_{N}\otimes\cdots\otimes H^{(k+1)}_{N}\\\ \vdots\\\ H^{(k+1)}_{N}\otimes
H^{(k+1)}_{N}\otimes\cdots\otimes
I^{0}_{N}\end{array}\right]\alpha\right\|_{1},$

where $I^{0}_{N}\in\mathbb{R}^{(N-k-1)\times N}$ denotes the last $N-k-1$ rows
of the identity $I_{N}$. We can rewrite the problem once more by parametrizing
the evaluations according to $f$ in (11), which we claim yields (12). The
equivalence between loss terms in the above problem and (12) is immediate (by
definition of $f$); to see the equivalence between penalty terms, it can be
directly checked that

$k!\Big{(}I^{0}_{N}\otimes H^{(k+1)}_{N}\otimes\cdots\otimes
H^{(k+1)}_{N}\Big{)}\alpha$

contains the differences of the function $\partial^{k}f/\partial x_{1}^{k}$
over all pairs of grid positions that are adjacent in the $x_{1}$ direction,
where $f$ is as in (11). This, combined with the fact that
$\partial^{k}f/\partial x_{1}^{k}$ is constant in between lattice positions,
means that

$k!\Big{\|}\Big{(}I^{0}_{N}\otimes H^{(k+1)}_{N}\otimes\cdots\otimes
H^{(k+1)}_{N}\Big{)}\alpha\Big{\|}_{1}=\sum_{x_{-1}}\mathrm{TV}\bigg{(}\frac{\partial^{k}f(\cdot,x_{-1})}{\partial
x_{1}^{k}}\bigg{)},$

the total variation of $\partial^{k}f/\partial x_{1}^{k}$ added up over all
slices of the lattice $Z_{n,d}$ in the $x_{1}$ direction. Similar arguments
apply to the penalty terms corresponding to dimensions $j=2,\ldots,d$, and
this completes the proof. ∎

###  A.3 Proof of Theorem 1

Denote by $f_{\epsilon}=\eta_{\epsilon}*f$ the mollified version of $f$, where
$\eta_{\epsilon}(x)=\epsilon^{-d}\eta(x/\epsilon)$, and
$\eta:\mathbb{R}^{d}\to\mathbb{R}$ is the standard mollifier, defined by

$\eta(x)=\begin{cases}\displaystyle
c\exp\bigg{(}\frac{1}{\|x\|_{2}^{2}-1}\bigg{)}&\text{if $\|x\|_{2}\leq 1$},\\\
0&\text{else},\end{cases}$

and $c>0$ is a normalization constant so that $\eta$ integrates to 1. By
construction, for any $\epsilon>0$, we have $f_{\epsilon}\in C^{\infty}(U)$.
Now from (18), simply exchanging the sum over absolute partial derivatives and
the integral, we have

$\displaystyle\mathrm{TV}(f_{\epsilon};U)$
$\displaystyle=\sum_{j=1}^{d}\int_{U}\bigg{|}\frac{\partial
f_{\epsilon}(x)}{\partial x_{j}}\bigg{|}\,dx$
$\displaystyle=\sum_{j=1}^{d}\int_{U_{-j}}\int_{I_{x_{-j}}}\bigg{|}\frac{\partial
f_{\epsilon}(x_{j},x_{-j})}{\partial x_{j}}\bigg{|}\,dx_{j}\,dx_{-j}$
$\displaystyle=\sum_{j=1}^{d}\int_{U_{-j}}\int_{I_{x_{-j}}}\bigg{|}\frac{\partial
f_{\epsilon}(x_{j},x_{-j})}{\partial x_{j}}\bigg{|}\,dx_{j}\,dx_{-j}$
$\displaystyle=\sum_{j=1}^{d}\int_{U_{-j}}\mathrm{TV}\big{(}f_{\epsilon}(\cdot,x_{-j});I_{x_{-j}}\big{)}\,dx_{-j},$

where in the last line we applied the representation of TV for smooth
functions in (18), but to the univariate function $x_{j}\mapsto
f_{\epsilon}(x_{j},x_{-j})$, for fixed $x_{-j}$. Recalling standard results on
approximation of BV functions by smooth functions (see, for example, Theorem
5.22 in Evans and Gariepy (2015)), by sending $\epsilon\to 0$, we have that
the left-most and right-most sides of the previous display approach those in
(20), completing the proof.

###  A.4 Proof of Theorem 3

Abbreviate $N^{\prime}=N-k-1$. Let $\beta_{i},u_{i},v_{i}$ be a triplet of
nonzero singular value, left singular vector, and right singular vector of
$D^{(k+1)}_{N,1}$, for $i\in[N^{\prime}]$ and let $p_{j}$, $j\in[k+1]$ form an
orthogonal basis for the null space of $D^{(k+1)}_{N,1}$. From Lemma 1 it
suffices to show incoherence of $u_{i},v_{i}$, $i\in[N^{\prime}]$, and
$p_{i}$, $i\in[k+1]$. Incoherence of $u_{i}$, $i\in[N^{\prime}]$ and $v_{i}$,
$i\in[N^{\prime}]$is established in Sadhanala et al. (2017) . Incoherence of
$p_{i}$, $i\in[k+1]$ may be seen by choosing, e.g., these vectors to be the
discrete Legendre orthogonal polynomials as in Neuman and Schonbach (1974).
Applying Lemma 1, we can see that $D_{n,d}^{(k+1)}$ satisfies the incoherence
property, as defined in Theorem 2.

From the incoherence property and Theorem 2, the KTF estimator $\hat{\theta}$,
satisfies

$\displaystyle\frac{1}{n}\|\hat{\theta}-\theta_{0}\|_{2}^{2}=O_{\mathbb{P}}\Bigg{(}\frac{\kappa}{n}+\frac{|I|}{n}+\frac{\mu}{n}\sqrt{\frac{\log
n}{n}\sum_{i\in[N]^{d}\setminus(I\cup[k+1]^{d})}\frac{1}{\xi_{i}^{2}}}\cdot\|\Delta\theta_{0}\|_{1}\,\Bigg{)},$
(51)

where we abbreviate $\Delta=D_{n,d}^{(k+1)}$, $\xi_{i},i\in[N]^{d}$ are
eigenvalues of $\Delta^{\mathsf{T}}\Delta$ with $\xi_{i}=0$ for
$i\in[k+1]^{d}$. We reindexed the eigenvalues so that they correspond to grid
positions.

Recall the shorthand $s=(k+1)/d$. For $s\leq 1/2$, set
$I=[k+2]^{d}\setminus[k+1]^{d}$. From Lemma 6,

$\displaystyle\sum_{i\in[N]^{d}\setminus(I\cup[k+1])^{d}}\frac{1}{\xi_{i}^{2}}\leq
c\begin{cases}n&s<1/2\\\ n\log n&s=1/2.\end{cases}$

Plugging this into (51) gives the desired bounds when $s\leq 1/2$:

$\displaystyle\frac{1}{n}\|\hat{\theta}-\theta_{0}\|_{2}^{2}=\begin{cases}O_{\mathbb{P}}\Big{(}\frac{(k+2)^{d}}{n}+\|\Delta\theta_{0}\|_{1}\sqrt{\log
n}\Big{)}&s<1/2\\\
O_{\mathbb{P}}\Big{(}\frac{(k+2)^{d}}{n}+\|\Delta\theta_{0}\|_{1}\log
n\Big{)}&s=1/2.\end{cases}$

Now consider the case $s>1/2$. Set
$I=\\{i\in[N]^{d}:\|(i-k-2)_{+}\|_{2}<r\\}\setminus[k+1]^{d}$ for an
$r\in[1,\sqrt{d}N]$ to be chosen later. $|I|\leq(r+k+1)^{d}$ because
$I\subseteq\big{[}\lfloor r\rfloor+k+1\big{]}^{d}$. Lemma 6 shows that for a
constant $c>0$ depending only on $k,d$,

$\displaystyle\sum_{i\in[N]^{d}\setminus(I\cup[k+1]^{d})}\frac{1}{\xi_{i}^{2}}\leq
cn^{2s}/r^{d(2s-1)}.$

Plug this bound in (51) and in order to minimize the resulting bound, choose
$r$ to balance

$(r+k+1)^{d}\quad\text{with}\quad\frac{C_{n}}{\sqrt{n}}\sqrt{n^{2s}/r^{d(2s-1)}\log
n}.$

This leads us to take

$(r+k+1)^{d}\asymp(C_{n}\sqrt{\log n})^{\frac{2}{2s+1}}n^{\frac{2s-1}{2s+1}}$

when $C_{n}\sqrt{\log n}/n=O_{\mathbb{P}}(1)$ and $r=1$ otherwise. With this
choice (51) gives the desired bound for $s>1/2$. This completes the proof. ∎

###  A.5 Incoherence of Kronecker product type operators

Let

$\Delta=\left[\begin{array}[]{c}D\otimes I\otimes\cdots\otimes I\\\ I\otimes
D\otimes\cdots\otimes I\\\ \vdots\\\ I\otimes I\otimes\cdots\otimes
D\end{array}\right]$ (52)

where each Kronecker product has $d$ terms. With $D=D^{(k+1)}_{N,1}$, we get
the KTF penalty operator $\Delta=D_{n,d}^{(k+1)}$.

######  Lemma 1.

Let $\Delta$ be as defined in (52) for a matrix
$D\in\mathbb{R}^{N^{\prime}\times N}$ with $N^{\prime}\leq N$. Let
$\gamma_{i},u_{i},v_{i},i\in[N]$ denote the singular values of $D$, its left
and right singular vectors. Note that
$\gamma_{i}=0,u_{i}=0,v_{i}\in\mathrm{null}(D)$ for $i\in[p]$ where
$p=\mathrm{nullity}(D)$. If these singular vectors are incoherent, that is
$\|v_{i}\|_{\infty}\leq\mu/\sqrt{N},\|u_{i}\|_{\infty}\leq\mu/\sqrt{N^{\prime}}$
for a constant $\mu\geq 1$, then the left singular vectors $\nu$ of $\Delta$
are incoherent with a constant $\mu^{d}$, that is,
$\|\nu\|_{\infty}\leq\mu^{d}/\sqrt{N^{d-1}N^{\prime}}$.

Note that $p=k+1$ when $\Delta$ is the KTF penalty operator with
$D=D_{N,1}^{(k+1)}.$

###### Proof of Lemma 1.

Abbreviate $\rho_{i}=\gamma_{i}^{2}$ for $i\in[N]$. We are looking for a total
of $N^{d}-p^{d}$ eigenvectors for $\Delta\Delta^{\mathsf{T}}$. Assume for
exposition that $d=3$. For any $(i,j,k)\in[N]^{d}\setminus[p]^{d}$ (where
$\setminus$ is the set difference operator), the vectors

$\displaystyle\nu_{i,j,k}:=\frac{1}{\sqrt{\rho_{i}+\rho_{j}+\rho_{k}}}\left[\begin{array}[]{c}\gamma_{i}\cdot
u_{i}\otimes v_{j}\otimes v_{k}\\\ \gamma_{j}\cdot v_{i}\otimes u_{j}\otimes
v_{k}\\\ \gamma_{k}\cdot v_{i}\otimes v_{j}\otimes u_{k}\end{array}\right]$
(56)

are eigenvectors of $\Delta\Delta^{\mathsf{T}}$ as verified below.

$\displaystyle\Delta\Delta^{\mathsf{T}}\left[\begin{array}[]{c}\gamma_{i}\cdot
u_{i}\otimes v_{j}\otimes v_{k}\\\ \gamma_{j}\cdot v_{i}\otimes u_{j}\otimes
v_{k}\\\ \gamma_{k}\cdot v_{i}\otimes v_{j}\otimes u_{k}\end{array}\right]$
$\displaystyle=\Delta\;\left(\gamma_{i}^{2}+\gamma_{j}^{2}+\gamma_{k}^{2}\right)v_{i}\otimes
v_{j}\otimes v_{k}$ (60)
$\displaystyle=\left(\rho_{i}+\rho_{j}+\rho_{k}\right)\left[\begin{array}[]{c}\gamma_{i}\cdot
u_{i}\otimes v_{j}\otimes v_{k}\\\ \gamma_{j}\cdot v_{i}\otimes u_{j}\otimes
v_{k}\\\ \gamma_{k}\cdot v_{i}\otimes v_{j}\otimes u_{k}\end{array}\right]$
(64)

We see all $N^{d}-p^{d}$ eigenvectors of $\Delta\Delta^{\mathsf{T}}$ here.
Notice that $\|z_{i,j,k}\|_{2}=1$ and the incoherence is readily available
given that the left and right singular vectors of $D$ are incoherent.

For general $d$, these $N^{d}-p^{d}$ eigenvectors are given by

$\displaystyle\nu_{i_{1},i_{2},.,i_{d}}=\frac{1}{\sqrt{\sum_{j=1}^{d}\rho_{i_{j}}}}\left[\begin{array}[]{c}\gamma_{i_{1}}\cdot
u_{i_{1}}\otimes v_{i_{2}}\otimes\ldots v_{i_{d}}\\\ \gamma_{i_{2}}\cdot
v_{i_{1}}\otimes u_{i_{2}}\otimes\ldots v_{i_{d}}\\\ \vdots\\\
\gamma_{i_{d}}\cdot v_{i_{1}}\otimes v_{i_{2}}\otimes\ldots
u_{i_{d}}\end{array}\right]$ (69)

with eigenvalues $\sum_{j=1}^{d}\rho_{i_{j}}$ and are easily seen to be
incoherent. ∎

###  A.6 Upper bound for continuous KTV class

Recalling the continuous analog of KTF penalty from (12), define the class

$\mathrm{KTV}_{n,d}^{k}(C)=\bigg{\\{}f:\sum_{j=1}^{d}\sum_{x_{-j}}\mathrm{TV}\bigg{(}\frac{\partial^{k}f(\cdot,x_{-j})}{\partial
x_{j}^{k}}\bigg{)}\leq C\bigg{\\}}$

for $C>0$. If the true signal $\theta_{0}$ on the grid is an evaluation of a
function $f\in\mathrm{KTV}_{k}^{d}(C)$, the rates in Theorem 3 hold with
$C_{n}$ replaced by $C,$ due to the following result.

######  Lemma 2.

Let $C>0$ and let $d\geq 1,k\geq 0$ be integers. For all
$f\in\mathrm{KTV}_{n,d}^{k}(C),$ if $\theta_{f}\in\mathbb{R}^{n}$ is the
evaluation of $f$ on the grid points $Z_{n,d}$, then

$\|D_{n,d}^{(k+1)}\theta_{f}\|_{1}\leq c_{1}C$

for a constant $c_{1}$ that depends only on $k$ and $d$.

###### Proof of Lemma 2.

Let $f$ be an arbitrary function from $\mathrm{KTV}_{n,d}^{k}(C)$. Pick a
$j\in[N]$ and an $x_{-j}$ and consider the function
$\phi(\cdot)=f(\cdot,x_{-j})$ ($f$ with all but its $j$ argument fixed to
elements of $x_{-j}$ appropriately in order). From Theorem 1 in Mammen (1991)
and its proof, there exists a spline $\tilde{\phi}$ such that

$\displaystyle\tilde{\phi}(i/N)$ $\displaystyle=\phi(i/N),\quad i\in[N]$
$\displaystyle\mathrm{TV}(\tilde{\phi}^{(k)})$
$\displaystyle\leq\mathrm{TV}(\phi^{(k)})$

Let $t_{1},\dots,t_{L}$ be the knots of $\tilde{\phi}$, which are not
necessarily in the set of input points. Because it is a spline, $\tilde{\phi}$
can be written as the sum of a polynomial and a linear combination of $k$th
degree truncated power basis functions $g_{t}:x\mapsto(x-t)_{+}^{k}/k!$

$\tilde{\phi}(u)=p(u)+\sum_{\ell=1}^{L}\beta_{\ell}g_{t_{\ell}}(u),\quad
u\in[0,1]$

where $p$ is a polynomial of degree $\leq k$ and
$\beta_{\ell}\in\mathbb{R},\ell\in[L]$. Let
$D_{\mathrm{1d}}^{(k+1)}=D^{(k+1)}_{N,1}$. Now

$\displaystyle\bigg{\|}D_{\mathrm{1d}}^{(k+1)}\begin{bmatrix}\phi(1/N)\\\
\vdots\\\ \phi(N/N)\end{bmatrix}\bigg{\|}_{1}$
$\displaystyle=\bigg{\|}D_{\mathrm{1d}}^{(k+1)}\begin{bmatrix}\tilde{\phi}(1/N)\\\
\vdots\\\ \tilde{\phi}(N/N)\end{bmatrix}\bigg{\|}_{1}$
$\displaystyle=\bigg{\|}D_{\mathrm{1d}}^{(k+1)}\begin{bmatrix}p(1/N)\\\
\vdots\\\ p(N/N)\end{bmatrix}+\sum_{\ell=1}^{L}\beta_{\ell}\cdot
D_{\mathrm{1d}}^{(k+1)}\begin{bmatrix}g_{t_{\ell}}(1/N)\\\ \vdots\\\
g_{t_{\ell}}(N/N)\end{bmatrix}\bigg{\|}_{1}$
$\displaystyle\leq\sum_{\ell=1}^{L}|\beta_{\ell}|\|D_{\mathrm{1d}}^{(k+1)}G_{\ell}^{(k)}\|_{1}$
(70)

where the vector $G_{\ell}^{(k)}$ is the evaluation of $g_{t_{\ell}}$ on
$1/N,\dots,N/N$, that is $(G_{\ell}^{(k)})_{i}=g_{t_{\ell}}(i/N),i\in[N]$.
Here we used the fact that $D_{\mathrm{1d}}^{(k+1)}$ times the evaluations of
a polynomial at the input points $1/N,\dots,N/N$ is $0$.

The terms in (70) can be bound as follows. For $\ell\in[L]$, let
$i_{\ell}=\max_{i\in[N]}\\{i/N\leq t_{\ell}\\}$, that is, let $i_{\ell}/N$ be
the largest input point that is not greater the knot $t_{\ell}$. For any
vector $v\in\mathbb{R}^{N}$, and $i\in[N-k-1]$,
$(D_{\mathrm{1d}}^{(k+1)}v)_{i}=0$ if $(v_{i},\dots,v_{i+k})$ is the
evaluation of a polynomial at $i/N,\dots,(i+k+1)/N$. $g_{t_{\ell}}$ is a
polynomial on $[0,t_{\ell}]$ and on $[t_{\ell},1]$ for $\ell\in[L]$.
Therefore, $D_{\mathrm{1d}}^{(k+1)}G_{\ell}^{(k)}$ is nonzero in at most $k+1$
elements. Letting $A_{i}$ denote the $i$th row of matrix $A$, we can write

$\displaystyle\|D_{\mathrm{1d}}^{(k+1)}G_{\ell}^{(k)}\|_{1}$
$\displaystyle=\sum_{i=1}^{N-k-1}\left|\big{(}D_{\mathrm{1d}}^{(k+1)}\big{)}_{i}G_{\ell}^{(k)}\right|$
$\displaystyle=\sum_{i=(i_{\ell}-k)\vee
1}^{i_{\ell}}\left|\big{(}D_{\mathrm{1d}}^{(k+1)}\big{)}_{i}G_{\ell}^{(k)}\right|$
$\displaystyle\leq\sum_{i=(i_{\ell}-k)\vee
1}^{i_{\ell}}\left\|\big{(}D_{\mathrm{1d}}^{(k+1)}\big{)}_{i}\right\|_{1}g_{t_{\ell}}\left(\frac{i_{\ell}+k+1}{N}\right)$
$\displaystyle\leq\sum_{i=(i_{\ell}-k)\vee
1}^{i_{\ell}}\left\|\big{(}D_{\mathrm{1d}}^{(k+1)}\big{)}_{i}\right\|_{1}\left(\frac{k+1}{N}\right)^{k}\frac{1}{k!}$
$\displaystyle\leq(k+1)\cdot\left(\frac{k+1}{N}\right)^{k}\frac{1}{k!}\max_{i\in[N-k-1]}\left\|\big{(}D_{\mathrm{1d}}^{(k+1)}\big{)}_{i}\right\|_{1}$
$\displaystyle=(k+1)\cdot\left(\frac{k+1}{N}\right)^{k}\frac{1}{k!}\cdot
2^{k+1}N^{k}$ $\displaystyle=b_{k}$

where $b_{k}$ is a constant depending only on $k$. Plugging this upper bound
in (70),

$\displaystyle\bigg{\|}D_{\mathrm{1d}}^{(k+1)}\begin{bmatrix}\phi(1/N)\\\
\vdots\\\ \phi(N/N)\end{bmatrix}\bigg{\|}_{1}$ $\displaystyle\leq
b_{k}\sum_{\ell=1}^{L}|\beta_{\ell}|=b_{k}\mathrm{TV}(\tilde{\phi}^{(k)})\leq
b_{k}\mathrm{TV}(\phi^{(k)}).$

This means,

$\displaystyle\|D_{n,d}^{(k+1)}\theta_{f}\|_{1}$
$\displaystyle=\sum_{j=1}^{d}\sum_{x_{-j}}\bigg{\|}D_{\mathrm{1d}}^{(k+1)}\begin{bmatrix}f(1/N,x_{-j})\\\
\vdots\\\ f(N/N,,x_{-j})\end{bmatrix}\bigg{\|}_{1}$
$\displaystyle\leq\sum_{j=1}^{d}\sum_{x_{-j}}b_{k}\mathrm{TV}\bigg{(}\frac{\partial^{k}f(\cdot,x_{-j})}{\partial
x_{j}^{k}}\bigg{)}$ $\displaystyle\leq b_{k}\cdot C.$

This completes the proof. ∎

###  A.7 Proof of Theorem 4

Here and henceforth, we use the notation $B_{p}(r)=\\{x:\|x\|_{p}\leq r\\}$
for the $\ell_{p}$ ball of radius $r$, where $p,r>0$ (and the ambient
dimension will be determined based on the context).

######  Lemma 3 (Lemma 7 in Sadhanala et al. (2016)).

Let $\mathcal{T}(r)=\\{\theta\in\mathbb{R}^{n}:\|D\theta\|_{1}\leq r\\}$ for a
matrix $D$ and $r>0$. Recall that
$\|D\|_{1,\infty}=\max_{i\in[n]}\|D_{i}\|_{1}$ where $D_{i}$ is the $i$th