RobustTrend: A Huber Loss with a Combined First and Second Order Difference Regularization for Time Series Trend Filtering

1 Introduction

The explosion of time series data generated by the Internet of Things (IoT) and many other sources has made the time series analytics more challenging. Trend filtering, which estimates the underlying trend from the time series, is an important task arising in a variety of real-world applications [Kim et al., 2009; Craigmile and Percival, 2006]. Specifically, we assume the time series consists of the trend component and a more rapidly changing random component. The goal of trend filtering is to extract the trend component accurately and efficiently. The extracted trend component can reveal certain traits of a time series and is crucial in many applications such as time series forecasting and anomaly detection. [?]

Many algorithms have been proposed for trend filtering, including moving average smoother, kernel smoothing, smoothing splines etc [Wu et al., 2007; Wu and Huang, 2009; Fried et al., 2018; Afanasyev and Fedorova, 2016]. Among them, one of the most classical and widely used trend filters is the Hodrick-Prescott (H-P) filtering [Hodrick and Prescott, 1997] and its variants. In the H-P trend filtering, the trend component is obtained by solving an optimization problem which minimizes the weighted sum of the residual size and the smoothness of the trend component. On the one hand, the extracted trend component is expected to close to the original time series, but meanwhile the smoothness constraint is imposed. Specifically, the sum-of-squares loss function is used to measure the difference between the trend and the original time series, and the smoothness of the trend component is measured using the second order difference of the trend. The drawback of H-P filtering is that it is not robust to outliers. Based on the H-P filtering, the ℓ1 trend filtering [Kim et al., 2009] is proposed, which basically substitutes the ℓ1norm for the ℓ2-norm in the regularization term in H-P trend filtering. Note that the ℓ1 trend filtering assumes that the underlying trend is piecewise linear, and the kinks, knots, or changes in slope of the estimated trend can be interpreted as abrupt changes or events.

Although the sum-of-squares loss function in both the H-P trend filtering and the ℓ1 trend filtering enjoys its popularity, it cannot handle heavy tailed error or outliers. One approach to deal with these situations is the least absolute deviation (LAD) loss [Wang et al., 2007] as it mitigates the loss incurred by big errors. However, the LAD loss is not adapted for small errors since it penalizes strongly on the small residuals. In the extreme case where the error is not heavy tailed or there is no outliers, it is less effective than the sum-ofsquares loss function. In practice, generally we expect the loss function can adaptively handle these cases. The Huber loss [Huber and Ronchetti, 2009] is a combination of the sumof-squares loss and the LAD loss, which is quadratic on small errors but grows linearly for large values of errors. As a result, the Huber loss is not only more robust against outliers but also more adaptive for different types of data.

In this paper, besides adopting Huber loss, we propose to use the combination of the first order and the second order difference regularization of the trend component in trend filtering. The regularization based on the first order difference encourages to preserve abrupt trend change, but meanwhile it tends to introduce the staircases for the extracted trend component when slow trend change exists. One way to reducing the staircasing effect is the introduction of higher order difference in the regularizier. In fact, the second order difference can effectively eliminates the artifacts (e.g., staircasing) introduced by the first regularization without smoothing the time series strongly. Also, it is proved recently that the introduction of the second order difference can achieve comparable results to higher order difference [Papafitsoros and Sch¨onlieb, 2014]. Thus, we consider the combination of the first and second order difference for the best trade-off.

To sum up, in this paper we propose a novel trend filtering method based on the Huber loss and the combination of the first and second order dynamics of the trend component. The resulting optimization problem is solved via an efficient MMADMM algorithm. The complexity of the Huber loss makes it challenging to directly apply the ADMM framework [Boyd et al., 2011]. Inspired by the majorization minimization [Sun et al., 2017], we minimize an upper bound of the Huber loss in the ADMM framework. In the experiments, we compare ten different trend filtering algorithms and results on both synthetic and real-world data sets demonstrate the good performance of our trend filtering algorithm. To the best of our knowledge, our empirical comparisons is one of the most extensive and comprehensive studies which includes common popular trend filtering methods, including H-P trend filter, ℓ1 trend filter, TV denoising filter [Chan et al., 2001], mixed trend filter [Tibshirani, 2014], wavelet trend filter [Craigmile and Percival, 2006], repeated median filter [Siegel, 1982], robfilter [Fried, 2004], EMD trend filter [Moghtaderi et al., 2011], and so on.

