4
$\begingroup$

I want to use a Savitzky-Golay filter to smooth some data. There is a right width to use based on the data that it is smoothing. A number of papers basically use "eyeball norm" on the parameters but that feels like voodoo.

How do I get a measure like AIC to tell me which implementation is the proper one to use?

When I think about AIC, I do so in terms of RSS, number parameters, and number of degrees of freedom. I can measure RSS for the smoothed data, but I don't know how to think about parameter count or number of degrees of freedom in this context. Hints or pointers would be appreciated. Clear explanation would be delightful.

I'm personally using either R, or LabVIEW, but the actual software package isn't important.

Work so far:

  • This post says "k is 2 parameters, n is samples, use RSS" without providing references or derivation.
  • This reference is about the 'PoMoS' library which claimes to use genetic algorithms and information criteria to find optimal polynomial structure of time-series.

References:

$\endgroup$
2
  • 1
    $\begingroup$ It may be worth noting that the AIC is often viewed as a proxy for Kullback-Liebler divergence. $\endgroup$
    – user78229
    Commented Dec 6, 2016 at 16:18
  • $\begingroup$ @DJohnson - AIC is derived from Kullback-Liebler, with assumptions. $\endgroup$ Commented Dec 6, 2016 at 17:28

3 Answers 3

4
$\begingroup$

The following findings are not restricted to the Savitzky-Golay filter, but can easily applied to other Finite Impulse Response (FIR) filters like the Modified Sinc Smoother proposed in (1) in TL;DR. I will first sketch the guideline for the most exact approach that also includes the boundaries, but only derive simplifications from it that I found to still work well in practice.

TL;DR

When viewing the non-parametric smoothing process as a linear projection, the $AICc$ (corrected for potentially small sample sizes) of the Savitzky-Golay filter (which should be replaced by newer filters (1)) can be approximated by

$AICc=AIC+\frac{2\cdot k\cdot\left(k+1\right)}{n-k-1}$

from

$AIC=n\cdot\ln\left(\frac{1}{n}\sum_{i=1}^{n}r_{i}^2\right)+2\cdot k$

$k=\left(p_{eff}+1\right)$

$p_{eff}\approx n\cdot c_{0}$

where

  • $n$ is the signal's size,
  • $c_{0}$ is the central coefficient of the filter (the one applied to the data point in the exact center of the moving window)
  • $r_{i}$ is the residual at the $i$-th data point, i.e., the difference between the original and the smoothed signal
  • $p_{eff}$ is the effective number of parameters, i.e, the generalization of the number of parameters for Linear Least Squares

This is based on the theory of non-parametric smoothing and thus not a subjective choice.
Cross-Validation is also possible in a similar fashion and besides that, the findings even justify eye-balling to a certain extent.
As an outlook, the results can be generalized easily to more advanced Finite Impulse Response (FIR) filter techniques.
Since the projection or hat matrix is "derived", outlier handling becomes accessible as an added bonus.
Python implementations are provided in the Appendix.

⚠️ Optimization of the filter parameters has to be done for the filter of derivative order 0! The optimized parameters can then be used for derivative computation, but only then! ⚠️

(1)see Schmid, Rath & Diebold, Why and How Savitzky–Golay Filters Should Be Replaced, ACS Meas. Sci. Au 2022, 2, 2, 185–196, DOI: 10.1021/acsmeasuresciau.1c00054

1) FIR-Smoothing as a linear projection

1.1) Central Part

Basically, the application of the Savitzky-Golay or any other Finite Impulse Response (FIR) filter is straightforward. The filter is aligned with the data and the filter coefficients are multiplied element-wise with the data. Assigning the sum of these products to the central data point in exact center of the window gives the smoothed data point. After this, the filter slides one step to the right and the process is repeated.

Of course, this process only works in the central part where the filter overlaps completely with the data. At the beginning and the end of the data, the filtering has to be handled differently, but this will be described later.

Actually, this process is better known as a so-called discrete convolution.
With a little linear algebra, it can be written as a matrix vector multiplication, i.e., as linear projection:

