---
abstract: 'Splines are a popular and attractive way of smoothing noisy data. Computing splines involves minimizing a functional which is a linear combination of a fitting term and a regularization term. The former is classically computed using a (weighted) L2 norm while the latter ensures smoothness. Thus, when dealing with grid data, the optimization can be solved very efficiently using the DCT. In this work we propose to replace the L2 norm in the fitting term with an L1 norm, leading to automatic robustness to outliers. To solve the resulting minimization problem we propose an extremely simple and efficient numerical scheme based on split-Bregman iteration combined with DCT. Experimental validation shows the high-quality results obtained in short processing times.'
address: ' Department of Electrical and Computer Engineering, Duke University.[^1]\'
author:
- 
- 
bibliography:
- 'mtepper.bib'
title: 'L1 Splines for Robust, Simple, and Fast Smoothing of Grid Data'
---

Introduction
============

Smoothing a dataset consists in finding an approximating function that captures important patterns in the data, while disregarding noise or other fine-scale structures. Let $y \in {{\mathbb R}}^{n_1 \times \dots \times n_m} \rightarrow {{\mathbb C}}$ be an $m$-dimensional discrete signal, where $n_j$ ($1 \leq j \leq d$) is the domain of $y$ along the $j$-th dimension. We can model $y$ by $$y = \hat{y} + r ,
    \label{eq:model}$$ where $r$ represents some noise and $\hat{y}$ is a smooth function. A very common regularization choice is to enforce $C^2$ continuity, in which case $\hat{y}$ is called a (cubic) spline. Smoothing $y$ relies upon finding the best estimate of $\hat{y}$ under the proper smoothness and noise assumptions. We can approximate $\hat{y}$ by minimizing an objective functional $${\mathrm{F}}(z) = {\mathrm{R}}_y (z) + s {\mathrm{P}}(z) ,$$ where ${\mathrm{R}}_y (z)$ is a data fitting term, defined by the distribution of $r$, ${\mathrm{P}}(z)$ is a regularization term, and $s$ is a scalar that determines the balance between both terms. Such scalar parameter can be automatically derived using Bayesian or MDL techniques, as will be later shown.

For clarity, we will describe in depth the case $m=1$, where we have $n=n_1$ samples. We extend the results later to the general $m$-dimensional case.

Let us begin by explaining the smoothing term. The $C^2$ continuity requirement leads to define ${\mathrm{P}}(z) = {\left\| Dz \right\|_{2}}^ 2$, $D$ being a discrete *second-order* differential operator, defined $\forall i, 2 \leq i \leq n-1$, by $$\begin{aligned}
    D_{i, i-1} &= \tfrac{2}{h_{i-1} (h_{i-1} + h_{i})} ,\\
    D_{i, i} &= \tfrac{-2}{h_{i-1} \ h_{i}} ,\\
    D_{i-1, i} &= \tfrac{2}{h_{i} (h_{i-1} + h_{i})} ,\end{aligned}$$ where $h_i$ represents the step, or sampling rate, between $y_i$ and $y_{i+1}$. Assuming repeating border elements, that is, $y_0 = y_1$ and $y_{n+1}  =y_n$, gives $D_{1, 1} = -D_{1, 2} = -1/h_1^2$, and $D_{n, n-1} = -D_{n, n} = -1/h_{n-1}^2$.

Regarding the fitting term, the classical assumption is that the noise $r$ in Equation () has Gaussian distribution with zero mean and unknown variance, which leads to setting ${\mathrm{R}}_y (z) = {\left\| z - y \right\|_{2}}^2$. Smoothing then can be formulated as the least-squares regression $$\hat{y} = {\underset{z}{\operatorname{argmin}}\ } {\left\| z - y \right\|_{2}}^2 + s {\left\| Dz \right\|_{2}}^ 2 .$$ For clarity, we call the estimate $\hat{y}$ obtained with this method an L2 spline.

It is a well known fact that least squares estimates for regression models are highly non-robust to outliers. Although there is no agreement on a universal and formal definition of an outlier, it is usually regarded as an observation that does not follow the patterns in the data. Notice that smoothing should produce an estimate $\hat{y}$ taking into account only important patterns, that is, the inliers, in the data. In this sense, the L2 formulation cannot correctly handle outliers by itself.