The rest of the paper is organized as follows. In Section 2, we briefly review the trend filtering and some related work. In Section 3, We introduce our robust trend filter and our efficient algorithm. We present the empirical studies of the proposed method in comparison with other state-of-the-art algorithms in Section 4, and conclude the paper in Section 5.

2 RelatedWork

There are various trend filtering methods proposed in the literature. As discussed in the previous section, the H-P trend filtering [Hodrick and Prescott, 1997] and ℓ1 trend filtering [Kim et al., 2009] extract the trend component by minimizing the weighted sum of the residual size and the smoothness of the trend component using the second order difference operator. The TV denoising filtering [Chan et al., 2001] adopts the same loss function as the H-P and ℓ1 trend filtering, but it utilizes the first order difference as the regularization to deal with abrupt trend changes. As a result, the TV denoising filtering can handle the abrupt trend change successfully, but it tends to introduce staircases for slow trend changes. In contrast, the H-P and ℓ1 trend filtering generally prone to trend delay for abrupt trend changes. The mixed trend filtering [Tibshirani, 2014] adopts two difference orders as regularization. It is assumed that the observed time series were drawn from an underlying function possessing different orders of piecewise polynomial smoothness. However, it is difficult to test this assumption adaptively. Thus, the selection of the two order remains a challenging problem, and limits its usage in practice. Besides, due to the sum-of-squares loss, all the aforementioned trend filtering methods cannot handle heavy tailed errors and outliers.

In [Fried, 2004; Fried et al., 2018], a robust trend filtering called robfilter is proposed. In robfilter, a robust regression functional for local approximation of the trend is used in a moving window. To further improve the robustness of robfilter, outliers in the original time series is removed by trimming or winsorization based on robust scale estimators. Empirical mode decomposition (EMD) [Wu et al., 2007] is a popular method for analyzing non-stationary time series. With iterative process, it finds local maxima and minima, and then interpolates them into upper and lower envelopes. The mean of the envelopes yields local trend, while the residual must comply two requirements to become an intrinsic mode function (IMF). This method assumes that the trend is "slowly varying" and the selection of "slowest" IMF sums up to be the trend. Since the local maxima and minima are sensitive to noise and outliers, the conventional EMD are not robust. Another variant, Ensemble EMD (EEMD), is a noiseassisted method. It is an ensemble of white noise-added data and treats mean of their output as the final true result. It resolves the mode mixing issue in original EMD but still suffer from outliers. The wavelet analysis [Afanasyev and Fedorova, 2016] is a transformation of the original time series into two types of coefficients, wavelet coefficients and scaling coefficients. Note that the wavelet coefficients are related to changes of average over specific scales, whereas scaling coefficients can be associated with averages on a specified scale. Since the scale associated with scaling coefficients is usually fairly large, the general idea behind the wavelet trend filter [Craigmile and Percival, 2006] is to associate the scaling coefficients with the trend, and the wavelet coefficients with the noise component. The benefits of trend filtering via wavelet transformation is that it can capture the abrupt change of trend, but how to choose orthogonal basis still remains a challenge.

In this paper, we focus on trend filtering for time series without distinguishing periodic and trend components. When the seasonality/periodic and trend components need to be separated, the seasonal-trend decomposition method can be utilized, such as STL [Cleveland et al., 1990], RobustSTL [Wen et al., 2019], and TBATS [De Livera et al., 2011].

3 Proposed Robust Trend Filter

3.1 Model Overview

We consider the time series of length N as y = [y0, y1, · · · , yN-1]T , which can be decomposed into trend and residual components [Alexandrov et al., 2012]