$\mathbf{\hat{y}}=\mathbf{H}\cdot\mathbf{y}$

with

  • $\mathbf{y}$ as the original data of length $n$
  • $\mathbf{\hat{y}}$ as the smoothed data of length $n$
  • $\mathbf{H}$ as the projection or hat matrix of shape $n\times n$ ("puts the hat on y")

For now, we will only consider the central part. With the window size $w$, the projection for the central part can be written as:

$\mathbf{\hat{y}}_{central}=\mathbf{H}_{central}\cdot\mathbf{y}$

with

  • $\mathbf{\hat{y}}_{central}$ as the central part of the smoothed data of length $n - w + 1$
  • $\mathbf{H}_{central}$ as the central part of the projection matrix of $(n - w + 1)\times n$

In fact, the central part goes from index $\lceil\frac{w}{2}\rceil$ to $n-\lfloor\frac{w}{2}\rfloor$. Given a filter like the Savitzky-Golay filter with $w=5$ and polynomial degree $d\lt w$, which looks like

$\begin{bmatrix}c_2&c_1&c_0&c_1&c_2\end{bmatrix}$

and a signal $\mathbf{y}$ of length $n=10$ the central hat matrix $\mathbf{H}_{central}$ is given by

$\mathbf{H}_{central}=\begin{bmatrix} c_2&c_1&c_0&c_1&c_2&0&0&0&0&0\\ 0&c_2&c_1&c_0&c_1&c_2&0&0&0&0\\ 0&0&c_2&c_1&c_0&c_1&c_2&0&0&0\\ 0&0&0&c_2&c_1&c_0&c_1&c_2&0&0\\ 0&0&0&0&c_2&c_1&c_0&c_1&c_2&0\\ 0&0&0&0&0&c_2&c_1&c_0&c_1&c_2\\ \end{bmatrix}$

It can compute the smoothed data $\mathbf{\hat{y}}_{central}$ which have the indices from $\lceil\frac{w}{2}\rceil=\lceil\frac{5}{2}\rceil=3$ to $n-\lfloor\frac{w}{2}\rfloor=10-\lfloor\frac{5}{2}\rfloor=10-2=8$.

This matrix is very sparse, which in theory allows for efficient computation. To summarize graphically how this matrix achieves the shift and multiply-sum operation, let's look at the following step-by-step instructions:

Central Savitzky-Golay convolution process depicted step by step

1.2) Boundaries

Now comes the tricky part. There is not enough data on the left and right side of the signal to support a full filter window. Let's take a look at how the Python package SciPy handles this (only the default). For the function scipy.signal.savgol_filter the boundary handling is specified by the mode:

SciPy savgol_filter mode description

So, how can this be cast into the two projections

$\mathbf{\hat{y}}_{left}=\mathbf{H}_{left}\cdot\mathbf{y}_{fit;left}$

$\mathbf{\hat{y}}_{right}=\mathbf{H}_{right}\cdot\mathbf{y}_{fit;right}$

with

  • $\mathbf{y}_{fit;left}$ and $\mathbf{y}_{fit;right}$ as the left and right boundary basis data of length $w$ to which the polynomial is fitted (according to the mode description)
  • $\mathbf{H}_{left}$ and $\mathbf{H}_{right}$ as the left and right boundary hat matrices of shape $\lfloor\frac{w}{2}\rfloor\times w$
  • $\mathbf{\hat{y}}_{left}$ and $\mathbf{\hat{y}}_{right}$ as the smoothed left and right boundary data of length $\lfloor\frac{w}{2}\rfloor$

From the SciPy-description, a polynomial needs to be fitted to the boundary data and its evaluation can be written as

$\mathbf{\hat{y}}_{left}=\mathbf{V}_{eval;left}\cdot\mathbf{b}_{left}$

$\mathbf{\hat{y}}_{right}=\mathbf{V}_{eval;right}\cdot\mathbf{b}_{right}$

