column of $D$. Then for any $r>0$, it holds that
$B_{1}(r/\|D\|_{1,\infty})\subseteq\mathcal{T}(r)$.

From Lemma 3 and the fact that $\|D_{n,d}^{(k+1)}\|_{1,\infty}=2^{k+1}d$

$B_{1}(r/(2^{k+1}d))\subseteq\mathcal{T}_{n,d}^{k}(r).$ (71)

for any $r>0$, and integers $d\geq 1,k\geq 0$.

To prove Theorem 4 we will use the following result from Birge and Massart
(2001), which gives a lower bound for the risk in a normal means problem, over
$\ell_{p}$ balls. We state the result in our notation.

######  Lemma 4 (Proposition 5 of Birge and Massart (2001)).

Assume i.i.d. observations $y_{i}\sim N(\theta_{0,i},\sigma^{2})$,
$i=1,\ldots,n$, and $n\geq 2$. Then the minimax risk over the $\ell_{p}$ ball
$B_{p}(r_{n})$, where $0<p<2$, satisfies

$n\cdot R(B_{p}(r_{n}))\geq
c\cdot\begin{cases}\displaystyle\sigma^{2-p}r_{n}^{p}\bigg{[}1+\log\bigg{(}\frac{\sigma^{p}n}{r_{n}^{p}}\bigg{)}\bigg{]}^{1-p/2}&\text{if
$\sigma\sqrt{\log n}\leq r_{n}\leq\sigma n^{1/p}/\sqrt{\rho_{p}}$}\\\
r_{n}^{2}&\text{if $r_{n}<\sigma\sqrt{\log n}$}\\\
\sigma^{2}n/\rho_{p}&\text{if $r_{n}>\sigma n^{1/p}/\sqrt{\rho}$}\end{cases}.$

Here $c>0$ is a universal constant, and $\rho_{p}>1.76$ is the unique solution
of $\rho_{p}\log\rho_{p}=2/p$.

###### Proof of Theorem 4.

It suffices to show that the minimax optimal risk
$R\big{(}\mathcal{T}_{n,d}^{k}(C_{n})\big{)}$ is lower bounded by the three
terms present in the statement’s lower bound separately:

$\displaystyle R\big{(}\mathcal{T}_{n,d}^{k}(C_{n})\big{)}$
$\displaystyle=\Omega\bigg{(}\frac{\kappa\sigma^{2}}{n}\bigg{)},$
$\displaystyle R\big{(}\mathcal{T}_{n,d}^{k}(C_{n})\big{)}$
$\displaystyle=\Omega\bigg{(}\frac{\sigma C_{n}}{n}\wedge\sigma^{2}\bigg{)},$
(72) $\displaystyle R\big{(}\mathcal{T}_{n,d}^{k}(C_{n})\big{)}$
$\displaystyle=\Omega\bigg{(}\bigg{(}\frac{C_{n}}{n}\bigg{)}^{\frac{2}{2s+1}}\sigma^{\frac{4s}{2s+1}}\wedge\sigma^{2}\bigg{)},$

where $\kappa=\mathrm{nullity}\big{(}D_{n,d}^{(k+1)}\big{)}=(k+1)^{d}$. First,
as the null space of $D_{n,d}^{(k+1)}$ has dimension $\kappa$, we get the
first lower bound:

$\displaystyle\inf_{\hat{\theta}}\sup_{\theta_{0}\in\mathcal{T}_{n,d}^{k}(C_{n})}\frac{1}{n}\mathbb{E}\|\hat{\theta}-\theta_{0}\|_{2}^{2}\;\geq\;\inf_{\hat{\theta}}\sup_{\theta_{0}\in\mathrm{null}\big{(}D_{n,d}^{(k+1)}\big{)}}\frac{1}{n}\mathbb{E}\|\hat{\theta}-\theta_{0}\|_{2}^{2}\;\geq\;\frac{\kappa\sigma^{2}}{n}.$

We get the second lower bound in (72) by using the $\ell_{1}$-ball embedding

$B_{1}\left(C_{n}/d_{\max}\right)\subset\mathcal{T}_{n,d}^{k}(C_{n})$

from (71) and then using Lemma 4. Finally, from Theorem 4 in Sadhanala et al.
(2017), it follows that

$R\big{(}\mathcal{C}_{n,d}^{k}(L_{n})\big{)}=\Omega\Big{(}\Big{(}\frac{\sigma^{2}}{n}\Big{)}^{\frac{2s}{2s+1}}L_{n}^{\frac{2}{2s+1}}\wedge\sigma^{2}\Big{)}$
(73)

with additional tracking for $\sigma^{2}.$ Taking $L_{n}=C_{n}/n^{1-s}$ and
applying Proposition 3 would then give the third lower bound in (72). This
completes the proof. ∎

###  A.8 Proof of Theorem 5 (Minimax linear rate)

We use the following shorthand for the risk of an estimator $\hat{\theta}$
over a class $\mathcal{K}$:

$\mathrm{Risk}(\hat{\theta})=\sup_{\theta_{0}\in\mathcal{K}}\frac{1}{n}\mathbb{E}\|\hat{\theta}-\theta_{0}\|_{2}^{2}.$

For a matrix $S\in\mathbb{R}^{n\times n}$ let $\mathrm{Risk}(S)$ also denote
the risk of the linear smoother $\hat{\theta}=Sy$.

###### Proof of Theorem 5.

For brevity, denote $D=D_{n,d}^{(k+1)}$ and let $S$ stand for a linear
smoother in the context of this proof. The minimax linear risk for the class
$\mathcal{T}_{n,d}^{k}(C_{n})$ is

$\displaystyle R_{L}(\mathcal{T}_{n,d}^{k}(C_{n}))$
$\displaystyle=\inf_{S\in\mathbb{R}^{n\times
n}}\sup_{\theta_{0}\in\mathcal{T}_{n,d}^{k}(C_{n})}\;\frac{1}{n}\mathbb{E}\|Sy-\theta_{0}\|_{2}^{2}$
$\displaystyle=\inf_{S}\sup_{\theta_{0}\in\mathcal{T}_{n,d}^{k}(C_{n})}\;\frac{1}{n}\mathbb{E}\|S(\theta_{0}+\epsilon)-\theta_{0}\|_{2}^{2}$
$\displaystyle=\frac{1}{n}\inf_{S}\sup_{\theta_{0}\in\mathcal{T}_{n,d}^{k}(C_{n})}\;\sigma^{2}\|S\|_{F}^{2}+\|(S-I)\theta_{0}\|_{2}^{2}$

where in the last line we used the assumption that $\epsilon_{i},i\in[n]$ are
i.i.d. with mean zero and variance $\sigma^{2}$ and used the notation
$\|A\|_{F}$ for the Frobenius norm of a matrix $A$. The infimum can be
restricted to the set of linear smoothers

$\mathbb{S}=\left\\{S:\mathrm{null}(S-I)\supseteq\mathrm{null}(D)\right\\}$

because if for a linear smoother $S$, if there exists
$\eta\in\mathrm{null}(D)$ such that $(S-I)\eta\neq 0$, then the inner supremum
above will be $\infty$, that is, its risk will be $\infty$. If the outer
infimum is over $\mathbb{S}$, then the supremum can be restricted to
$\\{\theta_{0}\in\mathrm{row}(D):\theta\in\mathcal{T}_{n,d}^{k}(C_{n})\\}$. We
continue to lower bound minimax linear risk as follows:

$\displaystyle R_{L}(\mathcal{T}_{n,d}^{k}(C_{n}))$
$\displaystyle=\frac{1}{n}\inf_{S\in\mathbb{S}}\sigma^{2}\|S\|_{F}^{2}+\sup_{\theta_{0}\in\mathrm{row}(D):\|D\theta_{0}\|_{1}\leq
C_{n}}\;\|(S-I)\theta_{0}\|_{2}^{2}$
$\displaystyle=\frac{1}{n}\inf_{S\in\mathbb{S}}\sigma^{2}\|S\|_{F}^{2}+\sup_{z:\|z\|_{1}\leq
C_{n}}\;\|(S-I)D^{+}z\|_{2}^{2}$
$\displaystyle=\frac{1}{n}\inf_{S\in\mathbb{S}}\sigma^{2}\|S\|_{F}^{2}+C_{n}^{2}\max_{i\in[m]}\left\|\left((S-I)D^{+}\right)_{i}\right\|_{2}^{2}$
(74)
$\displaystyle\geq\frac{1}{n}\inf_{S\in\mathbb{S}}\sigma^{2}\|S\|_{F}^{2}+\frac{C_{n}^{2}}{m}\sum_{i=1}^{m}\left\|\left((S-I)D^{+}\right)_{i}\right\|_{2}^{2}$
$\displaystyle=\inf_{S\in\mathbb{S}}\underbrace{\frac{\sigma^{2}}{n}\|S\|_{F}^{2}+\frac{C_{n}^{2}}{mn}\left\|(S-I)D^{+}\right\|_{F}^{2}}_{=:r(S)}$
(75)

In the third line, $(A)_{i}$ denotes the $i$th column of matrix $A$ and $m$
denotes the number of rows in $D$. In the fourth line, we used the fact that
the maximum of a set is at least as much as their average. In the last line —
within the context of this proof — we define the quantity $r(S)$ which is a
lower bound on the risk of a linear smoother $S\in\mathbb{S}$.

Notice that $r(\cdot)$ is a quadratic in the entries of $S$ and the constraint
$S\in\mathbb{S}$ translates to linear constraints on the entries of $S$.
Writing the KKT conditions, after some work, we see that $r(\cdot)$ is
minimized at

$\displaystyle S_{0}=a_{n}\left(\sigma^{2}L^{(k+1)}+a_{n}I\right)^{-1}$ (76)

where we denote $a_{n}=\frac{C_{n}^{2}}{m}$ and $L^{(k+1)}=D^{\mathsf{T}}D$
(the inverse is well defined because $a_{n}>0$). Further,
$S_{0}\in\mathbb{S}$. Therefore,

$\displaystyle R_{L}(\mathcal{T}_{n,d}^{k}(C_{n}))$ $\displaystyle\geq
r(S_{0}).$ (77)

We simplify the expression for $r(S_{0})$ now. Let $\lambda_{i},i\in[n]$ be
the eigenvalues of $L^{(k+1)}$. Then the eigenvalues of $S_{0}$ are

$\frac{a_{n}}{\sigma^{2}\xi_{i}+a_{n}},i\in[n]$

and the non-zero squared singular values of $(S_{0}-I)D^{+}$ are given by

$\frac{\sigma^{4}\xi_{i}}{(\sigma^{2}\xi_{i}+a_{n})^{2}},\quad\kappa<i\leq n.$

Using the fact that the squared Frobenius norm of a matrix is the sum of
squares of its singular values, substituting the above eigenvalues and
singular values in (75), we have

$\displaystyle r(S_{0})$
$\displaystyle=\frac{\sigma^{2}}{n}\sum_{i=1}^{n}\left(\frac{a_{n}}{\sigma^{2}\xi_{i}+a_{n}}\right)^{2}+\frac{a_{n}}{n}\sum_{i=1}^{n}\frac{\sigma^{4}\xi_{i}}{(\sigma^{2}\xi_{i}+a_{n})^{2}}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{\sigma^{2}a_{n}}{\sigma^{2}\xi_{i}+a_{n}}.$
(78)

Now we upper bound the risk $\mathrm{Risk}(S_{0})$ of the linear smoother
defined by $S_{0}$. From (74), we can write

$\displaystyle\mathrm{Risk}(S_{0})$
$\displaystyle=\frac{\sigma^{2}}{n}\|S_{0}\|_{F}^{2}+\frac{C_{n}^{2}}{n}\max_{i\in[m]}\left\|\left((S_{0}-I)D^{+}\right)_{i}\right\|_{2}^{2}.$

Let $D=U\Sigma V^{\mathsf{T}}$ be the singular value decomposition of $D$.
Also let the eigen-decomposition of $S_{0}-I=V\Lambda V^{\mathsf{T}}$. Then
using incoherence of columns of $U$, that is, the fact that there exists a
constant $c>1$ that depends only on $k,d$ such that
$U_{ij}^{2}\leq\frac{c}{m}$ for all $i\in[m],j\in[n]$, we can write