where τ = [τ0, τ1, · · · , τN-1]T denotes the trend component, and r = [r0, r1, · · · , rN-1]T denotes the residual component. Note that in this paper we consider the time series faced in real-world applications. Therefore, the trend component may contains both slow trend change and abrupt trend change, while the residual component may contains both noise and outliers.

3.2 Preliminary

In the H-P trend filtering [Hodrick and Prescott, 1997], where the trend estimation τt is chosen to minimize the objective function as follows

where the loss function (first term) measures the size of the residual signal after trend extraction while the regularization (second term) measures the smoothness of the extracted trend, and the regularization parameter λ > 0 controls the trade-off between them. The Eq. (2) can be equivalently formulated in the matrix form as

where D(2) is the second-order difference matrix, and its definition is given in Eq. (6) below.

By replacing the ℓ2-norm with ℓ1-norm in the regularization term of H-P trend filtering, we obtain the ℓ1 trend filtering [Kim et al., 2009] as follows:

The ℓ1 regularization in Eq. (4) promotes sparsity, which indicates ||D(2)τ ||1 would have many zero elements. As a result, the ℓ1 trend filter leads to a piecewise linear trend component.

As an extension of Eq. (4), piecewise polynomial trend signal can be obtained by using high-order difference operator in the regularization term [Kim et al., 2009], i.e.,

where D(k+1) ∈ R(n-k-1)×n is the discrete difference operator of order k + 1, which can be defined recursively as D(k+1) = D(1) · D(k) and the D(1) and D(2) are the firstorder and second-order difference, i.e.,

When k = 0 in Eq. (5), it corresponds to the total variation (TV) filter [Chan et al., 2001], which is suitable when the underlying trend signal is approximately piecewise constant.

3.3 Robust Trend Filtering with Sparse Model

All these trend filtering methods discussed in Section 3.2 adopt the sum-of-squares loss, which is appropriate when the residual part is near to Gaussian noise. However, in realworld applications, the residual often exhibits long-tailed distribution due to outliers. In this paper, we adopt the Huber loss from robust statistics [Huber and Ronchetti, 2009] to deal with possible long-tailed distributed residual. As discussed in Section 1, Huber loss combines the advantages of the sum-ofsquares loss and the LAD loss adaptively.

In order to simultaneously capture abrupt and slow trend changes, we consider two sparse regularization terms like the fussed LASSO [Tibshirani et al., 2005]. Specifically, we utilize first and second order difference operators in the regularization to deal with abrupt and slow trend changes, respectively. Since it is proved that the combination of the first and second order difference regularization can achieve comparable results to higher order difference regularization [Papafitsoros and Sch¨onlieb, 2014], in this paper we do not consider other higher order difference operators.

Based on the above discussion, we propose a Huber loss with a combined first and second order difference regularization as a robust trend filtering for time series (named RobustTrend in this paper), i.e.,

where D(1) , D(2) are the first-order and second-order difference matrix, respectively. And gγ(x) = ∑i gγ(xi) is the summation of elementwise Huber loss function with each element as

and its derivative can be obtained as

where sgn is the sign function.

3.4 Efficient MM-ADMMAlgorithm

In this section, we develop a specialized ADMM algorithm to solve the minimization problem of RobustTrend in Eq. (7). ADMM [Boyd et al., 2011] is a powerful framework in optimization and also utilized in trend filtering with sum-ofsquare loss [Ramdas and Tibshirani, 2016]. However, the standard ADMM is not efficient in the RobustTrend filter since no exact solution exists in the updating step with Huber loss (will be shown later). One may transform Eq. (7) into a series of ordinary fused lasso problems as in [Polson and Scott, 2016]. But it is still inefficient since each ordinary fused lasso problem needs to be solved iteratively as well. To alleviate this problem, we utilize the majorizationminimization (MM) framework [Sun et al., 2017] to approximate the Huber loss function in the ADMM updating step to obtain a closed form for efficient computation.

To simplify the ADMM formulation for Eq. (7), we rewrite it in an equivalent form as

where

Then we can obtain the augmented Lagrangian as