with

  • $\mathbf{V}_{eval;left}$ and $\mathbf{V}_{eval;right}$ as the Vandermonde matrices for the evaluation of the polynomials for the smoothed data
  • $\mathbf{b}_{left}$ and $\mathbf{b}_{right}$ as the polynomial coefficients of the fitted polynomials

For simplicity, only the left side is considered for now because the right side is easily derived from it.

Now, the polynomial evaluation can be related to $\mathbf{y}_{fit;left}$ by expressing the polynomial coefficients $\mathbf{b}_{left}$ as a function of $\mathbf{y}_{fit;left}$:

$\mathbf{b}_{left}=\mathbf{V^{+}}_{fit;left}\cdot\mathbf{y}_{fit:left}$

with

  • $\mathbf{V^{+}}_{fit;left}$ as the Moore-Penrose pseudo-inverse of the Vandermonde matrix $\mathbf{V}_{fit;left}$ for fitting the polynomial

Thanks to this, the polynomial evaluation can be written as

$\mathbf{\hat{y}}_{left}=\mathbf{V}_{eval;left}\cdot\mathbf{V^{+}}_{fit;left}\cdot\mathbf{y}_{left}=\mathbf{H}_{left}\cdot\mathbf{y}_{fit;left}$

so that the boundary hat matrix collapses to

$\mathbf{H}_{left}=\mathbf{V}_{eval;left}\cdot\mathbf{V^{+}}_{fit;left}$

Note that two different Vandermonde matrices are involved here, firstly the one for evaluation and secondly the one for fitting the polynomial.

Basically, the same applies to the right side, but it's not necessary to compute $\mathbf{H}_{right}$ explicitly because it can be derived from $\mathbf{H}_{left}$ by reversing its rows followed by reversing its columns.

1.3) Putting it all together

The three hat matrices $\mathbf{H}_{left}$, $\mathbf{H}_{central}$, and $\mathbf{H}_{right}$ can now be combined as

$\mathbf{H}=\begin{bmatrix} \begin{bmatrix}\mathbf{H}_{left}&0\end{bmatrix} \\\mathbf{H}_{central} \\\begin{bmatrix}0&\mathbf{H}_{right}\end{bmatrix} \end{bmatrix}$

The full smoothing process can be summed up in the animation that I uploaded to Wikimedia.

2) Cross-Validation and AIC

I'll first go into Cross-Validation because this is more natural to me.
The hat matrix has some nice properties. First of all, it gives the number of effective parameters of the smoothing process. There are different definitions out there like

$p_{eff}=trace\left(\mathbf{H}\right)$

or

$p_{eff}=trace\left(2\cdot\mathbf{H} - \mathbf{H}\cdot\mathbf{H}^{T}\right)$

where $trace(\mathbf{H})$ is the trace of the hat matrix, i.e., the sum of its main diagonal elements (taken from this Wikipedia article). According to Farhmeir, Kneib & Lang, Regression - Modelle, Methoden und Anwendungen (2009), the latter should be used when it comes to the estimation of the residual variance. For simplicity, let's take the first definition of $p_{eff}$.

With this, we basically got everything we need to perform a cross-validation. After smoothing, the residuals are given by

$\mathbf{r}=\mathbf{y}-\mathbf{\hat{y}}$

Making use of a further property of the hat matrix, the residual $r_{i}^{\left[-i\right]}$ that is observed at the $i$-th data point when the $i$-th data point is left out can be computed as

$r_{i}^{\left[-i\right]}=\frac{r_{i}}{1-h_{ii}}$

with $h_{ii}$ as the $i$-th main diagonal element of the hat matrix.

The cross-validation score can then be computed as the sum of the squared deleted residuals divided by the number of samples (according to Fahrmeir, et al., 2009):

$CV=\frac{1}{n}\sum_{i=1}^{n}\left(r_{i}^{\left[-i\right]}\right)^2$

To obtain the optimum smoothing parameters, the cross-validation score can be computed for different polynomial degrees and window lengths. The minimum of the cross-validation score will give the optimum parameters.