$\displaystyle\max_{i\in[m]}\left\|\left((S_{0}-I)D^{+}\right)_{i}\right\|_{2}^{2}$
$\displaystyle=\max_{i\in[m]}\left\|V\Lambda
V^{\mathsf{T}}V\Sigma^{+}(U^{\mathsf{T}})_{i}\right\|_{2}^{2}$
$\displaystyle=\max_{i\in[m]}\;(U^{\mathsf{T}})_{i}^{\mathsf{T}}(\Lambda\Sigma^{+})^{2}(U^{\mathsf{T}})_{i}$
$\displaystyle\leq\frac{c}{m}\mathrm{tr}\left((\Lambda\Sigma^{+})^{2}\right)$
$\displaystyle=\frac{c}{m}\sum_{i=1}^{n}\frac{\sigma^{4}\xi_{i}}{(\sigma^{2}\xi_{i}+a_{n})^{2}}.$

Plugging this back in the previous display and also using the fact that the
squared Frobenius norm of a matrix is equal to the sum of the squares of its
eigenvalues,

$\displaystyle\mathrm{Risk}(S_{0})$
$\displaystyle=\frac{\sigma^{2}}{n}\sum_{i=1}^{n}\left(\frac{a_{n}}{\sigma^{2}\xi_{i}+a_{n}}\right)^{2}+\frac{c\cdot
a_{n}}{n}\sum_{i=1}^{n}\frac{\sigma^{4}\xi_{i}}{(\sigma^{2}\xi_{i}+a_{n})^{2}}$
$\displaystyle\leq c\cdot r(S_{0})$

Combining this with the lower bound in (77), we have

$\displaystyle r(S_{0})\leq
R_{L}(\mathcal{T}_{n,d}^{k}(C_{n}))\leq\min\left\\{\sigma^{2},\mathrm{Risk}(S_{0})\right\\}\leq\min\left\\{\sigma^{2},\;\;c\cdot
r(S_{0})\right\\}.$ (79)

In other words, the minimax linear rate is essentially $r(S_{0})$ up to a
constant factor. Further, one of the estimators $\hat{y}=S_{0}y$, $\hat{y}=y$
achieves the minimax linear rate up to a constant factor.

Now we bound $r(S_{0}).$ Let $\kappa=(k+1)^{d}$ denote the nullity of $D$.
Recall from (78)

$\displaystyle
r(S_{0})=\frac{1}{n}\sum_{i=1}^{n}\frac{\sigma^{2}a_{n}}{\sigma^{2}\xi_{i}+a_{n}}=\frac{\kappa\sigma^{2}}{n}+\frac{1}{n}\sum_{i=\kappa+1}^{n}\frac{\sigma^{2}a_{n}}{\sigma^{2}\xi_{i}+a_{n}}.$
(80)

Lower bounding $r(S_{0})$. We give three lower bounds on $r(S_{0})$. By using
the fact that arithmetic mean of positive numbers is at least as large as
their harmonic mean, we have

$\displaystyle r(S_{0})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{\sigma^{2}a_{n}}{\sigma^{2}\xi_{i}+a_{n}}$
$\displaystyle\geq\frac{n\sigma^{2}a_{n}}{\sum_{i=1}^{n}(\sigma^{2}\xi_{i}+a_{n})}$
$\displaystyle=\frac{n\sigma^{2}a_{n}}{na_{n}+\sigma^{2}\|D\|_{F}^{2}}$
$\displaystyle=\frac{n\sigma^{2}a_{n}}{na_{n}+\sigma^{2}dn^{1-1/d}\|D_{\mathrm{1d}}^{(k+1)}\|_{F}^{2}}$
$\displaystyle=\frac{\sigma^{2}a_{n}}{a_{n}+\sigma^{2}dn^{-1/d}(n^{1/d}-k-1){{2k+2}\choose{k+1}}}$
$\displaystyle\geq\frac{\sigma^{2}a_{n}}{a_{n}+\sigma^{2}d4^{k+1}}$ (81)

Now we bound in $r(S_{0})$ in a second way. Let $n_{1}$ be the cardinality of
$\\{i\in[n]:\sigma^{2}\xi_{i}\leq a_{n}\\}$. Then

$\displaystyle
r(S_{0})=\frac{1}{n}\sum_{i=1}^{n}\frac{\sigma^{2}a_{n}}{\sigma^{2}\xi_{i}+a_{n}}\geq\frac{1}{n}\sum_{i=1}^{n_{1}}\frac{\sigma^{2}a_{n}}{a_{n}+a_{n}}=\frac{n_{1}\sigma^{2}}{2n}.$

Note that $n_{1}=\lfloor nF(a_{n}/\sigma^{2})\rfloor$ where $F$ is the
spectral distribution of $(D_{n,d}^{(k+1)})^{\mathsf{T}}D_{n,d}^{(k+1)}$
defined in Lemma 10. Applying Lemma 10, we get

$\displaystyle r(S_{0})$
$\displaystyle\geq\frac{\sigma^{2}}{2}\left(F\left(\frac{a_{n}}{\sigma^{2}}\right)-\frac{1}{n}\right)$
$\displaystyle\geq
c\sigma^{2}\min\big{\\{}1,(a_{n}/\sigma^{2})^{\frac{1}{2s}}\big{\\}}-\sigma^{2}/2n$
$\displaystyle=\min\big{\\{}c\sigma^{2},c\sigma^{2-\frac{1}{s}}a_{n}^{\frac{1}{2s}}-\sigma^{2}/2n\big{\\}}$
(82)

In the special case $s=1/2$, from Lemma 9 we get a third bound:

$\displaystyle r(S_{0})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{\sigma^{2}a_{n}}{\sigma^{2}\xi_{i}+a_{n}}\geq
c_{1}a_{n}\log\big{(}1+c_{2}/a_{n})$ (83)

where $c_{1},c_{2}$ constants that depend only on $k,d$.

From (80),(81), (82) and (83) we have the lower bound

$\displaystyle
r(S_{0})\geq\max\left\\{\frac{\kappa\sigma^{2}}{n},\frac{\sigma^{2}a_{n}}{a_{n}+\sigma^{2}d2^{2k+2}},\sigma^{2}\wedge
c\sigma^{2-\frac{1}{s}}a_{n}^{\frac{1}{2s}}-\frac{\sigma^{2}}{2n}\right\\}.$
(84)

and an additional lower bound of $c_{1}a_{n}\log\big{(}1+c_{2}/a_{n})$ in the
case $s=1/2$. Substituting $a_{n}=C_{n}^{2}/m$, using the assumption that
$C_{n}^{2}/n\leq 1$ and treating $k,d,\sigma$ as constants , we get the stated
lower bound.

Upper bounding $r(S_{0}).$ If $s<1/2$, then

$\displaystyle r(S_{0})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{\sigma^{2}a_{n}}{\sigma^{2}\xi_{i}+a_{n}}$
$\displaystyle\leq\frac{\kappa\sigma^{2}}{n}+\frac{1}{n}\sum_{i=\kappa+1}^{n}\frac{\sigma^{2}a_{n}}{\sigma^{2}\xi_{i}}$
$\displaystyle=\frac{\kappa\sigma^{2}}{n}+\frac{a_{n}}{n}\sum_{i=1}^{\kappa+1}\frac{1}{\xi_{i}}$
$\displaystyle\leq\frac{\kappa\sigma^{2}}{n}+\frac{a_{n}}{n}(c_{3}n)$
$\displaystyle=\frac{\kappa\sigma^{2}}{n}+c_{3}a_{n}$ (85)

We used Lemma 6 to control the second term in the third line. Similarly, if
$s=1/2$, $r(S_{0})\leq\kappa\sigma^{2}/n+c_{3}a_{n}\log n$. For the case
$s>1/2$, we can write

$\displaystyle r(S_{0})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{\sigma^{2}a_{n}}{\sigma^{2}\xi_{i}+a_{n}}$
$\displaystyle\leq\frac{1}{n}\sum_{i=1}^{n_{1}}\frac{\sigma^{2}a_{n}}{a_{n}}+\frac{1}{n}\sum_{i=n_{1}+1}^{n}\frac{\sigma^{2}a_{n}}{2\sigma^{2}\xi_{i}}$
$\displaystyle=\frac{n_{1}\sigma^{2}}{n}+\frac{a_{n}}{2n}\sum_{i=n_{1}+1}^{n}\frac{1}{\xi_{i}}$
$\displaystyle\leq
c\frac{\sigma^{2}}{n}+c\sigma^{2}\left(\frac{a_{n}}{\sigma^{2}}\right)^{\frac{1}{2s}}+c\frac{a_{n}}{2n}n^{2s}\left(n(a_{n}/\sigma^{2})^{\frac{1}{2s}}\right)^{1-2s}$
$\displaystyle\leq
c\frac{\sigma^{2}}{n}+c\sigma^{2-\frac{1}{s}}a_{n}^{\frac{1}{2s}}$ (86)

To get the fourth line, we used Lemma 10 to bound $n_{1}$ and Lemma 6 to bound
the summation.

#### Upper bound with the polynomial projection estimator
$\hat{\theta}^{\mathrm{poly}}$.

For brevity, let $\Pi$ denote the matrix that projects on to the null space of
$D$. Note that $(I-\Pi)D^{+}=D^{+}$. From bias variance decomposition similar
to that in (74),

$\displaystyle\frac{1}{n}\sup_{\theta_{0}\in\mathcal{T}_{n,d}^{k}(C_{n})}\mathbb{E}\big{[}\|\hat{\theta}^{\mathrm{poly}}-\theta_{0}\|_{2}^{2}\big{]}$
$\displaystyle=\frac{\sigma^{2}}{n}\|\Pi\|_{F}^{2}+\max_{i\in[m]}\|\big{(}(\Pi-I)D^{+}\big{)}_{i}\|_{2}^{2}$
$\displaystyle=\frac{\kappa\sigma^{2}}{n}+\max_{i\in[m]}\|D^{+}_{i}\|_{2}^{2}$

Then using incoherence of columns of $U$, that is, the fact that there exists
a constant $c>1$ that depends only on $k,d$ such that
$U_{ij}^{2}\leq\frac{c}{m}$ for all $i\in[m],j\in[n]$, we can write

$\displaystyle\max_{i\in[m]}\|D^{+}_{i}\|_{2}^{2}$
$\displaystyle=\max_{i\in[m]}\left\|V\Sigma^{+}(U^{\mathsf{T}})_{i}\right\|_{2}^{2}$
$\displaystyle=\max_{i\in[m]}\;(U^{\mathsf{T}})_{i}^{\mathsf{T}}(\Sigma^{+})^{2}(U^{\mathsf{T}})_{i}$
$\displaystyle\leq\frac{c}{m}\mathrm{tr}\left((\Sigma^{+})^{2}\right)$
$\displaystyle=\frac{c}{m}\sum_{i=\kappa+1}^{n}\frac{1}{\xi_{i}}$

Plugging this back in the above display and using the bound on
$\sum_{i=\kappa+1}^{n}\frac{1}{\xi_{i}}$ from Lemma 6, we get the desired
result.

#### Upper bound with the projection estimator (34) when $s>1/2$.

From (74), for the projection estimator $\hat{\theta}=S_{Q}y$ in (34),

$\mathrm{Risk}(\hat{\theta})=\frac{\sigma^{2}}{n}|Q|+\frac{C_{n}^{2}}{n}\max_{i\in[m]}\left\|\left((S_{Q}-I)D^{+}\right)_{i}\right\|_{2}^{2}.$
(87)

Set $Q=[\tau]^{d}$ for a $\tau$ to be chosen later from $(k+2,N]$. Also write
$S_{Q}-I=V\Lambda_{Q}V^{\mathsf{T}}$. Again using incoherence of columns of
$U$, we can write

$\displaystyle\max_{i\in[m]}\left\|\left((S_{Q}-I)D^{+}\right)_{i}\right\|_{2}^{2}$
$\displaystyle=\max_{i\in[m]}\left\|V\Lambda_{Q}V^{\mathsf{T}}V\Sigma^{+}(U^{\mathsf{T}})_{i}\right\|_{2}^{2}$
$\displaystyle=\max_{i\in[m]}\;(U^{\mathsf{T}})_{i}^{\mathsf{T}}(\Lambda_{Q}\Sigma^{+})^{2}(U^{\mathsf{T}})_{i}$
$\displaystyle\leq\frac{c}{m}\mathrm{tr}\left((\Lambda_{Q}\Sigma^{+})^{2}\right)$
$\displaystyle=\frac{c}{m}\sum_{i\in[N]^{d}\setminus Q}\frac{1}{\xi_{i}}$

The summation in the last line can be bound using Lemma 6 (recall $s>1/2$
here):

$\displaystyle\sum_{i\in[N]^{d}\setminus
Q}\frac{1}{\xi_{i}}\leq\sum_{\|(i-k-2)_{+}\|_{2}\geq\tau-k-2}\frac{1}{\xi_{i}}\leq
cn(n/(\tau-k-2)^{d})^{2s-1}$