where v is the dual variable, and ρ is the penalty parameter. Following the ADMM procedure in [Boyd et al., 2011], we can obtain the updating steps as

where u = (1/ρ)v is the scaled dual variable to make the formulation more convenient. The above z-minimization step can be efficiently solved by soft thresholding as

where the soft thresholding operator Sρ(a) is defined as

The τ-minimization step in Eq. (12) is nontrivial since no exact solution exists. One may use iterative algorithms, like standard gradient methods or conjugate gradient methods, to solve the τ-minimization step since it is convex. However, this would increase the overall computation since the ADMM itself is an iterative algorithm. Inspired by the work of [Eckstein and Bertsekas, 1992] that the ADMM can still converge when the minimization steps are carried out approximately, we utilize one-iteration MM framework like [Pham and Venkatesh, 2013] to efficiently solve the τ-minimization step in Eq. (12) with closed form solution. Specifically, we adopt the proved sharpest quadratic majorization for Huber loss as in [de Leeuw and Lange, 2009] to minimize the difference between the exact and approximate solution, i.e.,

where ηγ(xi|xki ) denotes the majorization of element xi's Huber loss and xki is the solution of xi at ADMM's kth iteration. Then, the summation of Huber loss can be obtained as

where C represents a constant number and A is a diagonal matrix as

Next, let define x = y - τ and replace the Huber loss in Eq. (12) with the quadratic majorization in Eq. (17), we can obtain an approximate τ-minimization step with oneiteration MM algorithm as

To obtain the closed form of the above approximate τ-minimization step, we take derivative of Eq. (19), set it to zero to obtain xk+1 , and put xk+1 into Eq. (20) to get the final τ-minimization step as

Therefore, to obtain the extracted trend signal of the proposed RobustTrend filter, our designed MM-ADMM algorithm sequentially calculates Eq. (21), Eq. (15), and Eq. (14) in each ADMM iteration until termination.

For the termination criteria, besides the maximum iteration number, we also check the primal and dual residuals in each ADMM iteration until they are small enough. Specifically, the values of primal residual rk and dual residual rk at kth iteration are calculated by

Then, ‖rk ‖2 and ‖sk ‖2 are checked if they are smaller than the corresponding tolerance thresholds as [Boyd et al., 2011]

where ϵabs > 0 is an absolute tolerance and ϵrel > 0 is a relative tolerance. A reasonable choice for the two tolerance is around 10-3 to 10-4 .

3.5 Online Extension for Data Streams

In order to support streaming data applications, the proposed RobustTrend filter is extended to online mode by applying sliding window scheme like robfilter and repeated median filter in [Fried, 2004; Fried et al., 2018]. During each sliding window, the trend signal is extracted using our designed MMADMM algorithm described in Section 3.4. To speed up, the ρDT D and Dy are computed once and cached for the τ-minimization step in Eq. (21). Furthermore, warm start is adopted in each new sliding window, where the initial values for τ , z, and u are from previous window's final solution by removing the first element and copying last element once. Since the trend signal usually exhibits slow change, this warm start trick often provides an appropriate approximation to result in fewer ADMM iterations.

4 Experiments and Comparisons

In this section, we conduct experiments on both synthetic data and real-world data to demonstrate the effectiveness of the proposed RobustTrend filtering algorithm.

4.1 Baseline Algorithms

We tested nine state-of-the-art trend filters as baseline algorithms to perform a relative comprehensive evaluation, including H-P filter, ℓ1 trend filter, TV denoising filter, Mixed trend filter [Tibshirani, 2014], Wavelet trend filter [Craigmile and Percival, 2006], Repeated median filter [Siegel, 1982], robfilter [Fried, 2004], EMD trend filter [Moghtaderi et al., 2011] and its extension EEMD trend filter.

4.2 Dataset

For synthetic data, we first generate the trend signal with 1000 data points, which contains a sin wave, a triangle wave, and a square wave with 1.0 amplitude to represent smooth slow trend change as well as abrupt trend change with sharp discontinuity. Next, we add Gaussian noise with 0.2 standard deviation. Then, we add 10 to 200 spikes and dips with 2.0 amplitude to represent 1% to 20% outlier ratio. As an example, the trend signal and trend with noise and 1% outlier are shown in Figure 1(a).