2.1) Simplifications for the CV-score

Since both the window length $w$ and the polynomial degree $d$ are integers, the optimization of the CV-score cannot be achieved by good continuous optimization algorithms.
Therefore, one has to rely on a grid search and then, the computation of the hat-matrix will turn out to be quite expensive.

Upon a closer look however, we may find that the only things required are the main diagonal elements of the hat matrix. Those can be obtained easily:

  • for the central hat matrix, the main diagonal elements are all identical, namely the central filter coefficient $c_{0}$, i.e, there is only a single calculation to make for getting $n - w + 1$ main diagonal elements at once.
  • for the left and right boundary hat matrices, the main diagonal elements can also be evaluated efficiently by taking $\sum _{j=1}^{d}{\left(\mathbf{V}_{eval;left}\circ\mathbf{V^{+}}_{fit;left}^{T}\right)_{ij}}$, i.e., the sum of each individual row of the element-wise product of the Vandermonde matrix for evaluation and the pseudo-inverse of the Vandermonde matrix for fitting.

If the boundaries are not considered in that detail and the assumption is made that they were still smoothed by a full filter window, everything simplifies a lot further because $c_{0}$ can be taken as the only main diagonal element of the hat matrix. This leaves

$p_{eff}=trace\left(\mathbf{H}\right)\approx n\cdot c_{0}$

$r_{i}^{\left[-i\right]}=\frac{r_{i}}{1-h_{ii}}=\frac{r_{i}}{1-c_{0}}$

$CV=\frac{1}{n}\sum_{i=1}^{n}\left(r_{i}^{\left[-i\right]}\right)^2\approx\frac{1}{n}\sum_{i=1}^{n}\left(\frac{r_{i}}{1-c_{0}}\right)^2=\frac{1}{n\cdot\left(1-c_{0}\right)^2}\sum_{i=1}^{n}r_{i}^2$

2.2) Link to the Akaike Information Criterion

The Akaike Information Criterion (AIC) for non-parametric regression is given by Fahrmeir, et al., 2009 as

$AIC=n\cdot\ln\left(\frac{1}{n}\sum_{i=1}^{n}r_{i}^2\right)+2\cdot k$

$k=\left(p_{eff}+1\right)$

With an adaption for small sample sizes, the AICc is given by

$AICc=AIC+\frac{2\cdot k\cdot\left(k+1\right)}{n-k-1}$

which I took from NIRPYResearch, The Akaike Information Criterion for model selection, 2021 (currently available here).

So, the same optimization can be performed with the AICc instead of the CV-score as well.

⚠️ A word of caution: The CV-Score and $AICc$ have to be computed for the filter of derivative order 0! If the first, second, or any other derivative are required from an optimum smooth, the CV-Score and $AICc$ have to first be optimized for derivative order 0 and the optimum parameters can then be used for the derivative computation. ⚠️

2.3) Results

Optimizing both the CV-score and the $AICc$ is carried out for the following random signal of size $n=501$. The width of the features in this signal is relatively uniform. Otherwise, the automated smoothing will be sub-optimum because the Savitzky-Golay filter is not spatially adaptive (in its basic version), which would have adverse effects on both criteria because the optimization would try to find a one-fits-all-solution which performs poor everywhere.

Test Signal

The grid search over multiple window-size-polynomial-degree-combinations gives the following optimization surfaces where the optimum window size for each polynomial degree is plotted as well:

Optimization surface for the Savitzky-Golay-Filter

The black area in the upper left corner was not computed because the window size $w$ was less than $d + 3$, i.e., the smoothing was either not possible or a perfect fit for $w = d + 1$.
For the upper right corner, the computation of the filter coefficients gets highly unstable. They are computed from polynomial matrices and these get badly conditioned in the SciPy implementation. Therefore, the coefficients are wrong and the smoothed series is terribly different from the original signal. However, it's nice to see that the criteria cover this well and exclude it naturally from the search space.
But what we actually see is that the difference between adjacent window sizes is quite small for higher polynomial degrees. Actually, there seems to be no such thing like a single optimum window size, but more of a continuum of good choices. This sort of justifies eye-balling to get the optimum smooth.