Tracing this back to (87),

$\mathrm{Risk}(\hat{\theta})\leq\frac{\sigma^{2}}{n}\tau^{d}+\frac{cC_{n}^{2}}{m}\cdot(n/(\tau-k-2)^{d})^{2s-1}$

Minimize this bound by setting $\tau$ such that
$\tau^{d}\asymp(C_{n}/\sigma)^{\frac{1}{s}}n^{1-\frac{1}{2s}}$ to get the
desired bound. ∎

######  Remark 14.

In Theorem 5, when $s\leq 1/2$, the lower bound may also be obtained by
embedding the $\ell_{1}$-ball $B_{1}(C_{n}/(2^{k+1}d))$ into
$\mathcal{T}_{n,d}^{k}(C_{n})$.

###  A.9 Proof of Theorem 6 (Discrete Sobolev classes)

#### Proof of upper bound.

Like in the proof of minimax linear rates for KTV class in Theorem 5, for the
projection estimator $\hat{\theta}=S_{Q}y$ where
$S_{Q}=V_{Q}V_{Q}^{\mathsf{T}}$, we can derive

$\frac{1}{n}\sup_{\theta_{0}\in\mathcal{W}_{n,d}^{k+1}(B_{n})}\mathbb{E}\big{[}\|\hat{\theta}-\theta_{0}\|_{2}^{2}\big{]}=\frac{\sigma^{2}}{n}|Q|+\frac{1}{n}\sup_{\theta_{0}\in\mathcal{W}_{n,d}^{k+1}(B_{n})}\|(I-S_{Q})\theta_{0}\|_{2}^{2}.$

Denote $D=D_{n,d}^{(k+1)}$ for brevity. Set $Q=[\tau]^{d}$, where
$\tau\in(k+2,N]$ is an integer (recall $N=n^{1/d}$) and analyze the maximum of
the second term:

$\displaystyle\sup_{\theta_{0}:\|D\theta_{0}\|_{2}\leq
B_{n}}\frac{1}{n}\|(I-S_{q})\theta_{0}\|_{2}^{2}$
$\displaystyle=\sup_{z:\|z\|_{2}\leq
C_{n}}\frac{1}{n}\|(I-S_{q})D^{\dagger}z\|_{2}^{2}$
$\displaystyle=\frac{B_{n}^{2}}{n}\sigma^{2}_{\max}\big{(}(I-S_{q})D^{\dagger}\big{)}$
$\displaystyle\leq\frac{B_{n}^{2}}{n}\frac{1}{4^{k+1}\sin^{2k+2}(\pi(\tau-k-2)/(2N))}$
$\displaystyle\leq\frac{B_{n}^{2}}{n}\frac{N^{2k+2}}{(\pi(\tau-k-2))^{2k+2}}.$

Here we denote by $\sigma_{\max}(A)$ the maximum singular value of a matrix
$A$. The last inequality above used the inequality $\sin(x)\geq x/2$ for
$x\in[0,\pi/2]$. The earlier inequality used that
$\sigma^{2}_{\max}((I-S_{Q})D^{\dagger})$ is the reciprocal of the smallest
eigenvalue $\rho_{Q}$ of $M=D^{\mathsf{T}}D$ with index in $[N]^{d}\setminus
Q$. That is,

$\displaystyle\rho_{Q}=\rho_{\tau+1,1,\dots,1}\geq\big{(}4\sin^{2}(\pi(\tau-k-2)/(2N))\big{)}^{k+1},$

where the last inequality is due to the relation in (100). Hence, we have
established

$\sup_{\theta_{0}:\|D\theta_{0}\|_{2}\leq
B_{n}}\frac{1}{n}\mathbb{E}\big{[}\|\hat{\theta}-\theta_{0}\|_{2}^{2}\big{]}\leq\frac{\sigma^{2}}{n}\tau^{d}+\frac{B_{n}^{2}}{n}\frac{N^{2k+2}}{(\pi(\tau-k-2))^{2k+2}}.$

Choosing $\tau$ to balance the two terms on the right-hand side above results
in
$\tau^{d}\asymp(k+2)^{d}+\big{(}B_{n}^{2}n^{2}s/\sigma^{2}\big{)}^{\frac{1}{2s+1}}$.
Also, in the edge case where $Q=[N]^{d}$, the risk is $\sigma^{2}$. Plugging
this choice of $\tau$ gives the upper bound result.

#### Proof of lower bound.

Similar to argument in the proof of Theorem 4, the nullity of
$D_{n,d}^{(k+1)}$ implies the lower bound

$\displaystyle
R(\mathcal{W}_{n,d}^{k+1}(C_{n}))=\Omega\big{(}\frac{\kappa\sigma^{2}}{n}\big{)}.$
(88)

The Holder ball embedding

$\mathcal{W}_{n,d}^{k+1}(C_{n})\supseteq\mathcal{C}_{n,d}^{k}(cC_{n}n^{s-\frac{1}{2}})$

implies that

$R(\mathcal{W}_{n,d}^{k+1}(C_{n}))\geq
R\big{(}\mathcal{C}_{n,d}^{k}(cC_{n}n^{s-\frac{1}{2}})\big{)}\gtrsim\bigg{(}\frac{C_{n}^{2}}{n}\bigg{)}^{\frac{1}{2s+1}}\sigma^{\frac{4s}{2s+1}}\wedge\sigma^{2},$

where the second inequality follows from (73). Putting these two bounds
together, we get the desired lower bound.

##  Appendix B Estimation theory for Graph trend filtering on grids

We recall the GTF operator from Wang et al. (2016) for convenience. Let
$G(V,E)$ be a graph with $n$ vertices and $m$ edges
$(u_{1},v_{1}),\dots,(u_{m},v_{m})\in[n]\times[n]$. Let
$D\in\mathbb{R}^{m\times n}$ be the incidence matrix of $G$ satisfying

$|(Dx)_{j}|=|x_{u_{j}}-x_{v_{j}}|\quad\text{ for all }x\in\mathbb{R}^{n}$

for all edges $(u_{j},v_{j})$ for $j\in[m]$. The graph Laplacian is
$L=D^{\mathsf{T}}D$. The GTF operators of all orders are defined by

$\displaystyle S_{n,d}^{(1)}=D,\quad$ $\displaystyle S_{n,d}^{(2)}=L,$
$\displaystyle S_{n,d}^{(2k+1)}=DL^{k},\quad$ $\displaystyle
S_{n,d}^{(2k)}=L^{k}\text{ for }k\geq 0,k\in\mathbb{Z}.$ (89)

###  B.1 Upper bounds for GTF

Wang et al. (2016) used Theorem 2 (their Theorem 6) in order to derive error
rates for GTF on 2d grids already; see their Corollary 8. Sadhanala et al.
(2017) refine this result using a tighter upper bound for the partial sum of
inverse eigenvalues. Here, we give a more general result that applies to not
just 2d grids, but all $d\geq 2$ and $k\geq 0$. We further show that these
rates are optimal by deriving a matching lower bound. Recall the abbreviation
$s=(k+1)/d$.

######  Theorem 7.

Assume that $d\geq 1$ and $k\geq 0$. Denote
$C_{n}=\|S_{n,d}^{(k+1)}\theta_{0}\|_{1}$. Then GTF defined by the estimator
in (30) with $D=S_{n,d}^{(k+1)}$ in (89) satisfies

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

with

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

With canonical scaling of $C_{n}$, we see the following error bound.

######  Corollary 2.

With canonical scaling $C_{n}=C_{n}^{*}=n^{1-s}$, the GTF estimator with
$\lambda$ scaling as in Theorem 7 satisfies

$\displaystyle\sup_{\theta_{0}\in\mathcal{S}_{n,d}^{k}(C_{n})}\frac{1}{n}\|\hat{\theta}-\theta_{0}\|_{2}^{2}=\begin{cases}O_{\mathbb{P}}\left(n^{-s}\sqrt{\log
n}\right)&s<1/2\\\ O_{\mathbb{P}}\left(n^{-s}\log n\right)&s=1/2\\\
O_{\mathbb{P}}\left(n^{-\frac{2s}{2s+1}}(\log
n)^{\frac{1}{2s+1}}\right)&s>1/2.\end{cases}$

Remarks following Theorem 3 for KTF apply for GTF as well. The proof is in
Appendix B.4.

###  B.2 Lower bound

Similar to the lower bound in Theorem 4 for KTV class, we give a bound for the
Graph total variation (GTV) class

$\displaystyle\mathcal{S}_{n,d}^{k}(C_{n})$
$\displaystyle=\\{\theta\in\mathbb{R}^{n}:\|S_{n,d}^{(k+1)}\theta\|_{1}\leq
C_{n}\\}.$ (90)

Due to the lower order discrete derivatives on the boundary of the grid
$Z_{n,d}$, the GTV class $\mathcal{S}_{n,d}^{k}(C_{n})$ cannot contain the
discrete Holder class with appropriate scaling
$\mathcal{C}_{n,d}^{k}(C_{n}n^{s-1})$; see Lemma 4 in Sadhanala et al. (2017).
However, by an alternative route Sadhanala et al. (2017) show a lower bound
for $\mathcal{S}_{n,d}^{k}(C_{n})$ that matches with the lower bound for the
Holder class $\mathcal{C}_{n,d}^{k}(C_{n}n^{s-1})$. We further tighten their
result by embedding an $\ell_{1}$ ball of appropriate size.

######  Theorem 8.

For any integers $k\geq 0$, $d\geq 1$, the minimax estimation error for GTV
class defined in (90) satisfies

$\displaystyle R\big{(}\mathcal{S}_{n,d}^{k}(C_{n})\big{)}$
$\displaystyle=\Omega\bigg{(}\frac{\sigma^{2}}{n}+\frac{\sigma
C_{n}}{n}+\left(\frac{C_{n}}{n}\right)^{\frac{2}{2s+1}}\sigma^{\frac{4s}{2s+1}}\wedge\sigma^{2}\bigg{)}.$

###### Proof of Theorem 8.

Similar to the proof of Theorem 4, it is sufficient to show three lower bounds
separately. We get the first two lower bounds just as in the proof of Theorem
4 by using the fact that $\mathrm{nullity}\big{(}S_{n,d}^{(k+1)}\big{)}=1$ and
the $\ell_{1}$-ball embedding

$B_{1}(C_{n}/(2^{k+1}d))\subseteq\mathcal{S}_{n,d}^{k}(C_{n})$

from Lemma 3 and the fact that $\|S_{n,d}^{(k+1)}\|_{1,\infty}\leq 2^{k+1}d$.
The third term is from Theorem 5 in Sadhanala et al. (2017). ∎

###  B.3 Minimax linear rate for GTV class

The minimax linear rate analysis for GTV class is very similar to that for KTV
class. So we simply state the result and skip the proof.

######  Theorem 9.

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}$ (91)

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 mean estimator,
$\hat{\theta}^{\mathrm{mean}}=\bar{y}\mathds{1}$, achieves the rate in (91).
When $s=1/2$, this estimator achieves the rate in (91) up to a log factor.
Lastly, if $C_{n}^{2}=O(n^{\alpha})$ for $\alpha<1$, and still $s=1/2$, then
the mean estimator achieves the rate in (91) without the additional log
factor.

###  B.4 Proof of Theorem 7

For $d=2$, it is shown in the proof of Corollary 8 in Wang et al. (2016) that
the GTF operator $S_{n,d}^{(k+1)}$ satisfies the incoherence property, as
defined in Theorem 2, with a constant $\mu=4$ when $k$ is even and $\mu=2$
when $k$ is odd. Here we extend this incoherence property for $d>2$ using
Lemma 1. We treat the cases where $k$ is odd and even separately.

If $k$ is odd we can extend the argument from Corollary 8 in Wang et al.
(2016) in a straightforward manner. The GTF operator is
$S_{n,d}^{(k+1)}=L^{(k+1)/2}$ where $L$ is the Laplacian of the
$d$-dimensional grid graph. Denoting the Laplacian of the chain graph of
length $N$ by $L_{\mathrm{1d}}$, $L$ is given by

$L=L_{\mathrm{1d}}\otimes I\otimes I+I\otimes L_{\mathrm{1d}}\otimes
I+I\otimes I\otimes L_{\mathrm{1d}}$

for $d=3$ and

$L=L_{\mathrm{1d}}\otimes I\dots\otimes I+I\otimes L_{\mathrm{1d}}\dots\otimes
I+\dots+I\otimes\dots I\otimes L_{\mathrm{1d}}$

