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

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.

• I don't think this question can be answered reasonably unless you tell us how accurately you want to know the locations of the break-points, what your guesstimate is for the shortest length of a linear segment and how many samples there are in a typical transition region. If the horizontal axis labels in your figure are sample numbers, then, with two transitions in the span from $x[-5]$ to $x[0]$, the task is more difficult than if the straight-line segments were of longer duration (in samples). Jan 25, 2012 at 16:54
• @DilipSarwate I updated the Question with the requirements(btw the xaxis is the magnetic field in tesla) Jan 25, 2012 at 17:59
• You can try this toolbox if you are working with MATLAB curve fitting toolbox
– Rhei
Jan 18, 2015 at 17:21

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:

• Looks like a great answer. Thanks for contributing. Jan 26, 2012 at 13:51

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.

• Totally stole my answer! = ) Jan 26, 2012 at 19:21
• Interresting Idea but sadly due to noise on the signal I don't get good results. Jan 28, 2012 at 6:11
• That expression whose magnitute is being compared to epsilon is actually an approximation to the second derivative of the data. There are other ways to compute this using more than three point which don't respond to noise as much. Look up Savitzky-Golay. Nov 11, 2012 at 2:32

(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 John Henry syndrome.

• The problem, at least with scipy is the positioning of the knots. scipy uses equaly spaced knots. Sep 17, 2013 at 8:49
• @P3trus, yes for a start, but then they can move -- see the plot. Anyway it targets total error, not knots. Sep 18, 2013 at 12:29
• @P3trus Have you tried using multivariate regression splines method that automatically selects the breakpoints iteratively? cs.rtu.lv/jekabsons/regression.html Oct 13, 2013 at 20:21
• @Atul Ingle, afaik breakpoint / knot selection is the same problem, from whatever spline fitter. If you know of different algorithms for that from R / regression people, could you post a link please ? Oct 25, 2013 at 10:11
• Are looking for packages in R/Matlab that do adaptive regression splines? Here: cran.r-project.org/web/packages/earth/index.html cran.r-project.org/web/packages/mda/index.html and also ARESLab in Matlab that I already posted the link for. Oct 30, 2013 at 23:28

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).

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.

• derivative might be offul noisy. i don't think i would recommend that. Mar 10, 2014 at 3:47

Using an l1 trend filter is another idea:

Paper

Online Example

• Your answer is a little bit too short to be constructive! Please consider making an effort to expand it in a pedagogical way. Apr 30, 2014 at 7:49