Note that we generate synthetic data with sin, triangle, and square wave as trend signal to simulate different types of trend changes in real-world data. If the trend signal only contains triangle/rectangle shape, the ℓ1/TV filter would be the optimal choice (under the case of white noise without outliers). However, these specific assumptions are of limited usage in practice. As the ground truth of trend is known, the mean squared error (MSE) and mean absolute error (MAE) are used for the performance evaluation later.

In addition to the synthetic data, we also adopt one realworld data from Yahoo's Anomaly Detection Dataset 1 , i.e., real 31 in A1Benchmark, which is shown in Figure 1(b). It can be seen that there is abrupt trend change around point 1220 and a outlier at point 1306.

4.3 Experiment Results on Synthetic Data

Figure 2 shows the trend filtering results from different algorithms when there is 1% outlier ratio. As can be seen in the figure, ℓ1 trend filter, wavelet trend filter, and EEMD trend filter are heavily affected by the abrupt level change around the square wave. TV denoising filter captures abrupt level change but exhibits staircasing around sine and triangle wave. Due to space limitation, we omit the results of H-P and EMD trend filters since their performances are worse than ℓ1 and EEMD trend filters, respectively. To better evaluate the performance quantitatively, we calculate MSE and MAE to quantify the extraction accuracy. In addition to extracting trends on 1% of outliers, we also increase the ratio of outliers to 5%, 10%, and 20%. Table 1 summarizes the performance of our proposed algorithm along with other nine state-of-the-art trending filtering algorithms. The best results for each case are highlighted in bold fonts. Overall, it can be observed that our algorithm outperforms others.

1 https://webscope.sandbox.yahoo.com/catalog.php?datatype=s

The trend recovery near change points is important as it measures how prompt we can capture the trend change in time. Thus, we calculate MSE and MAE around 9 change points and their neighbors (one before and one after each change point) in the synthetic data (total 27 points) when outlier ratio is 5%. The results are summarized in Table 2 where the best results are highlighted in bold fonts. Clearly, our algorithm achieves significantly more accurate trend near change points.

To evaluate the different components of our RobustTrend filter, we also compare TV denoising filter with Huber loss (i.e., RobustTrend filter without second order difference regularization), ℓ1 trend filter with Huber loss (i.e., RobustTrend filter without first order difference regularization), and RobustTrend filter with ℓ2-norm regularizations. The results are summarized in Table 3 with the outlier ratio setting to 5%, where the best results are highlighted in bold fonts. It can be seen that the Huber loss and the first order and second order difference ℓ1 regularization terms bring significant performance improvements in trend filtering.

4.4 Experiment Results on Real-World Data

We perform experiment on one real-world dataset from Yahoo depicted in Figure 1(b). In this experiment, we apply the online mode of the trend filtering to investigate how it performs on streaming data. The results of top-3 performance filters (our RobustTrend, robfilter, and repeated median filter) and the popular ℓ1 trend filter are summarized in Figure 3. Since the beginning of the data is almost constant, we only show the zoomed-in results after point 1200. It can be observed that our RobustTrend is able to follow the sudden increase of the trend and is not affected by the outlier. The ℓ1 filter and robfilter show a delay in capturing the trend change. The ℓ1 filter is also sensitive to the outliers in the original signal, leading to a dip in the trend. The repeated median algorithm overshoots the estimation of the trend and generates a higher trend than the actual one near the position where trend abruptly changes.

5 Conclusions

In this paper we propose a robust trend filtering method RobustTrend based on the Huber loss. The Huber loss can adaptively deal with residual outliers and different types of data. In order to impose the smoothness of the trend and meanwhile capture the abrupt trend change, we propose to use the combination of the first order and the second order difference of the trend component as regularization in trend filtering. To efficiently solve the resulting optimization problem, we design an MM-ADMM algorithm which applies majorizationminimization to facilitate the updating step for the Huber loss effectively.