for general $d$ where each term in the summation is a Kronecker product of $d$
matrices. Let $\alpha_{i},u_{i},i\in[N]$ be the eigenvalues and eigenvectors
of $L_{\mathrm{1d}}$. As shown in Wang et al. (2016), in 1d, we have the
incoherence property $\|u_{i}\|_{\infty}\leq\sqrt{2/N}$ for all $i\in[N]$. The
eigenvalues of $L$ are $\sum_{j=1}^{d}\alpha_{i_{j}}$ and the corresponding
eigenvectors are $u_{i_{1}}\otimes\dots\otimes u_{i_{d}}$ for
$i_{1},\dots,i_{d}\in[N]$. Clearly, incoherence holds for the eigenvectors of
$L$ with constant $\mu=2^{d/2}$.

If $k$ is even, then the left singular vectors of $S_{n,d}^{(k+1)}$ are the
same as those of $S_{n,d}^{(1)}$. We know that both the left and right
singular vectors of $D_{\mathrm{1d}}^{(1)}$ satisfy the incoherence property
with constant $\mu=\sqrt{2}$ (see the proof of Corollary 7 in Wang et al.
(2016)). Setting $D=D_{\mathrm{1d}}^{(1)}$ in Lemma 1, we see that the left
singular vectors of $S_{n,d}^{(1)}$ and hence those of $S_{n,d}^{(k+1)}$
satisfy incoherence property with constant $2^{d/2}$. Therefore, for all
integers $k\geq 0$, the left singular vectors of $S_{n,d}^{(k+1)}$ are
incoherent with constant $2^{d/2}$.

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

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

where $\rho_{i},i\in[N]^{d}$ are the eigenvalues
${S_{n,d}^{(k+1)}}^{\mathsf{T}}S_{n,d}^{(k+1)}$ and $\mu=2^{d/2}$.

Consider the set $I=\\{i\in[N]^{d}:\|i-1\|_{2}<r\\}\setminus\\{1\\}^{d}$ for
an $r\in[1,\sqrt{d}N]$ chosen later. Lemma 5 gives the key calculation where
it is shown that for large enough $n$,

$\displaystyle\sum_{\|i-1\|\geq r}\frac{1}{\rho_{i}^{2}}=\sum_{\|i-1\|\geq
r}\frac{1}{\lambda_{i}^{k+1}}\leq c\begin{cases}n&s<1/2\\\
n\log(2\sqrt{d}N/r)&s=1/2\\\ n(n/r^{d})^{2s-1}&s>1/2\end{cases}$

where $\lambda_{i},i\in[N]^{d}$ are eigenvalues of the Laplacian $L$ and $c>0$
is a constant that depends only on $k,d$ .

For $s\leq 1/2$, to minimize the upper bound in (92), set $r=1$ so that $I$ is
empty and apply the above inequality. This gives the desired bound. Now
consider $s>1/2$. Note that $|I|\leq r^{d}$ because $I\subseteq[r]^{d}$.
Therefore (92) reduces to

$\frac{1}{n}\|\hat{\theta}-\theta_{0}\|_{2}^{2}=O_{\mathbb{P}}\Big{(}\frac{r^{d}}{n}+\frac{\mu}{n}\sqrt{\log
n(n/r^{d})^{2s-1}}\|\Delta\theta_{0}\|_{1}\Big{)}$ (93)

To minimize the upper bound in (93) balance

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

This leads us to take

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

and plugging this in (93) gives the desired bound for $s>1/2$. This completes
the proof. ∎

##  Appendix C Technical lemmas

######  Lemma 5.

Consider the eigenvalues $\\{\lambda_{i}:i=(i_{1},\cdots,i_{d})\in[N]^{d}\\}$
of the $d$-dimensional grid graph Laplacian with $n=N^{d}$ nodes. Let $k$ be a
non-negative integer and $r_{0}\in[1,\sqrt{d}N]$. Then,

$\sum_{i\in[N]^{d}:\|i-1\|_{2}^{2}\geq r_{0}^{2}}\frac{1}{\lambda_{i}^{k}}\leq
c\begin{cases}n&2k<d\\\ n\log(2\sqrt{d}N/r_{0})&2k=d\\\
N^{2k}r_{0}^{d-2k}&2k>d\end{cases}$

for a constant $c>0$ that depends on $k,d$ but not on $N,r_{0}$.

###### Proof of Lemma 5.

Let $I$ denote the summation on the left. Then

$\displaystyle I=\sum_{i\in[N]^{d}:\|i-1\|_{2}\geq
r_{0}}\frac{1}{\lambda_{i}^{k}}$ $\displaystyle=\sum_{\|i-1\|_{2}\geq
r_{0}}\Big{(}\sum_{j=1}^{d}4\sin^{2}\frac{\pi(i_{j}-1)}{2N}\Big{)}^{-k}$
$\displaystyle\leq\sum_{\|i-1\|_{2}\geq
r_{0}}\Big{(}\sum_{j=1}^{d}\frac{\pi^{2}(i_{j}-1)^{2}}{4N^{2}}\Big{)}^{-k}$
$\displaystyle=cN^{2k}\sum_{\|i-1\|_{2}\geq
r_{0}}\Big{(}\sum_{j=1}^{d}(i_{j}-1)^{2}\Big{)}^{-k}$ $\displaystyle\leq
cN^{2k}\sum_{i\in\\{0,1,..,N-1\\}^{d}:\|i\|_{2}\geq r_{0}}\|i\|_{2}^{-2k}$
(94)

In the second line, we used the fact that $\sin x\geq x/2$ for
$x\in[0,\pi/2]$.

#### Case $r_{0}\geq 2\sqrt{d}$.

In the last expression, upper bound $\|i\|_{2}^{-2k}$ with the integral of
$f(x)=\|x\|_{2}^{-2k},f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ over the unit
length cube whose top right corner is at $i$. Note that, the norm of any point
in this cube is at least
$\|i-\mathds{1}\|_{2}\geq\|i\|_{2}-\|\mathds{1}\|_{2}=\|i\|_{2}-\sqrt{d}\geq
r_{0}-\sqrt{d}\geq r_{0}/2.$ Therefore, we can continue to bound

$\displaystyle I$ $\displaystyle\leq
cN^{2k}\int_{r_{0}/2\leq\|x\|_{2}\leq\sqrt{d}N}\|x\|_{2}^{-2k}\;dx$
$\displaystyle\leq cN^{2k}\int_{r_{0}/2\leq
r\leq\sqrt{d}N}(r^{2})^{-k}r^{d-1}\;dr$

The last line is obtained by changing to polar coordinates and integrating out
the angles. Recall that the constants $c$ may change from line to line and
they may depend on $k,d$, but not on $N,r_{0}$.

If $d=2k,$ then the integral

$I\leq cN^{2k}\log(2\sqrt{d}N/r_{0})=cn\log(2\sqrt{d}N/r_{0}).$

If $2k<d,$ then

$I\leq cN^{2k}\big{(}(N\sqrt{d})^{d-2k}-(r_{0}/2)^{d-2k}\big{)}\leq cN^{d}.$

If $2k>d,$ then

$I\leq cN^{2k}\big{(}(r_{0}/2)^{d-2k}-(N\sqrt{d})^{d-2k}\big{)}.$

Treating $d,k$ as constants, we write

$I\leq cN^{2k}r_{0}^{d-2k}.$

#### Case $r_{0}<2\sqrt{d}$.

Continuing from (94), write

$\displaystyle I\leq
cN^{2k}\sum_{i\in\\{0,1,..,N-1\\}^{d}:\|i\|_{2}\in[r_{0},2\sqrt{d})}\|i\|_{2}^{-2k}+cN^{2k}\sum_{i\in\\{0,1,..,N-1\\}^{d}:\|i\|_{2}\geq
2\sqrt{d}}\|i\|_{2}^{-2k}$ (95)

From the previous case, the second summation can be upper bound with $cn$ if
$2k<d$, $cn\log n$ if $2k=d$ and $cN^{2k}$ if $2k>d$. In the first summation
(in the above display), the number of entries $i$ is at most $(2\sqrt{d})^{d}$
and each entry is at most $r_{0}^{-2k}$. Therefore the first term is at most
$cN^{2k}(2\sqrt{d})^{d}r_{0}^{-2k}$. Putting the two sums together, we can
verify the stated bounds. ∎

Following lemma provides a result analogous to Lemma 5 for KTF.

######  Lemma 6.

Let $\\{\xi_{i}:i=(i_{1},\cdots,i_{d})\in[N]^{d}\\}$ be the eigenvalues of
${D_{n,d}^{(k+1)}}^{\mathsf{T}}D_{n,d}^{(k+1)}$ and suppose
$r_{0}\in[1,\sqrt{d}N]$. Then,

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

In the case $2k+2>d$,

$\sum_{i\in[N]^{d}:\|(i-k-2)_{+}\|_{2}\geq r_{0}}\frac{1}{\xi_{i}^{2}}\leq
cN^{2k+2}r_{0}^{d-2k-2}$

Here $c>0$ is a constant that depends on $k,d$ but not on $N,r_{0}$.

###### Proof.

Using Lemma 7, we can write

$\sum_{i\in[N]^{d}:\|(i-k-2)_{+}\|_{2}\geq r_{0}}\frac{1}{\xi_{i}^{2}}\leq
d^{2k}\sum_{i\in[N]^{d}:\|(i-k-2)_{+}\|_{2}\geq
r_{0}}\frac{1}{\lambda_{i-k-1}^{2k+2}}\leq
d^{2k}\sum_{i\in[N]^{d}:\|i-1\|_{2}\geq r_{0}}\frac{1}{\lambda_{i}^{2k+2}}.$
(96)

Applying Lemma 5 directly gives the desired result in the case $2k+2>d$. In
the case $2k+2\leq d$, we get the bound by setting $r_{0}=1$ in (96) and then
applying Lemma 5. ∎

######  Lemma 7.

Let $\\{\xi_{i}:i=(i_{1},\cdots,i_{d})\in[N]^{d}\\}$ be the eigenvalues of
${D_{n,d}^{(k+1)}}^{\mathsf{T}}D_{n,d}^{(k+1)}$ for $k\geq 0,d\geq 1,N\geq
1,n=N^{d}$. Let $\alpha_{i},i\in[N]$ be the eigenvalues of $L$, the Laplacian
of chain graph of length $N$. Let
$\lambda_{i_{1},\dots,i_{d}}=\sum_{j=1}^{d}\alpha_{i_{j}},\quad i\leq N\text{
elementwise}$ with the convention that $\alpha_{\ell}=0$ for $\ell\leq 0$.
Then

$\xi_{i}\geq d^{-k}\lambda_{i-k-1}^{k+1}\text{ for }i\in[N]^{d}.$

###### Proof.

Abbreviate $D=D_{N,1}^{(k+1)}$, and let $G$ be the $k$th order GTF operator
defined over a 1d chain of length $N$. Also let $N^{\prime}=N-k-1$, and
$k^{\prime}=\lfloor(k+1)/2\rfloor$. Let

  * •

$\beta_{\ell},\ell\in[N^{\prime}]$ be the eigenvalues of $DD^{\mathsf{T}}$

  * •

$\gamma_{\ell}$, $\ell\in[N^{\prime\prime}]$ be the eigenvalues of
$GG^{\mathsf{T}}$ where $N^{\prime\prime}=N-1\\{k\text{ is even}\\}$

$GG^{\mathsf{T}}$ and $G^{\mathsf{T}}G$ should have the same nonzero
eigenvalues. From the definition of $G$, $G^{\mathsf{T}}G=L^{k+1}$. The first
eigenvalue of $L$ is 0 and the rest are nonzero. Putting these facts together,
we see that

$\gamma_{\ell}=\alpha_{\ell+N-N^{\prime\prime}}^{k+1}\text{ and
}\alpha_{\ell}^{k+1}\leq\gamma_{\ell}\quad\text{ for
}\ell\in[N^{\prime\prime}].$ (97)

Removing the top $k^{\prime}$ and bottom $k^{\prime}$ rows of $G$ yields $D$,
i.e.,

$D=PG,\quad\text{where}\;\,P=\left[\begin{array}[]{ccc}0_{N^{\prime}\times
k^{\prime}}&I_{N^{\prime}}&0_{N^{\prime}\times k^{\prime}}\end{array}\right].$

As $DD^{\mathsf{T}}=PGG^{\mathsf{T}}P^{\mathsf{T}}$ and
$PP^{\mathsf{T}}=I_{N^{\prime}}$, Cauchy interlacing theorem (Lemma 8) tells
us that

