Fit a Piecewise Linear Model to Data with Unknown Knots and Number of Segments

Question

What is a robust way to fit piecewise linear but noisy data?

I'm measuring a signal, which consists of several almost linear segments. I'd like to atomatically fit several lines to the data to detect the transitions.

The dataset consists of a few thousand points, with 1-10 segments and I know the number of segments.

This is an example of what I'd like to do automatically.

Answer 1

I tried two approaches, naively (using only 3 segments). Surely there would be fancier methods out there.

RANSAC, supposed to be a robust fitting mechanism. It's easy to stop the algorithm after a number of segments. However it may be difficult to enforce continuity between segments--as seems required in your application-- at least with a simple implementation. As a proof of concept, I created an image from the data points so that I could use the RANSAC engine available in $ImageLines$, the line detection function of Mathematica.

Fit a piecewise linear model using a general purpose minimizer. It's easy to enforce segments continuity. Interestingly, testing for residuals and other properties may provide enough information to determine automatically the number of segments--I've not tried it though. That's how it looks in Mathematica:

Answer 2

I don't claim the following method is robust, but it might work for you. With thousands of points $x[n]$ and perhaps ten or so straight-line segments, proceed as follows.

• Process the points $x[n]$ to create a bit array $y[n]$ as follows. $$y[n] = \begin{cases}1, &\text{if} ~ |(x[n+1]-x[n]) - (x[n]-x[n-1])| < \epsilon,\\ 0, &\text{otherwise.}\end{cases}$$ Here $\epsilon$ is a small number chosen to suit your notion of how close to a straight line you want points $x[n-1],x[n], x[n+1]$ to hew to. The criterion will be recognized by the cognoscenti as demanding that the straight line through $(n-1, x[n-1])$ and $(n,x[n])$ has nearly the same slope as the straight line through $(n,x[n])$ and $(n+1,x[n+1])$.

• If $y[n]$ is an array of ten or so longish runs of $1$s separated by runs of $0$s with occasional stray $1$s here and there to mar the beauty, relax, you are on the right track. Else, if there are too few runs or too many runs of $1$s, repeat the previous step with a different $\epsilon$.

• Use linear least-mean-square-error curve-fitting to fit straight lines to the points identified by $y[n]$ as belonging to the same straight-line segment. You now have ten straight lines fitting points, say, Line A fits points $x[3]$ through $x[88]$; line B fits points $x[94]$ through $x[120]$, Line C fits points $x[129]$ through $\cdots$, and so on. Extend A rightwards and B leftwards to find out where they intersect; extend B rightwards and C leftwards to find out where they intersect, etc. Congratulations, you now have a continuous and piecewise linear model for your data.

Answer 3

(Years later) piecewise-linear functions are splines of degree 1, which most spline fitters can be told to do. scipy.interpolate.UnivariateSpline for example can be run with k=1 and a smoothing parameter s, which you'll have to play with -- see scipy-interpolation-with-univariate-splines .
In Matlab, see how-to-choose-knots .

Added: finding optimal knots is not easy, because there can be many local optima. Instead, you give UnivariateSpline a target s, sum of error^2, and let it determine the number of knots. After fitting, get_residual() will get the actual sum of error^2, and get_knots() the knots. A small change in s may change the knots a lot, especially in high noise -- ymmv.
The plot shows fits to a random piecewise-linear function + noise for various s.

For fitting piecewise constants, see Step detection. Can that be used for pw linear ? Don't know; starting off by differentiating noisy data will increase the noise, wrong.

Other testfunctions, and/or links to papers or code, would be welcome. A couple of links:
piecewise-linear-regression-with-knots-as-parameters
$\qquad$ Linear splines are very sensitive to where the knots are placed
knot-selection-for-cubic-regression-splines
$\qquad$ This is a tricky problem and most people just select the knots by trial and error.
$\qquad$ One approach which is growing in popularity is to use penalized regression splines instead.

---

Added March 2014: Dynamic programming is a general method for problems with nested subproblems like this:

optimal k lines
    = optimal k - 1 lines up to some x
    + cost of the last line x to the end
over x  (all x in theory, nearby x in practice)

Dynamic programming is very clever, but can it beat brute force + heuristics for this task ?
See the excellent course notes by Erik Demaine under MIT 6.006 Intro to algorithms
also google segmented linear regression
also John Henry syndrome.

---

Answer 4

I think a convex, robust optimization can be formulated to this problem based on a denoising model. It deals with the main challenge, estimating the location of the segment joins, the knots.

There are 2 assumptions to be made:

1. The model is piece wise linear.
2. The number of knots is sparse compared to the number samples.

The combination of the assumption means that the number of cases where the 2nd derivative of the estimated signal is not zero is sparse.
By using the ${L}_{1}$ norm to promote sparsity the problem can be formulated as:

$$ \arg \min_{\boldsymbol{x}} \underbrace{\frac{1}{2} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2}}_{\text{Denoising}} + \lambda \underbrace{\sum_{i = 2}^{n - 1} \left| {x}_{i - 1} - 2 {x}_{i} + {x}_{i + 1} \right|}_{\text{Sparse 2nd derivative}} = \frac{1}{2} {\left\| \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{D} \boldsymbol{x} \right\|}_{1} $$

This is very similar to the Total Variation (TV) Denoising model. The difference is in the $\boldsymbol{D}$ matrix. Where in the TV Denoising case it represents the 1st order forward finite differences operator and in this case it represents the central 2nd order finite differences operator.

The data:

With both the noise level and the knots unknown, this is the result of the model:

The problem is solved using ADMM which gives the same results as the DCP solver.

The full code is available on my StackExchange Signal Processing GitHub Repository (Look at the SignalProcessing\Q1227 folder).

Answer 5

Take the derivative and look for areas of nearly constant value. You would need to create the algorithm to search for those areas with ideally some level of +/- slope and that would give you the slope of the line for that section. You might want to perform some smoothing, like a sliding mean, prior to doing the sectional classification. The next step would be to get the y-intersection, which should be trivial at that point.

Answer 6

Using an l1 trend filter is another idea:

Paper

Online Example