In order to solve this problem, we propose to take a different assumption on the distribution of the noise $r$ in Equation (). By choosing a distribution with fatter tails than the Gaussian distribution, the derived estimator will correctly handle outliers. We thus assume that $r$ follows a Laplace distribution with zero mean and unknown scale parameter, a common practice in other problems as we will further discuss below. This leads to the fitting term ${\mathrm{R}}_ y(z) = {\left\| z - y \right\|_{1}}$ and the regression then becomes $$\hat{y} = {\underset{z}{\operatorname{argmin}}\ } \ {\left\| z - y \right\|_{1}} + s {\left\| Dz \right\|_{2}}^2 .$$ We call the estimate $\hat{y}$ obtained with this formulation an L1 spline.

Let us point out that the use of L1 fitting terms for solving inverse problems is not new. For example, in 2001, Nikolova proved the theoretic pertinence of using L1 fitting terms for image denoising [@nikolova02]. Other interesting works have addressed this approach for total variation image denoising [@alliney97; @chan2004; @aujol06; @nikolova12] or total variation optical flow [@zach07; @wedel09; @raket11] (this robustness-to-outliers type of ideas was previously introduced in the context of optical flow by Black and Anandan [@black91]). Recall that the total-variation regularization term involves first-order derivatives, while the proposed L1 splines, on the other hand, involve second-order derivatives.

Regularly sampled signals are extremely common in practice, and their analysis becomes easier and faster. In particular, we follow the most common choice when dealing with discrete $m$-dimensional data, which is assuming a “rectangular” Cartesian sampling pattern. When the sampling is isotropic, i.e., “square,” we refer to this type of data as grid data.

We developed an iterative algorithm for computing L1 splines, based on split-Bregman iteration [@goldstein09], that is specially suited for the case of grid data. This algorithm is extremely fast, both in running time and in the number of iterations until convergence. It is also outstandingly simple, making the implementation completely straightforward.

The remainder of the paper is structured as follows. In Section \[sec:L2splines\], we overview a fast algorithm to compute L2 splines and robust L2 splines, a modification of the least squares regression that allows to handle outliers. In Section \[sec:L1splines\], we present the proposed algorithm for computing L1 splines. Then, in Section \[sec:results\], we show results obtained with L1 splines, systematically outperforming its L2 and robust L2 counterparts in the presence of outliers. We also show that the proposed computational algorithm is very efficient. Finally, in Section \[sec:conclusions\] we provide some concluding remarks.