$\gamma_{i}\leq\beta_{i}\leq\gamma_{i+N^{\prime\prime}-N^{\prime}},\quad\text{
for }i\in[N^{\prime}].$ (98)

Thanks to the Kronecker sum structure, the eigenvalues of
$(D_{n,d}^{(k+1)})^{\mathsf{T}}D_{n,d}^{(k+1)}$ are

$\xi_{i_{1},\dots,i_{d}}=\sum_{j=1}^{d}\rho_{i_{j}},\quad i\in[N]^{d},$

where $\rho_{1},\dots,\rho_{N}$ denote the eigenvalues of $D^{\mathsf{T}}D$,
i.e., $\rho_{1}=\cdots=\rho_{k+1}=0$ and
$\rho_{\ell+k+1}=\beta_{\ell}\ell\in[N^{\prime}]$. Similarly, we can write the
eigenvalues of the Laplacian of the $d$-dimensional grid graph as

$\lambda_{i_{1},\dots,i_{d}}=\sum_{j=1}^{d}\alpha_{i_{j}},\quad i\in[N]^{d}.$

For arbitrary $i\in[N]^{d}$, we can write

$\xi_{i_{1},\dots,i_{d}}=\sum_{j=1}^{d}\beta_{i_{j}-k-1}\geq\sum_{j=1}^{d}\gamma_{i_{j}-k-1}\geq\sum_{j=1}^{d}\alpha_{i_{j}-k-1}^{k+1}\geq
d^{-k}\lambda_{i_{1}-k-1,\dots,i_{d}-k-1}^{k+1},$

with the convention $\alpha_{\ell}=\beta_{\ell}=\gamma_{\ell}=0$ for $\ell\leq
0$. The first inequality is due to (98), the second is due to (97), and the
third is due to a simple application of Jensen’s inequality:
$(\frac{1}{d}\sum_{j=1}^{d}a_{i})^{k}\leq\frac{1}{d}\sum_{j=1}^{d}a_{i}^{k}$
if $k\geq 1$ and $a\geq 0$ elementwise. ∎

######  Lemma 8 (Cauchy Interlacing theorem).

Let $A$ be an $n\times n$ symmetric matrix, $P\in\mathbb{R}^{m\times n}$ be an
orthogonal projection matrix (satisfying $PP^{\mathsf{T}}=I_{m}$) with $m\leq
n$ and define $B=PAP^{\mathsf{T}}$. Let
$\alpha_{1}\leq\alpha_{2}\leq\dots\leq\alpha_{n}$ be the eigenvalues of $A$
and $\beta_{1}\leq\beta_{2}\leq\dots\leq\beta_{m}$ be the eigenvalues of $B$.
Then

$\alpha_{i}\leq\beta_{i}\leq\alpha_{i+n-m},\quad\text{ for }i\in[m].$

######  Lemma 9.

Let $\\{\xi_{i}:i=(i_{1},\cdots,i_{d})\in[N]^{d}\\}$ be the eigenvalues of
${D_{n,d}^{(k+1)}}^{\mathsf{T}}D_{n,d}^{(k+1)}$ for $k\geq 0,d\geq 1,N\geq
1,n=N^{d}$. Suppose $s=1/2$ and $a>0$. Then

$\sum_{i\in[N]^{d}}\frac{1}{\xi_{i}+a}\geq
cn\log\big{(}1+\pi^{2k+2}a^{-1}\big{)}$

for a constant $c$ that depends only on $k,d$.

###### Proof of Lemma 9.

From (100) and the inequality $\sin x\leq x$ for $x\geq 0$, for any
$i\in[N]^{d}$,

$\displaystyle\xi_{i}=\sum_{j=1}^{d}\rho_{i_{j}}\leq\sum_{j=1}^{d}4^{k+1}\sin^{2k+2}\frac{\pi(i_{j}-1)}{2N}\leq\pi^{2k+2}n^{-2s}\|i-1\|_{2k+2}^{2k+2}\leq\pi^{2k+2}n^{-2s}\|i-1\|_{2}^{2k+2}.$

With this inequality,

$\displaystyle\sum_{i\in[N]^{d}}\frac{1}{\xi_{i}+a}$
$\displaystyle\geq\sum_{i\in[N]^{d}}\frac{1}{\pi^{2k+2}n^{-2s}\|i-1\|_{2}^{2k+2}+a}$
$\displaystyle\geq
c\int_{r=0}^{N}\frac{1}{\pi^{2k+2}n^{-2s}r^{2k+2}+a}r^{d-1}\;dr.$ (99)

In the second inequality is obtained as follows. Consider axis-parallel unit
cubes with corners located at integer coordinates. Let
$A_{i}\subset\mathbb{R}^{d}$ be the cube with its farthest corner from origin
located at $i$, for $i\in[N]^{d}$. Clearly,

$\frac{1}{\pi^{2k+2}n^{-2s}\|i-1\|_{2}^{2k+2}+a}\geq\int_{A_{i}}\frac{1}{\pi^{2k+2}n^{-2s}\|x\|_{2}^{2k+2}+a}\;dx.$

Next observe that the set $\\{x\in\mathbb{R}^{d}:\|x\|_{2}\leq N,x\geq 0\\}$
is contained in the cube $\\{x\in\mathbb{R}^{d}:\|x\|_{\infty}\leq N,x\geq
0\\}$ The former set is the non-negative orthant of the $\ell_{2}$ ball of
radius $N$ in $\mathbb{R}^{d}$. So, for radially symmetric functions $f$,
integral of $f$ over this set is $2^{-d}$ times its integral over the
$\ell_{2}$ ball. This justifies (99) after a change to polar coordinates. In
the integral (99), noting that $s=1/2,2k+2=d$, put $u=r^{d}$ to get

$\sum_{i\in[N]^{d}}\frac{1}{\xi_{i}+a}\geq
c\int_{0}^{n}\frac{1}{\pi^{2k+2}u/n+a}\;du=cn\log\big{(}1+\pi^{2k+2}a^{-1}\big{)}.\qed$

######  Lemma 10.

Let $\\{\xi_{i}:i=(i_{1},\cdots,i_{d})\in[N]^{d}\\}$ be the eigenvalues of
${D_{n,d}^{(k+1)}}^{\mathsf{T}}D_{n,d}^{(k+1)}$ for $k\geq 0,d\geq 1,N\geq
1,n=N^{d}$. Let $\alpha_{i},i\in[N]$ be the eigenvalues of $L$, the Laplacian
of chain graph of length $N$. Define

$F(t)=\frac{1}{n}\sum_{i\in[N]^{d}}1\\{\lambda_{i}\leq t\\},\quad\text{ for
}t\in[0,\lambda_{n}].$

Then there exist constants $c_{1},c_{2},c_{3}>0$ independent of $n,t$ such
that

$\displaystyle c_{1}t^{\frac{d}{2k+2}}\leq F(t)\leq
c_{2}+c_{3}t^{\frac{d}{2k+2}}$

for all $t\in[0,\lambda_{n}]$.

###### Proof of Lemma 10.

Using the notation in the proof of Lemma 6,

$\lambda_{i_{1},\dots,i_{d}}=\rho_{i_{1}}+\dots+\rho_{i_{d}},\quad\text{ for
}(i_{1},\dots,i_{d})\in[N]^{d}$

where $\rho_{\ell}=\beta_{\ell-k-1},\ell\in[N]$ with the convention that
$\beta_{\ell}=0$ for $\ell\leq 0$. From (98), (97) and the fact that the
eigenvalues of chain Laplacian are given by $4\sin^{2}\frac{\pi(\ell-1)}{2N}$
for $\ell\in[N]$, we have

$\displaystyle\left(4\sin^{2}\frac{\pi(\ell-k-2)_{+}}{2N}\right)^{k+1}\leq\rho_{i}\leq\left(4\sin^{2}\frac{\pi(\ell-1)}{2N}\right)^{k+1},\quad\text{
for }\ell\in[N]$ (100)

where $(x)_{+}=\max\\{x,0\\}$ for $x\in\mathbb{R}$. The upper bound can be
argued as follows.

$\displaystyle nF(t)$
$\displaystyle=\sum_{i\in[N]^{d}}1\Big{\\{}\lambda_{i_{1},\dots,i_{d}}\leq
t\Big{\\}}$
$\displaystyle=\sum_{i\in[N]^{d}}1\Big{\\{}\sum_{j=1}^{d}\rho_{i_{j}}\leq
t\Big{\\}}$
$\displaystyle\leq\sum_{i\in[N]^{d}}1\Big{\\{}\sum_{j=1}^{d}4^{k+1}\sin^{2k+2}\frac{\pi(i_{j}-k-2)_{+}}{2N}\leq
t\Big{\\}}$
$\displaystyle\leq\sum_{i\in[N]^{d}}1\Big{\\{}\sum_{j=1}^{d}\left(\frac{\pi}{2}\right)^{2k+2}(i_{j}-k-2)_{+}^{2k+2}\leq
tN^{2k+2}\Big{\\}}$

In the third line, we use (100) and in the fourth line, we use the fact that
$\sin x\geq x/2$ for $x\in[0,\pi/2]$. Observe that

$\displaystyle\left(\frac{\pi}{2}\right)^{2k+2}\sum_{j=1}^{d}(i_{j}-k-2)_{+}^{2k+2}\leq
tN^{2k+2}$
$\displaystyle\Rightarrow\|(i-k-2)_{+}\|_{\infty}\leq\frac{2}{\pi}Nt^{\frac{1}{2k+2}}$
$\displaystyle\Rightarrow\|i\|_{\infty}\leq
k+2+\frac{2}{\pi}Nt^{\frac{1}{2k+2}}$

Applying this fact to the previous bound on $nF(t)$,

$\displaystyle nF(t)\leq\sum_{i\in[N]^{d}}1\Big{\\{}\|i\|_{\infty}\leq
k+2+\frac{2}{\pi}Nt^{\frac{1}{2k+2}}\Big{\\}}\leq\Big{(}k+2+\frac{2}{\pi}Nt^{\frac{1}{2k+2}}\Big{)}^{d}\leq
2^{d-1}(k+2)^{d}+2^{2d-1}\pi^{-d}nt^{\frac{d}{2k+2}}.$

where in the last inequality we used the fact that $(a+b)^{d}\leq
2^{d-1}(a^{d}+b^{d})$ for $a,b\geq 0,d\geq 1$.

The lower bound can be derived as follows. Certainly, $F(t)\geq F(0)=\kappa/n$
for $t\geq 0$. We can write

$\displaystyle nF(t)$
$\displaystyle=\sum_{i\in[N]^{d}}1\big{\\{}\lambda_{i_{1},\dots,i_{d}}\leq
t\big{\\}}$
$\displaystyle=\sum_{i\in[N]^{d}}1\Big{\\{}\sum_{j=1}^{d}\rho_{i_{j}}\leq
t\Big{\\}}$
$\displaystyle=\sum_{i\in[N]^{d}}1\Big{\\{}\sum_{j=1}^{d}4^{k+1}\sin^{2k+2}\frac{\pi(i_{j}-1)}{2N}\leq
t\Big{\\}}$
$\displaystyle\geq\sum_{i\in[N]^{d}}1\Big{\\{}\sum_{j=1}^{d}\pi^{2k+2}(i_{j}-1)^{2k+2}\leq
tN^{2k+2}\Big{\\}}$
$\displaystyle=\sum_{i\in[N]^{d}}1\Big{\\{}\sum_{j=1}^{d}\|i-1\|_{2k+2}\leq
r\Big{\\}}$

where $r=\frac{1}{\pi}Nt^{\frac{1}{2k+2}}$. In the third line, we used (100)
and in the fourth line, we used the fact that $\sin x\leq x$ for $x\geq 0$.
Note that, we can inscribe a cube $\\{i:\|i-1\|_{\infty}\leq
rd^{\frac{-1}{2k+2}}\\}$ in the $\ell_{2k+2}$ body $\\{i:\|i-1\|_{2k+2}\leq
r\\}$ and the cube contains $\big{(}1+\lfloor
rd^{\frac{-1}{2k+2}}\rfloor\big{)}^{d}$ lattice points in $[N]^{d}$.
Therefore, continuing to bound from the previous display,

$nF(t)\geq\big{(}1+\lfloor\frac{1}{\pi}d^{\frac{-1}{2k+2}}Nt^{\frac{1}{2k+2}}\rfloor\big{)}^{d}\geq\frac{1}{\pi^{d}}d^{\frac{-d}{2k+2}}nt^{\frac{d}{2k+2}}.$

