Multivariate Trend Filtering for Lattice Data

TL;DR

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 is studied, revealing a number of interesting phenomena, including the dominance of KTF over linear smoothers in estimating heterogeneously smooth functions.

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^{\mathrm{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$).

Figures (6)

• 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$.
• 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$): KTF achieves the minimax rate (up to log factors) over $\mathcal{T}_{n,d}^k(C_n^*)$; andno 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.
• 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 curve is hidden behind the red curve).
• Figure 5: Illustration of multivariate interpolation from Algorithm \ref{['alg:interp']}, 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 \ref{['alg:interp_1d']}. 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 \ref{['alg:interp_1d']}.
• Figure 7: MSE of various methods for estimating the signal function $f_0$ in Figure \ref{['fig:intro']}, over a square lattice with $n=256^2=65536$ points, and subject to Gaussian noise that yields 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. As their tuning parameters lie on different scales, their values are scaled to fit onto a single x-axis. The takeaway is that KTF and GTF, in particular for $k=2$, achieve clearly the best performance.

Theorems & Definitions (49)

• Proposition 1
• Proposition 2
• Remark 1
• Remark 2
• Remark 3
• Theorem 1
• Remark 4
• Remark 5
• Proposition 3: sadhanala2017higher
• Theorem 2: wang2016trend