How do the smooths look like? Well, an eye-balled search could not be better for all of the polynomial degrees involved. The higher degrees look even smoother and also have the lower CV scores and $AICc$ when zoomed in.

Optimized smoothing with Savitzky-Golay

Zoom-In of the optimized smooths

3) Outlook

As mentioned in the beginning, the results apply to all Finite Impulse Response Filters (FIR), not only the Savitzky-Golay filter.
They also easily generalize for adaptive smoothing where each data point has its own polynomial degree and window length assigned if the main diagonal elements are not taken as a generic $h_{ii}=c_{0}$ for all, but as $h_{ii}=c_{0}\left(w_{i},d_{i}\right)$.
On top of that, access to the hat matrix allows for outlier handling, e.g., by making use of Cook's distance or the Peña sensitivity.

Edit
With the simplified formulas, the Cross-Validation approach - which was a Leave-One-Out-Cross-Valdiation - effectively becomes Generalized Cross-Validation (GCV).

Appendix - Implementation example

In Python, the simplified solutions could be implemented as follows:

### Imports ###

import numpy as np
from scipy.signal import savgol_coeffs, savgol_filter

### Function Definitions ###

def get_cv_score_simple(
    signal: np.ndarray,
    window_length: int,
    poly_degree: int,
) -> float:
    """
    Computes the simplified CV-score for a Savitzky-Golay filter applied to a signal.
    """

    # the central filter coefficient is computed
    central_coeff = savgol_coeffs(
        window_length=window_length,
        polyorder=poly_degree,
        use="dot",
    )[window_length // 2]

    # then, the smooth signal is computed
    smooth_signal = savgol_filter(
        x=signal,
        window_length=window_length,
        polyorder=poly_degree,
    )

    # finally, the residuals and cross-validation score are computed
    residuals = signal - smooth_signal

    return np.square(residuals).sum() / (((1.0 - central_coeff) ** 2) * signal.size)

def get_aicc_simple(
    signal: np.ndarray,
    window_length: int,
    poly_degree: int,
) -> float:
    """
    Computes the simplified AICc for a Savitzky-Golay filter applied to a signal.
    """

    # the central filter coefficient is computed
    central_coeff = savgol_coeffs(
        window_length=window_length,
        polyorder=poly_degree,
        use="dot",
    )[window_length // 2]

    # then, the smooth signal is computed
    smooth_signal = savgol_filter(
        x=signal,
        window_length=window_length,
        polyorder=poly_degree,
    )

    # from this, the residuals and Mean Squared Error are computed
    residuals = signal - smooth_signal
    mse = np.square(residuals).sum() / signal.size

    # the number of parameters is computed
    num_params = signal.size * central_coeff
    k = num_params + 1

    # then, the AICc is fully specified
    return (
        signal.size * np.log(mse)
        + 2.0 * k
        + ((2.0 * k * (k + 1)) / (signal.size - k - 1))
    )

$\endgroup$
6
  • 1
    $\begingroup$ This is a very solid answer! I thought it was an @whuber answer. Nicely done. So, please tell me about your favorite filters for robustly estimating derivatives, ones that are not variations on Savitzky-Golay. $\endgroup$ Commented May 30 at 14:51
  • 1
    $\begingroup$ @EngrStudent Thanks for the corrections, the upvote, and making it the accepted answer. I am planning to also comment on the proposal of the Hodrick-Prescott smoothing because it requires a lot of tweaking to get it fast. But once working, this is a good filter for doing the job when it comes to derivation. $\endgroup$
    – MothNik
    Commented May 30 at 17:12
  • 1
    $\begingroup$ Besides, I also found the Modified Sinc smoother mentioned to work like a charm for signals with uniform features. As for the Savitzky-Golay the hat-matrix is directly given by the coefficients, but this one can actually suppress high-frequency noise. Unlike the Hodrick-Prescott/Whittaker-Henderson it's not really spatially adaptive - which does not mean that the former is truly, they are just better in this respect. This can be alleviated when the curvature of the signal can be estimated because then the polynomial degree can be increased when curvature is high to preserve features. $\endgroup$
    – MothNik
    Commented May 30 at 17:16
  • 1
    $\begingroup$ I think this discussion easily gets too big for a comment section, though. In my opinion Hodrick-Prescott (smoothing splines in general) and the Modified Sinc smoother have their pros and cons, but the cons - be it runtime, hat matrix computation, or little flexibility - can be removed for either of them and then it becomes a matter of personal preference imo. $\endgroup$
    – MothNik
    Commented May 30 at 17:19
  • 1
    $\begingroup$ Another important point to consider: Do the datapoints have equal variance or not? The Savitzky-Golay filter can in theory be applied to unequal variance by rewriting the polynomial fit problem, but this is terribly slow. For the Modified Sinc, this gets way more invested and I would not know where to start implementing this because there is no such tracable origin like "Let's fit polynomials". In such cases, the Hodrick-Prescott/Whittaker-Henderson handles this all naturally and should be preferred. $\endgroup$
    – MothNik
    Commented May 30 at 17:46
6
$\begingroup$

Economists have some rules of thumb they use ... but they don't call it Savitzky-Golay (who wrote their paper in 1964), they call it Hodrick-Prescott (1997). [Actuaries call it Whittaker-Henderson graduation (the usual reference for Whittaker being his 1923 book, and for Henderson, papers in 1916 and 1924), so they do somewhat better at naming it for the right people. But it's basically just discrete smoothing splines, where derivatives of some order are replaced with differences of some order.]

However, those rules that are used for Hodrick-Prescott filtering are based on some assumptions and approximations and don't necessarily compare well with smoothing parameters chosen in other ways.

One thing you might do to choose a smoothing parameter is some form of crossvalidation (if you have an implementation that deals with missingness); if your primary aim is to forecast you could substitute mean square one-step ahead prediction error.

However, it's not hard to estimate degrees of freedom consumed by the fit; taking the approach in Ye (1998) [1] of summing the partial derivatives of the fit with respect to the observations $\sum_i \frac{\partial \hat{y}_i}{\partial {y}_i}$ to obtain model d.f. (a definition which is consistent with our usual notions when we apply that to cases where our usual notions work).

In this particular case, this calculation is easily written as the trace of a matrix. If you can write $\hat{Y} = HY$ - which is easy in this case - then model df is the trace of $H$.

(You could argue for it from more basic considerations but this is a fairly intuitive way to look at it to my mind, and generalizes to a very wide array of methods for getting a some fit.)

[1] Ye, J. (1998),
"On measuring and correcting the effects of data mining and model selection."
J Am Statist Assoc., 93, pp120–31.

$\endgroup$
2
  • $\begingroup$ Your help with nomenclature, with good references, and with the point on how to determine degrees of freedom are quite helpful. Thank you. $\endgroup$ Commented Apr 7, 2016 at 12:34
  • $\begingroup$ This leads to the Kalman filter. It is interesting that the Kalman filter and the AIC have a clean and direct bridge. How fun. $\endgroup$ Commented Apr 8, 2016 at 11:07
2
$\begingroup$

I don't think you can really avoid some subjectivity here. The whole concept of smoothing is based on the idea that there is some slowly-varying signal that will be more apparent or meaningful if you remove variation around it, and bandwidth is the parameter that determines where in frequency terms the line is drawn between signal and noise. That means that the ideal bandwidth is not simply a property of the data, but also of what you want to achieve.

If you can make explicit your objective for smoothing it may help to determine a sensible bandwith. For example if you just want "to make the chart look smooth" then eyeballing a selection and picking the one you like best does exactly what you want! Or perhaps you could pick a function to maximise such as at association between the smoothed series and some other series, but while this adds an objective element to the process you have still chosen both the series and measure of association.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.