where in the last inequality we used the fact $(1+\lfloor x\rfloor)^{d}\geq
x^{d}$ for $x\geq 0$. ∎

##  Appendix D Fast algorithm for degrees of freedom

From (16), given a KTF estimate $\hat{\theta}$, $\mathrm{nullity}(D_{-A})$ is
an unbiased estimate of degrees of freedom of KTF where $A$ denotes the set of
rows $r$ in $D$ for which $(D\hat{\theta})_{r}\neq 0.$ We give an algorithm to
compute $\mathrm{nullity}(D_{-A})$ in $O(ndk)$ time.

Input: fitted values $\theta\in\mathbb{R}^{n}$, trend filtering order $k\geq
0$
Output: estimate of degrees of freedom

  1. 1.

struct Piece: set = false, knowns = 0, start, end $\in[N]^{d}$

  2. 2.

pieces : Piece[], pieces-containing[i]: Piece[], set[i] : bool for
$i\in[N]^{d}$

  3. 3.

for $d^{\prime}\in[d]$ for $i^{\prime}\in[N]^{d-1}$:

    1. (a)

$i\mathord{\leftarrow}i^{\prime}$ with an extra $1$ inserted before
$d^{\prime}$th entry:
$i_{j}\mathord{=}i^{\prime}_{j}$ for $j\mathord{<}d^{\prime}$,
$i_{j}\mathord{=}1$ for $j\mathord{=}d^{\prime}$,
$i_{j}\mathord{=}i^{\prime}_{j\mathord{-}1}$ for $j\mathord{>}d^{\prime}$

    2. (b)

Make-Polynomial-Pieces$(\theta,i,d^{\prime})$

  4. 4.

df = 0

  5. 5.

for $p$ in pieces:

    1. (a)

if $p.set:$ continue

    2. (b)

df
$\mathord{+}\mathord{=}\max\\{0,\min\\{p.length,k\mathord{+}1\\}-p.\mathrm{knowns}\\}$

    3. (c)

Spread$(p)$

  6. f.

return df

Spread-Vertex$(i)$

for piece $q$ in pieces-containing$[i]:$

  1. 1.

if $q.set$ continue

  2. 2.

$q.\mathrm{knowns}\mathord{+}\mathord{+}$

  3. 3.

if $q.\mathrm{knowns}\mathord{>}k:$ Spread$(q)$

Algorithm 3 Degrees-of-freedom$(\theta,k)$

Spread$(p:\text{Piece})$

$p.set$ = true

for vertex $i$ on the line $[p.start,p.end]$:

  1. 1.

if set$[i]:$ continue

  2. 2.

set$[i]\leftarrow$ true

  3. 3.

Spread-Vertex($i$)

Make-Polynomial-Pieces$(\theta,i\in[N]^{d},d^{\prime}\in[d])$ makes polynomial
pieces on the line containing $i$ along axis $d^{\prime}$

  1. 1.

$a_{j}\leftarrow\theta[i$ with $i_{d^{\prime}}=j]$ for $j\in[N]$

  2. 2.

while $j\leq N$

    1. (a)

start = $j$, end = $j$

    2. (b)

while ($\smash{j\mathord{\leq}N^{\prime}}$ and $\smash{\langle
w,a[j\mathord{:}j\mathord{+}k\mathord{+}1]\rangle\mathord{=}0})$
$j\mathord{+}\mathord{+}$

    3. (c)

if $j\mathord{\neq}$ start: end $\leftarrow j\mathord{+}k$

    4. (d)

Piece p(i with $i_{d^{\prime}}$= start, i with $i_{d^{\prime}}$ = end)

    5. (e)

for $j^{\prime}$ in [start, end]: pieces-containing$[i$ with
$i_{d^{\prime}}=j^{\prime}]$.add(p)

    6. (f)

pieces.add(p)

Notation: $N^{\prime}=N-k-1,$ $w\in\mathbb{R}^{k+2}$ is the $(k+1)$th order
difference vector.

###  D.1 Time complexity and correctness of the algorithm

Let $A$ denote the set of rows $r$ in $D$ for which $(D\hat{\theta})_{r}\neq
0$. Denote the null space of $D_{-A}$ with $\mathcal{N}$.

In step 3 of Degrees-of-freedom, we find line segments on the lattice where
$\theta$ is a $k$th degree polynomial. In a bit more detail, for each straight
line in the lattice between opposing faces, we find segments along the line
where $\theta$ is a $k$th degree polynomial. We call these line segments
(polynomial) _pieces_ in the algorithm. In 2d, Make-Polynomial-Pieces is
called on rows and columns separately, and a piece is a part of a row or a
column. An important characterization of the null space $\mathcal{N}$ is the
following:

A $\theta\in\mathcal{N}$ iff it is a $k$th degree polynomial on all the pieces
found in step 3.

If a piece has fewer than $k+2$ elements, then any $\theta$ is trivially a
$k$th degree polynomial on the piece. Otherwise, $k+1$ values on the piece
determine a polynomial on the piece.

We pretend to build a vector $\eta\in\mathcal{N}$ by making sure that $\eta$
is a $k$th degree polynomial on the pieces. In step 5, we pretend to set the
values of $\eta$ in a piece $p$. The number of new entries in $\eta$ required
to determine a polynomial on piece $p$ is shown in 5(b). Once the new entries
are picked arbitrarily, all the values on the piece are determined via the
constraints in $D_{-A}\eta=0.$ Then we propagate the values from this piece to
other adjoining pieces in a depth-first fashion. By the end of the procedure,
df accumulates the total number of free parameters that we can use to build
such a $\eta$. The dimension of $\mathcal{N}$ is equal to the number of free
parameters in the algorithm.

#### Time Complexity.

It takes $O(nd(k+1))$ time to make the polynomial pieces in line 3 of Degrees-
of-freedom. The number of pieces is at most $nd(k+1).$ Therefore the for loop
in line 5 is run at most as many times. Spread is called on a piece exactly
once and Spread-Vertex is called on a node exactly once. A node is contained
in a maximum of $(k+2)d$ pieces. Therefore, the total time complexity is
$O(nd(k+2)).$

#### Correctness.

Suppose the number of free parameters returned by the algorithm is $f$. Given
the values at the $f$ free nodes $F\subset[n]$, the values are determined at
all $n$ nodes. Further this mapping from $\mathbb{R}^{f}\mapsto\mathbb{R}^{n}$
is linear. Therefore there exists a matrix $C$ with size $n\times f$ such that
$Cb\in\mathcal{N}$ for any $b\in\mathbb{R}^{f}.$ Further, $(Cb)_{F}$ is a
permutation of $b$, because the values at free nodes are not modified by $C$.
Therefore, there are $f$ rows in $C$ corresponding to the free nodes $F$,
which when vertically stacked together form a permutation of $f\times f$
identity matrix. Therefore, the column span of $C$ has dimension $f$. Hence
$f\leq\mathrm{dim}(\mathcal{N}).$ Conversely, consider any
$\eta\in\mathcal{N}$. Given the entries $b$ of $\eta$ at the locations of free
parameters, then the rest of the entries of $\eta$ are determined by
$\eta=Cb.$ Therefore $\eta$ must lie in the column span of $C$. Therefore
$f=\mathrm{dim}(\mathcal{N}).$

##  Appendix E More details on optimization algorithms

#### Generic Quadratic Programming on the Dual

Recall that KTF solves the following convex optimization problem:

$\hat{\theta}=\mathop{\mathrm{argmin}}_{\theta\in\mathbb{R}^{n}}\;\frac{1}{2}\|y-\theta\|_{2}^{2}+\lambda\|D\theta\|_{1}.$
(101)

with $D=D_{n,d}^{(k+1)}$. The corresponding Lagrange dual problem is

$\displaystyle\max_{u}$
$\displaystyle\;\;-\frac{1}{2}uDD^{\mathsf{T}}u+y^{\mathsf{T}}D^{\mathsf{T}}u$
(102) subject to $\displaystyle\;\;-\lambda\leq u\leq\lambda.$

Note that the dual problem is a standard quadratic program (QP) and can be
solved using the interior point method (IPM) to high precision. Then the
primal solution can be constructed using $\hat{\theta}=y-D^{\mathsf{T}}u^{*}$
using the optimal solution $u^{*}$. In practice, IPM takes only a few
iterations to converge444In theory it could take up to $O(n^{1/2})$ iterations
to converge., but each iteration involves solving a linear system. This linear
system is sparse since $D$ is sparse — it has only $O(dkn)$ non-zero elements.
However, the condition number of the linear system grows as the weights on the
barrier function increase, which makes it difficult to exploit the sparsity
using methods such as preconditioned conjugates gradient method. On the other
hand, direct solvers such as Gaussian elimination and Cholesky decomposition
can take up to $O(n^{3})$. Sometimes this can be improved by exploiting the
banded structure of the linear system, we will describe a particular version
of the interior point method using logarithmic-barrier function.

#### Primal-Dual Interior Point method

The primal-dual version of the interior point solver proposed by (Kim et al.,
2009) for $\ell_{1}$ trend filtering can be straightforwardly applied to any
generalized lasso problem, including KTF. The main idea is to trace a “central
path” using Newton’s method with an increasing weights $t$ on the logarithmic
barrier functions. The computation is dominated by computing the search
direction of the Newton step, which boils down to solving the following system
of linear equations

$\begin{bmatrix}DD^{\mathsf{T}}&I&-I\\\ I&J_{1}&0\\\
-I&0&J_{2}\end{bmatrix}\begin{bmatrix}\Delta u\\\ \Delta\mu_{1}\\\
\Delta\mu_{2}\end{bmatrix}=-\begin{bmatrix}DD^{T}u-Dy+\mu_{1}-\mu_{2}\\\
f_{1}+(1/t)\mu_{1}^{-1}\\\ f_{2}+(1/t)\mu_{2}^{-1}\end{bmatrix}$ (103)

where $\mu_{1},\mu_{2}\in\mathbb{R}^{m}$ are the dual variables of the dual
problem (102), $f_{1}=u-\lambda\mathbf{1}$, $f_{2}=-u-\lambda\mathbf{1}$,
$J_{i}=\mathrm{diag}{\mu_{i}}^{-1}\mathrm{diag}(f_{i})$ are diagonal matrices
and $\mu_{i}^{-1}$ denotes entrywise inversion. Following the derivation of
(Kim et al., 2009), we can further eliminate $\Delta\mu_{1}$ and
$\Delta\mu_{2}$ and solve a linear system of the form

$(DD^{\mathsf{T}}-J_{1}^{-1}J_{2}^{-1})\Delta
u=-(DD^{\mathsf{T}}u-D^{(k+1)}y-(1/t)f_{1}^{-1}+(1/t)f_{2}^{-1}).$ (104)

and then construct the remainder of the solutions using

$\displaystyle\Delta\mu_{1}=-(\mu_{1}+(1/t)f_{1}^{-1}+J_{1}^{-1}\Delta u),$
$\displaystyle\Delta\mu_{2}=-(\mu_{2}+(1/t)f_{2}^{-1}-J_{2}^{-1}\Delta u).$

Unlike in the trend filtering problems in 1D where (104) is a banded linear
system with a bandwidth at most $2k+3$, in a $d$-dimensional grid, the linear
system is the following:

$\begin{bmatrix}\big{(}D_{\mathrm{1d}}D_{\mathrm{1d}}^{\mathsf{T}}\big{)}\otimes
I\otimes...\otimes I,&D_{\mathrm{1d}}\otimes
D_{\mathrm{1d}}^{\mathsf{T}}\otimes I\otimes...\otimes
I,&\cdots,&D_{\mathrm{1d}}\otimes I\otimes...\otimes I\otimes
D_{\mathrm{1d}}^{\mathsf{T}},\\\ \ D_{\mathrm{1d}}^{\mathsf{T}}\otimes
D_{\mathrm{1d}}\otimes I\otimes...\otimes
I,&I\otimes\big{(}D_{\mathrm{1d}}D_{\mathrm{1d}}^{\mathsf{T}}\big{)}\otimes
I\otimes...\otimes I,&\cdots,&I\otimes D_{\mathrm{1d}}\otimes
I\otimes...\otimes I\otimes D_{\mathrm{1d}}^{\mathsf{T}}\\\
\vdots&\vdots&\ddots&\vdots\\\ D_{\mathrm{1d}}^{\mathsf{T}}\otimes
I\otimes...\otimes I\otimes D_{\mathrm{1d}}&I\otimes
D_{\mathrm{1d}}^{\mathsf{T}}\otimes I\otimes...\otimes I\otimes
D_{\mathrm{1d}}&\cdots,&I\otimes...\otimes
I\otimes\big{(}D_{\mathrm{1d}}D_{\mathrm{1d}}^{\mathsf{T}}\big{)}\end{bmatrix}-J_{1}^{-1}J_{2}^{-1}$