Smoothing splines {#sec:L2splines}
=================

As aforementioned, the classical assumption is that the noise $r$ in Equation () follows a Gaussian distribution with zero mean and unknown variance. This leads to solve the least-squares regression $$\hat{y} = {\underset{z}{\operatorname{argmin}}\ } {\left\| z - y \right\|_{2}}^2 + s {\left\| Dz \right\|_{2}}^ 2 .
    \label{eq:LSformulation}$$ Since both terms are differentiable, we obtain $$\hat{y} = (I + s {D^\mathrm{T}} D)^{-1} y .
     \label{eq:LSsolution}$$

Garcia proposed a very efficient method for dealing with regularly sampled data [@garcia10]. Assuming that the data are equally spaced, that is, without loss of generality $\forall i, h_i = 1$, we obtain $$D = 
\begin{pmatrix}
-1 & 1 \\
1 & -2 & 1 \\
& \ddots & \ddots & \ddots \\
&& 1 & -2 & 1 \\
&&& 1 & -1
\end{pmatrix} .$$ An eigendecomposition of $D$ yields $D = U \varLambda {U^\mathrm{T}}$, where $\varLambda$ is a diagonal matrix containing the eigenvalues of $D$, given by [@yueh05] $$\varLambda_{i,j} =
    \begin{cases}
        -2 + 2 \cos((i-1) \pi / n) , & \text{if $i=j$;} \\
        0 , & \text{otherwise.}
    \end{cases}$$ Since $U$ is a unitary matrix, we can write Equation () as $$\hat{y} = U (I + s \varLambda^2)^{-1} {U^\mathrm{T}} y .
     \label{eq:LSsolutionDecomposed}$$ Let us define the matrix $\varGamma = (I + s \varLambda^2)^{-1}$. Trivially, $$\varGamma_{i, j} =
    \begin{cases}
     \left[ 1 + s (-2 + 2 \cos((i-1) \pi / n))^2 \right]^{-1} , & \text{if $i=j$;} \\
     0 , & \text{otherwise.}
     \end{cases}
     \label{eq:varGamma}$$ Following Strang [@strang99] and Garcia [@garcia10], let us observe that ${U^\mathrm{T}}$ is a DCT-II matrix and $U$ is an inverse DCT-II matrix. Then, Equation () can be expressed as $$\hat{y} = {\mathrm{DCT^{-1}}} (\varGamma \ {\mathrm{DCT}} (y)) ,
    \label{eq:LSsolutionDCT}$$ where ${\mathrm{DCT}}(\cdot)$ and ${\mathrm{DCT^{-1}}}(\cdot)$ stand for the DCT-II and inverse DCT-II functions. Equation () provides a fast and simple algorithm for computing L2 splines.

Robust estimation. {#sec:robustL2}
------------------

Often in practice there are in $y$ some values $y_i$ that could not be observed (or recorded) for some reason. We would like to be able to handle such cases in such a way that the missing values are inferred from the ones that can be observed. Let $W$ be an $n \times n$ diagonal matrix such that $W_{i, i}$ represents a weight assigned to observation $i$. $W$ is defined by $$W_{i, i} = 
    \begin{cases}
        0 &\text{if datapoint $i$ is missing;} \\
        \rho &\text{otherwise.}
    \end{cases}$$ where $\rho$ is some arbitrary constant in $(0, 1]$; in practice, and without loss of generality, we set $\rho=1$. We can then solve $$\hat{y} = {\underset{z}{\operatorname{argmin}}\ } {\left\|  W^{1/2} (z - y) \right\|_{2}}^2 + s {\left\| Dz \right\|_{2}}^ 2 ,
    \label{eq:WLSformulation}$$ which will simply omit the missing points from the computation of the residual while the regularizer will still have a smoothing effect over both present and missing points. Equation () acts as an impainting algorithm, filling the missing values in such a way that continuity between filled values and smoothed ones is preserved. The minimization of Equation () gives $$(I + s {D^\mathrm{T}}D) \hat{y} = W (y - \hat{y}) + \hat{y} .$$ This leads to the iterative procedure $$\hat{y}^{k+1} = (I + s {D^\mathrm{T}}D)^{-1} \left( W \left( y - \hat{y}^k  \right) + \hat{y}^k \right) ,$$ which, similarly to Equation (), becomes $$\hat{y}^{k+1} = {\mathrm{DCT}}^{-1} \left( \varGamma \ {\mathrm{DCT}} \left( W \left(y - \hat{y}^k \right) + \hat{y}^k \right) \right) .
    \label{eq:WLSsolutionDCT}$$

On a different note, real data often present observations that lie abnormally far from their “true” value, i.e., that do not appear to follow the pattern of the other data points. The main drawback of the penalized least squares formulation Equation () is its sensitivity to these outliers. To address this issue, weights can be assigned to every point, as in Equation (), such that outliers exert less influence during the estimation process. In this case, the weights are iteratively refined during the estimation process using robust estimators for the mean and variance of the data. Defining these estimators is a complex problem by itself. For details about how $W$ can be set and updated for added robustness to outliers, refer to Garcia’s work [@garcia10].

L1 splines {#sec:L1splines}
==========

In this section we introduce a different splines formulation in order to handle outliers in the data. We assume that the noise $r$ in Equation () follows a Laplace distribution with zero mean and unknown scale parameter, which leads to ${\mathrm{R}}_ y(z) = {\left\| z - y \right\|_{1}}$. The regression then becomes $$\min_{z} \ {\left\| z - y \right\|_{1}} + s {\left\| Dz \right\|_{2}}^2 .
    \label{eq:L1formulation}$$

Goldstein and Osher [@goldstein09] proposed a very elegant and efficient algorithm for solving the L1 constrained problem (related to a number of very efficient optimization algorithms, e.g., see [@combettes11]) $$\min_u {\left\| \Phi(u) \right\|_{1}} + {\mathrm{H}}(u) .$$ For this, they consider the equivalent problem $$\min_u {\left\| d \right\|_{1}} + {\mathrm{H}}(u) \quad \text{s.t.} \quad d = \Phi(u) ,$$ which they first convert it into the unconstrained problem $$\min_u {\left\| d \right\|_{1}} + {\mathrm{H}}(u) + \tfrac{\lambda}{2} {\left\| d - \Phi(u) \right\|_{2}}^2 .$$ In this form, the penalty function does not accurately enforce the constraint for small $\lambda$. The constraint is enforced by letting $\lambda \rightarrow \infty$. However, another solution for this new formulation is found by using the following two-phase algorithm $$\begin{aligned}
    \left( u^{k+1}, d^{k+1} \right) &= {\underset{u}{\operatorname{argmin}}\ } {\left\| d \right\|_{1}} + {\mathrm{H}}(u) +
    \tfrac{\lambda}{2} {\left\| d - \Phi(u) - b^k \right\|_{2}}^2 ,\\
    b^{k+1} &= b^k + \left( \Phi(u^{k+1}) - d^k \right) .\end{aligned}$$ This algorithm is often denoted in the literature as split-Bregman iteration. This class of algorithms has several nice theoretical properties and has successfully been applied to several problems in practice such as image restoration [@osher05], image denoising [@xu07], compressed sensing [@yin08], and image segmentation [@goldstein10]; see also [@combettes11] and references therein.

We use this technique for solving Equation (). We begin by setting $\Phi(z) = z-y$ and $H(z) = s {\left\| Dz \right\|_{2}}^2$, which leads to the problem $$\min_{z, d} \ {\left\| d \right\|_{1}} + s {\left\| Dz \right\|_{2}}^ 2  \quad \text{s.t.} \quad d = z - y .$$ We then transform it into the unconstrained form $$\min_{z, d} \ {\left\| d \right\|_{1}} + s {\left\| Dz \right\|_{2}}^ 2 + \tfrac{\lambda}{2} {\left\| d - z + y \right\|_{2}}^2 ,
    \label{eq:L1formulationExtended}$$ and the Bregman iteration simply takes the form $$\begin{aligned}
    \left( z^{k+1}, d^{k+1} \right) & = {\underset{z, d}{\operatorname{argmin}}\ } \ {\left\| d \right\|_{1}} + s {\left\| Dz \right\|_{2}}^ 2 +
    \tfrac{\lambda}{2} {\left\| d - z + y - b^k \right\|_{2}}^2 ,
    \label{eq:splitBregman1}\\
    b^{k+1} & = b^k + (z^{k+1} - y - d^{k+1}) . \label{eq:splitBregman2}\end{aligned}$$

Because of the splitting of the L1 and L2 components in the functional (), we can perform this minimization efficiently by iteratively minimizing with respect to $z$ and $d$ separately, $$\begin{aligned}
    z^{k+1} &= {\underset{z}{\operatorname{argmin}}\ } \ s {\left\| Dz \right\|_{2}}^ 2 + \tfrac{\lambda}{2} {\left\| d^k - z + y - b^k \right\|_{2}}^2 \label{eq:splitBregman1_1} , \\
    d^{k+1} &= {\underset{d}{\operatorname{argmin}}\ } \ {\left\| d \right\|_{1}} + \tfrac{\lambda}{2} {\left\| d - z^{k+1} + y - b^k \right\|_{2}}^2 \label{eq:splitBregman1_2} .\end{aligned}$$ For minimizing Equation () we set $\tilde{y} = d^k + y - b^k$ and $\tilde{s} = 2s / \lambda$. We obtain $$\min_{z} \ \tilde{s} {\left\| Dz \right\|_{2}}^ 2 + {\left\| \tilde{y} - z \right\|_{2}}^2 .$$ This is a classical L2 spline and can be minimized using Equation (), as already explained. The optimal value of $d$ in Equation () can be explicitly computed using shrinkage operators, $$\begin{aligned}
    d^{k+1} &= {\underset{d}{\operatorname{argmin}}\ } \ {\left\| d \right\|_{1}} + \tfrac{\lambda}{2} {\left\| d - z^{k+1} + y - b^k \right\|_{2}}^2 \nonumber \\
    &= {\mathrm{Shrink}} (z^{k+1} - y + b^k , 1/\lambda) ,\end{aligned}$$ where $$\begin{aligned}
    {\mathrm{Shrink}}(v, \gamma) = 
    \begin{pmatrix}
    {\mathrm{shrink}}(v_1, \gamma) \\
    \vdots \\
    {\mathrm{shrink}}(v_j, \gamma) \\
    \vdots \\
    {\mathrm{shrink}}(v_m, \gamma)
    \end{pmatrix}
    \intertext{and}
    {\mathrm{shrink}}(x, \gamma) = \frac{x}{|x|} \max(|x|-\gamma, 0) .
    \label{eq:shrinkage}\end{aligned}$$ We thus obtain a very efficient algorithm for computing L1 splines, combining DCT and shrinkage operators.

On a different note let us mention that Equation () can also be interpreted [@mateos12] as a relaxation of $$\min_{z, d} {\left\| d \right\|_{0}} + s {\left\| Dz \right\|_{2}}^2 + \tfrac{\lambda}{2} {\left\| d-z+y \right\|_{2}}^2 ,
    \label{eq:L0formulation}$$ where the L0 norm is replaced by its (convex) L1 counterpart. In this case the underlying model for $y$ is $y = \hat{y} + r + d$, where $r$ is zero-mean Gaussian noise and $d$ represents the “oulier” noise. Under this assumptions, $d$ practically becomes an “indicator function” of the presence of (sparse) outliers (see also [@black91]). Besides the different angle in the derivation of the model, our approach differs from [@mateos12] in two very important points. First, we use split-Bregman iteration by introducing the variable $b^k$ in the optimization procedure, see equations () and (). In [@mateos12] Equation () is first solved using direct alternate minimization over $z$ and $d$, and then Equation () is solved via non-convex minimization using the previous solution as a starting point. Second, considering the grid structure, we use the DCT approach to solve Equation (), instead of the classical Cholesky decomposition. The combination between split-Bregman and DCT results in a sound and fast algorithm for computing L1 splines on grid data.

Handling missing data
---------------------

In the classical L2 formulation, a diagonal binary weighting matrix $W$ is used to cope with missing values (see Section \[sec:robustL2\] for details). Let us denote by $w$ the diagonal of $W$. Let us first define $$\begin{aligned}
    {\left\| z \right\|_{1, w}} &= \sum_{\substack{i=1 \\ w_i = 1}}^{m} |z_i|
    & \text{and} &&
    {\left\| z \right\|_{2, w}} &= \left( \sum_{\substack{i=1 \\ w_i = 1}}^{m} z_i^2 \right)^{1/2} .\end{aligned}$$ We then pose Equation () as $$\begin{gathered}
    \hat{y} = {\underset{z}{\operatorname{argmin}}\ } {\left\| (z - y) \right\|_{2, w}}^2 + s {\left\| Dz \right\|_{2}}^ 2 ,\end{gathered}$$ and equivalently extend Equation () as $$\begin{gathered}
    \hat{y} = {\underset{z}{\operatorname{argmin}}\ } {\left\| (z - y) \right\|_{1, w}} + s {\left\| Dz \right\|_{2}}^ 2 .\end{gathered}$$ This leads to $$\min_{z, d} \ {\left\| d \right\|_{1, w}} + s {\left\| Dz \right\|_{2}}^ 2 + \tfrac{\lambda}{2} {\left\| d - z + y \right\|_{2, w}}^2 .$$ Then the split-Bregman iteration can be written as $$\begin{aligned}
    z^{k+1} &= {\underset{z}{\operatorname{argmin}}\ } \, s {\left\| Dz \right\|_{2}}^ 2 + \tfrac{\lambda}{2} {\left\| d^k - z + y - b^k \right\|_{2,w}}^2 , \label{eq:WsplitBregman1_1} \\
    d^{k+1} &= {\underset{d}{\operatorname{argmin}}\ } \, {\left\| d \right\|_{1,w}} + \tfrac{\lambda}{2} {\left\| d - z^{k+1} + y - b^k \right\|_{2,w}}^2 , \label{eq:WsplitBregman1_2} \\
    b^{k+1} &= b^k + (z^{k+1} - y - d^{k+1}) . \label{eq:WsplitBregman2}\end{aligned}$$ Equation () can be solved using Equation (). Solving Equation () amounts to performing a shrinkage operation on the dimensions where $w$ equals 1. In Equation (), it suffices to update the dimensions of $b^k$ where $w$ equals 1.

Handling multidimensional data
------------------------------

Let us now return to the general case of $m$-dimensional data. Following Garcia [@garcia10], we extend Equation () as $$\hat{y} = {\mathrm{DCT}}^{-1}_m \left( \varGamma^m \circ {\mathrm{DCT}}_m (y) \right) ,
    \label{eq:LSsolutionDCTd}$$ where ${\mathrm{DCT}}_m(\cdot)$ and ${\mathrm{DCT}}^{-1}_m(\cdot)$ stand for the $m$-dimensional DCT-II and inverse DCT-II functions, and $\circ$ denotes the Schur (element-wise) product. Notice that the multidimensional DCT is simply a composition of one-dimensional DCTs along each dimension. Extending Equation (), $\varGamma^m$ is an $m$-th order tensor defined by $$\varGamma^m = 1^m  \div \left( 1^m + s \varLambda^m \circ \varLambda^m \right) ,
    \label{eq:varGammad}$$ where $1^m$ is an $m$-th order tensor of ones, and $\div$ denotes the element-wise division. Finally, $\varLambda^m$ is an $m$-th order tensor, defined by $$\varLambda^m_{i_1, \dots, i_m} = \sum_{j=1}^{d} \left(  -2 + 2 \cos \frac{(i_j - 1) \pi}{n_j} \right) .
    \label{eq:varLambdad}$$ where $n_j$ denotes the size of $\varLambda^m$ along the $j$-th dimension.

**The algorithm and its complexity.** The pseudocode for the general $m$-dimensional case is presented in Algorithm \[algo:l1spline\]. Let us analyze its complexity. The DCT and inverse DCT require $O(n \log n)$ operations, where $n=\prod_{1 \leq j \leq m} \, n_j$. The remaining operations are linear in $m$. The overall complexity of the algorithm is then $O \left( N_o N_i (m + n \log n) \right)$, where $N_o$ and $N_i$ are, respectively, the number of outer-loop and inner-loop iterations in Algorithm \[algo:l1spline\]. Notice that Goldstein and Osher [@goldstein09] recommend to perform only one inner-loop iteration for achieving optimal efficiency. Thus, we set $N_i = 1$ for all experiments. We will later see that in many cases the algorithm converges quickly ($N_o$ can be very small then). The algorithm’s complexity is thus dominated by the computation of the DCT and inverse DCT. Of course, these standard operations can be easily computed using GPU, speeding-up the execution by several orders of magnitude.

compute $\varGamma^m$ according to equations () and (). $d^1 \gets 0$, $b^1 \gets 0$, $k \gets 1$

$z^{k+1} \gets {\mathrm{DCT}}^{-1}_m (\varGamma^m \circ {\mathrm{DCT}}_m (d^k + y - b^k))$ $\displaystyle d^{k+1} \gets {\mathrm{Shrink}} (z^{k+1} - y + b^k , 1/ \lambda)$ $b^{k+1} \gets b^k + (z^{k+1} - y - d^{k+1})$ $k \gets k+1$ $z^k$

Experimental results {#sec:results}
====================

For all experiments we adhere to the following setup:

1.  using generalized cross validation, we find the best estimate $\hat{s}$ for $s$ for the robust L2 formulation (problem ());

2.  we then find L2 splines (problem ()), robust L2 splines (problem ()),[^2] and/or L1 splines (problem ()), setting $s = \hat{s}$.

This protocol allows us to show that, even when $s$ is chosen to fit optimally the robust L2 formulation, the proposed method provides better estimates. For the L1 formulation, in Equation (), we simply set $\lambda=\min(s,1)$ for all examples. We recall that $N_i$ is set to 1. We also set $\varepsilon = 10^{-3}$ (see Algorithm \[algo:l1spline\]) and additionally limit the maximum number of outer iterations to a hundred. The algorithm stops when any of the two conditions is met.

Fig. \[fig:spline\_1D\_16\] presents two one-dimensional examples. We depict the original signal $\hat{y} \in [1, \dots, n] \rightarrow {{\mathbb R}}$, where $n=2^{16}$. We observe the signal $y = \hat{y} + r_1$, where $r_1$ is Gaussian noise. Some points $y_j$ ($1 \leq j \leq n$) are further contaminated with uniform noise $r_2$, where $r_2 \in [ a, \dots, b ]$, such that $y_j = \min (\max (\hat{y}_j + r_1 + r_2, a), b)$. The points affected by $r_2$ are depicted in red and the remaining ones in green. In the top row, $a = -5, b = 5$; and in the bottom row, $a = 0, b = 5$. In both examples, only the L1 spline is correct. The classical and robust L2 splines are both unable to correctly recover the original data in the corrupted part.

Fig. \[fig:errorPlot\] shows the evolution of the relative error as the number of iterations increases for the example in Fig. \[fig:spline\_1D\_16\] (top row). As we can observe, the proposed algorithm is able to converge quickly, reaching a precision of $10^{-3}$ in less than 20 iterations.

![The relative error (in logarithmic scale) as the iterations progress when computing the first example in Figure \[fig:spline\_1D\_16\]. The algorithm is able to quickly decrease the error during the first 20 iterations.[]{data-label="fig:errorPlot"}](errorPlot.pdf){width=".5\columnwidth"}

Fig. \[fig:timeDistribution\] depicts the relative time-cost of each operation during the execution of the proposed method (Fig. \[fig:spline\_1D\_16\], top row). Computing the DCT and the inverse DCT covers more than 84% of the total running time. Implementing these standard operations in GPU would boost the performance of the algorithm by orders of magnitude.

![Percentage of the execution time spent in each operation when computing the first example in Figure \[fig:spline\_1D\_16\]. Clearly, the vast majority of time is spent in DCT or inverse DCT operations.[]{data-label="fig:timeDistribution"}](timeChart.pdf){width=".5\columnwidth"}

In Fig. \[fig:spline\_2D\_n30\] we present a two-dimensional example. We depict the original signal $\hat{y} \in [1, \dots, 256]^2\rightarrow [-6.5497, 8.1054]$ in Fig. \[fig:spline\_2D\_n30\_original\], and we add two types of noise: first, Gaussian noise $r_1$ with zero mean and variance $\sigma^2=2$ (Fig. \[fig:spline\_2D\_n30\_noisy\]), and then uniform noise $r_2$ in the interval $[ -5 \cdot \max (\hat{y} + r_1), \dots, 5 \cdot \max (\hat{y} + r_1) ]$ (Fig. \[fig:spline\_2D\_n30\_corrupted\]). Again, only the L1 spline correctly recovers the original signal.

  -- -- --
        
        
  -- -- --

We next test the proposed algorithm with a climate time-series provided by the Met Office Hadley Centre [@brohan06].[^3] The dataset contains the evolution of global average land temperature anomaly (in ) with respect to the 1961-1990 average temperature. The results, which confirm an upward trend in the second half of the 20th century, are shown in Fig. \[fig:temp\].

We also test on this dataset the effect of varying the parameters $s$ and $\lambda$, see Fig. \[fig:params\]. As in the classical L2 formulation, $s$ has a direct impact on the obtained result, controlling the degree of smoothness of the solution, see Fig. \[fig:params\_s\]. On the contrary, Fig. \[fig:params\_lambda\] shows that the newly introduced parameter $\lambda$ is very stable and provides very similar results in a wide range $(\lambda \in [0.1, 100])$. This stability allows us to fix its value to $\lambda=1$ for all the experiments in this work.

Another interesting example is presented in [@mateos12]. The dataset consists of power consumption measurements (in kW) for a government building, collected every fifteen minutes from July 2005 to October 2010. As in [@mateos12], we downsample the data by a factor of four, yielding one measurement per hour, and use only a subset of the whole data. The results are displayed in Fig. \[fig:load\].

![Power consumption measurements (in kW) for a government building [@mateos12] with zoom-in details on the bottom.[]{data-label="fig:load"}](spline_1D_load-a.pdf "fig:"){width=".2\columnwidth"} ![Power consumption measurements (in kW) for a government building [@mateos12] with zoom-in details on the bottom.[]{data-label="fig:load"}](spline_1D_load-b.pdf "fig:"){width=".2\columnwidth"} ![Power consumption measurements (in kW) for a government building [@mateos12] with zoom-in details on the bottom.[]{data-label="fig:load"}](spline_1D_load-c.pdf "fig:"){width=".2\columnwidth"}

We also test in a synthetic example the ability to recover signals with sharp transitions, see Fig. \[fig:square\]. In this case we use a simple piece-wise constant function. We can observe clear overshoot (plus ringing) effects on the L2 and robust L2 splines. The robust L2 spline also results in transitions with less vertical slopes, creating a bluring effect. With the L1 spline we obtain a much better reconstruction, with almost non-existent overshooting.

This very same effect can be observed in real examples, see Fig. \[fig:119082\]. When approximating images with splines, some structure is lost by blur and some structure is artificially created by overshooting and ringing. This can be observed in Figs. \[fig:119082\_original\] and \[fig:119082\_noise\], were the difference between the original image and the image estimated by robust L2 splines exhibits structure. Observe, however, that almost no structure in the difference is visible when the reconstruction is performed using L1 splines.

[@c@c@]{}

Application to range data
-------------------------

In this section we perform smoothing of depth data obtained with a Kinect camera. This kind of data is particularly challenging because:

-   it presents relatively smooth areas separated by sharp transitions,

-   edges are highly noisy, that is, edge pixels oscillate over time between foreground and background, and

-   it contains missing data, which appear for two different reasons: (1) the disparity between the IR projector and the IR camera produces “shadows,” and (2) the depth cannot be recovered in areas where the IR pattern is not clearly observable (e.g., because they receive direct sunlight or interference from another Kinect).

We use splines to interpolate and denoise these data, showing the advantage of L1 splines over its robust L2 counterpart. The displayed images are part of the LIRIS human activities dataset [@harl2012].

In the first example, shown in Fig. \[fig:kinect2D\], we use a single depth frame (with standard Kinect resolution of $640 \times 480$). The missing data are represented in black, while depth data goes from red to yellow as depth increases. Both, the L1 spline and the robust L2 spline are able to interpolate the missing data with reasonable values. Notice, however, that the latter exhibits, as aforementioned, overshooting and ringing (clearly perceived in the 1D profile). These effects are much milder in the L1 reconstruction.

[ccc]{}

In Fig. \[fig:kinect3D\_original\] we can clearly observe that the position of the missing data is not consistent across frames. We can integrate data from several frames to achieve more accurate interpolations, by performing 3D reconstructions. Thus, in this example, we treat depth data as a 3D signal (2D + time), by considering three consecutive frames. The data dimensionality is then $640 \times 480 \times 3$. A full depth video can be smoothed by using 3D splines as a sliding-window type of filter. The robust L2 spline again presents a noisier behavior and with significative overshooting. On the other hand, the L1 spline is much smoother in smooth areas while correctly preserving abrupt transitions.

**Running times.** We present in Table \[tab:runningTime\] the running-time and number of iterations until convergence for every example in this work. The time of the robust L2 and the L1 splines is comparable. All code is written in pure Matlab, with no C++ or mex optimizations. All experiments were run on a MacBook Pro with a 2.7GHz Intel Core i7 processor. Finally, note that in most cases the algorithm converges in less than twenty outer-iterations (recall that $\varepsilon=10^{-3}$). In the example in Fig. \[fig:spline\_2D\_n30\], the maximum number of iterations (100) is reached with a final error of $10^{-2.8}$.

  -- ------------------------------- --------------------------- ------------------ ------- --------
                                                                  Robust L2 spline          
                                                                        Time         Time    Iters.
     Fig. \[fig:spline\_1D\_16\]              $2^{20}$                 4.931         3.590     7
     Fig. \[fig:spline\_2D\_n30\]         $256 \times 256$             0.541         0.982    100
     Fig. \[fig:temp\]                           163                   0.007         0.045     72
     Fig. \[fig:load\]                           501                   0.010         0.010     11
     Fig. \[fig:square\]                      $2^{10}$                 0.035         0.007     6
     Fig. \[fig:119082\_original\]        $321 \times 481$             1.167         0.758     17
     Fig. \[fig:119082\_noise\]           $321 \times 481$             1.143         0.873     20
     Fig. \[fig:kinect2D\]                $480 \times 640$             0.911         0.266     2
     Fig. \[fig:kinect3D\]            $480 \times 640 \times 3$        7.511         2.202     3
  -- ------------------------------- --------------------------- ------------------ ------- --------

  : Execution times (in seconds) and number of iterations until convergence of the proposed algorithm for the different experiments performed in this work.[]{data-label="tab:runningTime"}

Conclusions {#sec:conclusions}
===========

We have presented a new method for robustly smoothing regularly sampled data. We do this with modified splines, where we replace the classical L2-norm in the fitting term by an L1-norm. This automatically handles outliers, thus obtaining a robust approximation.

We also presented a new technique, using split-Bregman iteration, for solving the resulting optimization problem. The algorithm is extremely simple and easy to code. The method converges very quickly and has a small memory footprint. It also makes extensive use of the DCT, thus being straightforward to implement in GPU. These characteristics make this method very suitable for large-scale problems.

Acknowledgment {#acknowledgment .unnumbered}
==============

Work partially supported by NSF, ONR, NGA, ARO, DARPA, and NSSEFF. We thank Dr. Gonzalo Mateos for kindly providing the power consumption dataset.

[^1]: This work was partially done while the authors were with the Department of Electrical and Computer Engineering, University of Minnesota.

[^2]: Code available at <http://www.mathworks.com/matlabcentral/fileexchange/25634-robust-spline-smoothing-for-1-d-to-n-d-data>.

[^3]: Data are available in <http://hadobs.metoffice.com/crutem3/diagnostics/global/nh+sh/annual>.