where $D_{\mathrm{1d}}=D_{N,1}^{(k+1)}$. This is still sparse, structured, but
the bandedness is on the order of $O(n^{1-1/d}+k^{2})$. Moreover, the above
matrix is not full rank, and the condition number of the linear system blows
up as the dual variables $\mu_{1},\mu_{2}$ converges to $0$ with
$t\rightarrow\infty$.

#### Proximal Dykstra’s algorithm

Proximal Dykstra’s algorithm is an operator-splitting method for solving
problems of the form

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;\frac{1}{2}\|y-\theta\|_{2}^{2}+r_{1}(\theta)+r_{2}(\theta)+...+r_{d}(\theta)$
(105)

where $r_{1},...,r_{d}$ are convex but possibly non-smooth functions. We can
clearly see that the regularizer in KTF decomposes into this form with

$r_{i}(\theta)=\sum_{j_{1},j_{2},\dots,j_{i-1},j_{i+1},\dots,j_{d}}\|\tilde{D}^{(k+1)}_{\mathrm{1d}}\theta[j_{1},\dots,j_{i-1},:,j_{i+1},\dots,j_{d}]\|_{1}$

where we reshape $\theta$ into the $d$-dimensional tensor format and use the
“numpy.ndarray” notation.

The proximal Dykstra algorithm (see, e.g., Tibshirani, 2017) initializes
$\theta^{(0)}=y,z^{(-d+1)}=...=z^{(0)}=0$ and then iteratively applies the
following update rule for $t=1,2,3,...$:

$\displaystyle\theta^{(t)}$
$\displaystyle=\mathsf{prox}_{r_{t\;\mathsf{mod}\;d+1}}(\theta^{(t-1)}+z^{(t-d)})$
$\displaystyle z^{(t)}$ $\displaystyle=\theta^{(t-1)}+z^{(t-d)}-\theta^{(t)}.$

where $\cdot\;\mathsf{mod}\;\cdot$ is the modulo operator, and the proximal
operator

$\mathsf{prox}_{r}(u)=\mathop{\mathrm{argmin}}_{\theta}\frac{1}{2}\|u-\theta\|^{2}+r(\theta).$

Note that this is equivalent to a cyclic block coordinate descent in the dual.

For KTF, each proximal problem can be parallelized (Barbero and Sra, 2014).
Specifically, on a $d$-dim regular grid, the proximal operator of $r_{i}$
further splits into solving $O(n^{1-1/d})$ 1D-trend filters of size $n^{1/d}$
in parallel. Each subproblem can be solved efficiently in $O(n^{1.5/d})$ time
with the primal-dual interior point method for $k\geq 1$ (Tibshirani, 2014)
and in linear time when $k=0$ using dynamic programming (Johnson, 2013).

#### Douglas-Rachford Splitting

Another operator-splitting method for solving KTF is through the Douglas-
Rachford (DR) algorithm (Eckstein and Bertsekas, 1992). For simplicity, we
will focus our discussion on the case of 2D grids. The DR algorithm
generically solves the following unconstrained problem:

$\mathop{\mathrm{minimize}}_{\theta}f(\theta)+g(\theta)$ (106)

for convex functions $f,g$. The update rules include initializing an auxiliary
variable $z^{(0)}=y$ and applying the following for $t=0,1,2,...$:

$\displaystyle\theta^{(t+1)}$ $\displaystyle=\mathsf{prox}_{f}(z^{(t)})$
$\displaystyle z^{(t+1)}$
$\displaystyle=z^{(t)}+\mathsf{prox}_{g}(2\theta^{(t+1)}-z^{(t)})-\theta^{(t+1)}.$

There are multiple ways of applying this to our problem. We apply the DR
algorithm to the dual of the following reformulation according to (Barbero and
Sra, 2014, Algorithm 9):

$\mathop{\mathrm{minimize}}_{\theta\in\mathbb{R}^{n}}\;\frac{1}{2}\|\theta\|_{2}^{2}+\left(\lambda\|\tilde{D}^{(k+1)}_{\mathrm{1d}}\otimes
I\theta\|_{1}-\langle\theta,y\rangle\right)+\left(\lambda\|I\otimes\tilde{D}^{(k+1)}_{\mathrm{1d}}\theta\|_{1}\right).$
(107)

We refer interested readers to (Barbero and Sra, 2014) for the derivation of
the dual and the conversion of the problem into one that resembles (106).
Ultimately, the proximal operator of the conjugate function (an indicator on a
certain polytope) can be evaluated using the proximal operator of the $r_{1}$
and $r_{2}$ as in the proximal Dykstra updates via the Moreau decomposition:

$\mathsf{prox}_{r_{i}}(u)+\mathsf{prox}_{r_{i}^{*}}(u)=u.$

In other words, the Douglas-Rachford algorithm enjoys the same computational
benefits of the proximal Dykstra’s algorithm as each proximal operator
evaluation involves only solving 1D trend filtering problems in parallel.

## References

  * Almansa et al. (2008) Andrés Almansa, Coloma Ballester, Vicent Caselles, and Gloria Haro.  A TV based restoration model with local constraints.  _Journal of Scientific Computing_ , 34(3):209–236, 2008.
  * Barbero and Sra (2014) Alvaro Barbero and Suvrit Sra.  Modular proximal optimization for multidimensional total-variation regularization.  arXiv: 1411.0589, 2014.
  * Bassett and Sharpnack (2019) Robert Bassett and James Sharpnack.  Fused density estimation: Theory and methods.  _Journal of Royal Statistical Society, Series B_ , 81, 2019.
  * Bibaut and van der Laan (2019) Aurélien F Bibaut and Mark J van der Laan.  Fast rates for empirical risk minimization over càdlàg functions with bounded sectional variation norm.  arXiv: 1907.09244, 2019.
  * Birge and Massart (2001) Lucien Birge and Pascal Massart.  Gaussian model selection.  _Journal of the European Mathematical Society_ , 3(3):203–268, 2001.
  * Boyd et al. (2011) Steve Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein.  Distributed optimization and statistical learning via the alternative direction method of multipliers.  _Foundations and Trends in Machine Learning_ , 3(1):1–122, 2011.
  * Candès and Guo (2002) Emmanuel J. Candès and Franck Guo.  New multiscale transforms, minimum total variation synthesis: Applications to edge-preserving image reconstruction.  _Signal Processing_ , 82(11):1519–1543, 2002\.
  * Chambolle (2004) Antonin Chambolle.  An algorithm for total variation minimization and applications.  _Journal of Mathematical Imaging and Vision_ , 20(1):89–97, 2004.
  * Chambolle (2005) Antonin Chambolle.  Total variation minimization and a class of binary MRF models.  In David Hutchison, Takeo Kanade, Josef Kittler, Jon M. Kleinberg, Friedemann Mattern, John C. Mitchell, Moni Naor, Oscar Nierstrasz, C. Pandu Rangan, Bernhard Steffen, Madhu Sudan, Demetri Terzopoulos, Dough Tygar, Moshe Y. Vardi, Gerhard Weikum, Anand Rangarajan, Baba Vemuri, and Alan L. Yuille, editors, _Energy Minimization Methods in Computer Vision and Pattern Recognition_ , pages 136–152. Springer, 2005.
  * Chambolle and Lions (1997) Antonin Chambolle and Pierre-Louis Lions.  Image recovery via total variation minimization and related problems.  _Numerische Mathematik_ , 76(2):167–188, 1997\.
  * Chan et al. (2000) Tony Chan, Antonio Marquina, and Pep Mulet.  High-order total variation-based image restoration.  _SIAM Journal on Scientific Computing_ , 22(2):503–516, 2000.
  * Chan and Esedoglu (2005) Tony F. Chan and Selim Esedoglu.  Aspects of total variation regularized $L^{1}$ function approximation.  _SIAM Journal on Applied Mathematics_ , 65(5):1817–1837, 2005.
  * Chatterjee and Goswami (2021) Sabyasachi Chatterjee and Subhajit Goswami.  Adaptive estimation of multivariate piecewise polynomials and bounded variation functions by optimal decision trees.  _Annals of Statistics_ , 49(5):2531–2551, 2021\.
  * Chui et al. (1992a) Charles K. Chui, Jeffrey M. Lemm, and Sahra Sedigh.  _An Introduction to Wavelets_.  Academic Press, 1992a.
  * Chui et al. (1992b) Charles K. Chui, Joachim Stöckler, and Joseph D. Ward.  Compactly supported box-spline wavelets.  _Approximation Theory and its Applications_ , 8(3):77–100, 1992b.
  * Daubechies (1992) Ingrid Daubechies.  _Ten Lectures on Wavelets_.  Society for Industrial and Applied Mathematics, 1992.
  * del Álamo et al. (2021) Miguel del Álamo, Housen Li, and Axel Munk.  Frame-constrained total variation regularization for white noise regression.  _Annals of Statistics_ , 49(3), 2021.
  * DeVore and Lorentz (1993) Ronald DeVore and George Lorentz.  _Constructive Approximation_.  Springer, 1993.
  * DeVore and Lucier (1992) Ronald A. DeVore and Bradley J. Lucier.  Wavelets.  _Acta Numerica_ , 1:1–56, 1992.
  * DeVore et al. (1998) Ronald A. DeVore, Sergei V. Konyagin, and Vladimir N. Temlyakov.  Hyperbolic wavelet approximation.  _Constructive Approximation_ , 14(1):1–26, 1998\.
  * Dong et al. (2011) Yiqiu Dong, Michael Hintermüller, and M. Monserrat Rincon-Camacho.  Automated regularization parameter selection in multi-scale total variation models for image restoration.  _Journal of Mathematical Imaging and Vision_ , 40(1):82–104, 2011.
  * Donoho and Johnstone (1998) David L. Donoho and Iain M. Johnstone.  Minimax estimation via wavelet shrinkage.  _Annals of Statistics_ , 26(8):879–921, 1998\.
  * Eckstein and Bertsekas (1992) Jonathan Eckstein and Dimitri P. Bertsekas.  On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators.  _Mathematical Programming_ , 55(1):293–318, 1992\.
  * Efron (1986) Bradley Efron.  How biased is the apparent error rate of a prediction rule?  _Journal of the American Statistical Association_ , 81(394):461–470, 1986.
  * Evans and Gariepy (2015) Lawrence C. Evans and Ronald F. Gariepy.  _Measure Theory and Fine Properties of Functions_.  CRC Press, 2015.  Revised edition.
  * Fang et al. (2021) Billy Fang, Adityanand Guntuboyina, and Bodhisattva Sen.  Multivariate extensions of isotonic regression and total variation denoising via entire monotonicity and Hardy-Krause variation.  _Annals of Statistics_ , 49(2):769–792, 2021\.
  * Gavish et al. (2010) Matan Gavish, Boaz Nadler, and Ronald Coifman.  Multiscale wavelets on trees, graphs and high dimensional data: theory and applications to semi supervised learning.  In _Proceedings of the Annual Conference on Learning Theory_ , 2010\.
  * Göbel et al. (2018) Franziska Göbel, Gilles Blanchard, and Ulrike von Luxburg.  Construction of tight frames on graphs and application to denoising.  In _Handbook of Big Data Analytics_ , pages 503–522. Springer, 2018\.
  * Goldenshluger and Lepski (2008) Alexander Goldenshluger and Oleg Lepski.  Universal pointwise selection rule in multivariate function estimation.  _Bernoulli_ , 14(4):1150–1190, 2008.
  * Goldenshluger and Lepski (2009) Alexander Goldenshluger and Oleg Lepski.  Structural adaptation via $L_{p}$-norm oracle inequalities.  _Probability Theory and Related Fields_ , 143(1–2):41–71, 2009.
  * Goldenshluger and Lepski (2011) Alexander Goldenshluger and Oleg Lepski.  Bandwidth selection in kernel density estimation: Oracle inequalities and adaptive minimax optimality.  _Annals of Statistics_ , 39(3):1608–1632, 2011\.
  * Goldenshluger and Lepski (2013) Alexander Goldenshluger and Oleg Lepski.  General selection rule from a family of linear estimators.  _Theory of Probability & Its Applications_, 57(2):209–226, 2013.
  * Goldenshluger and Lepski (2014) Alexander Goldenshluger and Oleg Lepski.  On adaptive minimax density estimation on $R^{d}$.  _Probability Theory and Related Fields_ , 159(3):479–543, 2014.
  * Haris et al. (2019) Asad Haris, Noah Simon, and Ali Shojaie.  Generalized sparse additive models.  _arXiv:1903.04641 [math, stat]_ , March 2019.
  * Hastie and Tibshirani (1990) Trevor Hastie and Robert Tibshirani.  _Generalized Additive Models_.  Chapman & Hall, 1990.
  * Hayfield and Racine (2008) Tristen Hayfield and Jeffrey S. Racine.  Nonparametric econometrics: The np package.  _Journal of Statistical Software_ , 27(5), 2008.  URL https://www.jstatsoft.org/v27/i05/.
  * Hutter and Rigollet (2016) Jan-Christian Hutter and Philippe Rigollet.  Optimal rates for total variation denoising.  In _Proceedings of the Annual Conference on Learning Theory_ , 2016\.
  * Johnson (2013) Nicholas Johnson.  A dynamic programming algorithm for the fused lasso and $l_{0}$-segmentation.  _Journal of Computational and Graphical Statistics_ , 22(2):246–260, 2013.
  * Johnstone (2015) Iain M. Johnstone.  _Gaussian Estimation: Sequence and Wavelet Models_.  Cambridge University Press, 2015.  Draft version.
  * Kerkyacharian et al. (2001) Gérard Kerkyacharian, Oleg V. Lepski, and Dominique Picard.  Nonlinear estimation in anisotropic multi-index denoising.  _Probability Theory and Related Fields_ , 121(2):137–170, 2001.
  * Kerkyacharian et al. (2008) Gérard Kerkyacharian, Oleg V. Lepski, and Dominique Picard.  Nonlinear estimation in anisotropic multi-index denoising. Sparse case.  _Theory of Probability & Its Applications_, 52(1):58–77, 2008.
  * Ki et al. (2021) Dohyeong Ki, Billy Fang, and Adityanand Guntuboyina.  MARS via LASSO.  arXiv: 2111.11694, 2021.
  * Kim et al. (2009) Seung-Jean Kim, Kwangmoo Koh, Stephen Boyd, and Dimitry Gorinevsky.  $\ell_{1}$ trend filtering.  _SIAM Review_ , 51(2):339–360, 2009.
  * Koenker et al. (1994) Roger Koenker, Pin Ng, and Stephen Portnoy.  Quantile smoothing splines.  _Biometrika_ , 81(4):673–680, 1994.
  * Korostelev and Tsybakov (2003) Aleksandr P. Korostelev and Alexandre B. Tsybakov.  _Minimax Theory of Image Reconstructions_.  Springer, 2003.
  * Lepski (2015) Oleg V. Lepski.  Adaptive estimation over anisotropic functional classes via oracle approach.  _Annals of Statistics_ , 43(3):1178–1242, 2015\.
  * Lepski and Spokoiny (1997) Oleg V. Lepski and Vladimir G. Spokoiny.  Optimal pointwise adaptive methods in nonparametric estimation.  _Annals of Statistics_ , 25(6):2512–2546, 1997\.
  * Lepski et al. (1997) Oleg V. Lepski, Enno Mammen, and Vladimir G. Spokoiny.  Optimal spatial adaptation to inhomogeneous smoothness: An approach based on kernel estimates with variable bandwidth selectors.  _Annals of Statistics_ , 25(3):929–947, 1997\.
  * Lepskii (1991) Oleg V. Lepskii.  On a problem of adaptive estimation in Gaussian white noise.  _Theory of Probability & Its Applications_, 35(3):454–466, 1991.
  * Lepskii (1992) Oleg V. Lepskii.  Asymptotically minimax adaptive estimation. I: Upper bounds. Optimally adaptive estimates.  _Theory of Probability & Its Applications_, 36(4):682–697, 1992.
  * Lepskii (1993) Oleg V. Lepskii.  Asymptotically minimax adaptive estimation. II. Schemes without optimal adaptation: Adaptive estimators.  _Theory of Probability & Its Applications_, 37(3):433–448, 1993.
  * Lin et al. (2017) Kevin Lin, James Sharpnack, Alessandro Rinaldo, and Ryan J. Tibshirani.  A sharp error analysis for the fused lasso, with application to approximate changepoint screening.  In _Advances in Neural Information Processing Systems_ , 2017.
  * Lorentz and Madych (1992) Rudolph A. H. Lorentz and Wolodymyr R. Madych.  Wavelets and generalized box splines.  _Applicable Analysis_ , 44(1–2):51–76, 1992\.
  * Mallat (2009) Stephane Mallat.  _A Wavelet Tour of Signal Processing_.  Academic Press, 2009.  Third edition.
  * Mallat (1989a) Stephane G. Mallat.  Multiresolution approximations and wavelet orthonormal bases of $L^{2}(R)$.  _Transactions of the American Mathematical Society_ , 315(1):69–87, 1989a.
  * Mallat (1989b) Stephane G. Mallat.  A theory for multiresolution signal decomposition: The wavelet representation.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 11(7):674–693, 1989b.
  * Mammen (1991) Enno Mammen.  Nonparametric regression under qualitative smoothness assumptions.  _Annals of Statistics_ , 19(2):741–759, 1991\.
  * Mammen and van de Geer (1997) Enno Mammen and Sara van de Geer.  Locally apadtive regression splines.  _Annals of Statistics_ , 25(1):387–413, 1997\.
  * Meyer (1987) Yves Meyer.  Principe d’incertitude, bases Hilbertiennes et algebres d’operateurs.  _Séminaire Bourbaki_ , 145–146(662):209–223, 1987.
  * Meyer (1990) Yves Meyer.  _Ondelettes et Opérateurs_.  Hermann, 1990.
  * Meyer and Roques (1993) Yves Meyer and Sylvie Roques.  _Progress in Wavelet Analysis and Applications_.  Atlantica Séguier Frontières, 1993.
  * Nason (2016) Guy Nason.  _wavethresh: Wavelets Statistics and Transforms_ , 2016.  URL https://CRAN.R-project.org/package=wavethresh.  R package version 4.6.8.
  * Neuman and Schonbach (1974) Charles P. Neuman and Dave I. Schonbach.  Discrete (Legendre) orthogonal polynomials—a survey.  _International Journal for Numerical Methods in Engineering_ , 8(4):743–770, 1974.
  * Neumann (2000) Michael H. Neumann.  Multivariate wavelet thresholding in anisotropic function spaces.  _Statistica Sinica_ , 10(2):399–431, 2000.
  * Neumann and von Sachs (1997) Michael H. Neumann and Rainer von Sachs.  Wavelet thresholding in anisotropic function classes and application to adaptive estimation of evolutionary spectra.  _Annals of Statistics_ , 25(1):38–76, 1997.
  * Ortelli and van de Geer (2021) Francesco Ortelli and Sara van de Geer.  Tensor denoising with trend filtering.  arXiv: 2101.10692, 2021.
  * Padilla and Scott (2016) Oscar Hernan Madrid Padilla and James G. Scott.  Nonparametric density estimation by histogram trend filtering.  arXiv: 1509.04348, 2016.
  * Padilla et al. (2018) Oscar Hernan Madrid Padilla, James Sharpnack, James G. Scott, and Ryan J. Tibshirani.  The DFS fused lasso: Linear-time denoising over general graphs.  _Journal of Machine Learning Research_ , 18:176–1, 2018\.
  * Padilla et al. (2020) Oscar Hernan Madrid Padilla, James Sharpnack, Yanzhen Chen, and Daniela Witten.  Adaptive non-parametric regression with the k-nn fused lasso.  _Biometrika_ , 107(2):293–310, 2020.
  * Ramdas and Tibshirani (2016) Aaditya Ramdas and Ryan J. Tibshirani.  Fast and flexible ADMM algorithms for trend filtering.  _Journal of Computational and Graphical Statistics_ , 25(3):839–858, 2016.
  * Riemenschneider and Shen (1992) Sherman D. Riemenschneider and Zuowei Shen.  Wavelets and pre-wavelets in low dimensions.  _Journal of Approximation Theory_ , 71(1):18–38, 1992.
  * Rudin and Osher (1994) Leonid I. Rudin and Stanley Osher.  Total variation based image restoration with free local constraints.  In _Proceedings of 1st International Conference on Image Processing_ , volume 1, pages 31–35, Austin, TX, USA, 1994. IEEE Comput. Soc. Press.
  * Rudin et al. (1992) Leonid I. Rudin, Stanley Osher, and Emad Faterni.  Nonlinear total variation based noise removal algorithms.  _Physica D: Nonlinear Phenomena_ , 60(1):259–268, 1992.
  * Sadhanala (2019) Veeranjaneyulu Sadhanala.  _Nonparametric Methods with Total Variation Type Regularization_.  PhD thesis, Machine Learning Department, Carnegie Mellon University, 2019\.
  * Sadhanala and Tibshirani (2019) Veeranjaneyulu Sadhanala and Ryan J. Tibshirani.  Additive models via trend filtering.  _Annals of Statistics_ , 47(6):3032–3068, 2019\.
  * Sadhanala et al. (2016) Veeranjaneyulu Sadhanala, Yu-Xiang Wang, and Ryan J. Tibshirani.  Total variation classes beyond 1d: Minimax rates, and the limitations of linear smoothers.  In _Advances in Neural Information Processing Systems_ , volume 29, 2016.
  * Sadhanala et al. (2017) Veeranjaneyulu Sadhanala, Yu-Xiang Wang, James Sharpnack, and Ryan J. Tibshirani.  Higher-total variation classes on grids: Minimax theory and trend filtering methods.  In _Advances in Neural Information Processing Systems_ , volume 30, 2017.
  * Sharpnack et al. (2013) James Sharpnack, Aarti Singh, and Akshay Krishnamurthy.  Detecting activations over graphs using spanning tree wavelet bases.  In _Proceedings of the International Conference on Artificial Intelligence and Statistics_ , 2013.
  * Steidl et al. (2006) Gabriel Steidl, Stephan Didas, and Julia Neumann.  Splines in higher order TV regularization.  _International Journal of Computer Vision_ , 70(3):214–255, 2006.
  * Stein (1981) Charles Stein.  Estimation of the mean of a multivariate normal distribution.  _Annals of Statistics_ , 9(6):1135–1151, 1981\.
  * Tansey and Scott (2015) Wesley Tansey and James Scott.  A fast and flexible algorithm for the graph-fused lasso.  arXiv: 1505.06475, 2015.
  * Tibshirani (2014) Ryan J. Tibshirani.  Adaptive piecewise polynomial estimation via trend filtering.  _Annals of Statistics_ , 42(1):285–323, 2014\.
  * Tibshirani (2017) Ryan J. Tibshirani.  Dykstra’s algorithm, ADMM, and coordinate descent: Connections, insights, and extensions.  In _Advances in Neural Information Processing Systems_ , volume 30, 2017.
  * Tibshirani (2020) Ryan J. Tibshirani.  Divided differences, falling factorials, and discrete splines: Another look at trend filtering and related problems.  arXiv: 2003.03886, 2020.
  * Tibshirani and Taylor (2011) Ryan J. Tibshirani and Jonathan Taylor.  The solution path of the generalized lasso.  _Annals of Statistics_ , 39(3):1335–1371, 2011\.
  * Tibshirani and Taylor (2012) Ryan J. Tibshirani and Jonathan Taylor.  Degrees of freedom in lasso problems.  _Annals of Statistics_ , 40(2):1198–1232, 2012\.
  * Tsybakov (2009) Alexandre B. Tsybakov.  _Introduction to Nonparametric Estimation_.  Springer, 2009.
  * Vogel and Oman (1996) Curtis R. Vogel and M. E. Oman.  Iterative methods for total variation denoising.  _SIAM Journal on Scientific Computing_ , 17(1):227–238, 1996.
  * Wang et al. (2014) Yu-Xiang Wang, Alexander Smola, and Ryan J. Tibshirani.  The falling factorial basis and its statistical applications.  In _Proceedings of the International Conference on Machine Learning_ , 2014.
  * Wang et al. (2016) Yu-Xiang Wang, James Sharpnack, Alex Smola, and Ryan J. Tibshirani.  Trend filtering on graphs.  _Journal of Machine Learning Research_ , 17(105):1–41, 2016.
  * Ye and Padilla (2021) Steven Siwei Ye and Oscar Hernan Madrid Padilla.  Non-parametric quantile regression via the k-nn fused lasso.  _Journal of Machine Learning Research_ , 22(111):1–38, 2021.