Fast Algorithms for Segmented Regression

arXiv:1607.03990v1 [cs.LG] 14 Jul 2016

Jayadev Acharya
MIT
jayadev@csail.mit.edu

Ilias Diakonikolas
University of Southern California
ilias.d@ed.ac.uk

Jerry Li
MIT
jerryzli@csail.mit.edu

Ludwig Schmidt
MIT
ludwigs@mit.edu

July 15, 2016

Abstract
We study the fixed design segmented regression problem: Given noisy samples from a piecewise linear function f , we want to recover f up to a desired accuracy in mean-squared error.
Previous rigorous approaches for this problem rely on dynamic programming (DP) and, while
sample efficient, have running time quadratic in the sample size. As our main contribution, we
provide new sample near-linear time algorithms for the problem that – while not being minimax
optimal – achieve a significantly better sample-time tradeoff on large datasets compared to the
DP approach. Our experimental evaluation shows that, compared with the DP approach, our
algorithms provide a convergence rate that is only off by a factor of 2 to 4, while achieving
speedups of three orders of magnitude.

1

Introduction

We study the regression problem – a fundamental inference task with numerous applications that
has received tremendous attention in machine learning and statistics during the past fifty years (see,
e.g., [MT77] for a classical textbook). Roughly speaking, in a (fixed design) regression problem,
we are given a set of n observations (xi , yi ), where the yi ’s are the dependent variables and the xi ’s
are the independent variables, and our goal is to model the relationship between them. The typical
assumptions are that (i) there exists a simple function f that (approximately) models the underlying
relation, and (ii) the dependent observations are corrupted by random noise. More specifically, we
assume that there exists a family of functions F such that for some f ∈ F the following holds:
yi = f (xi ) + i , where the i ’s are i.i.d. random variables drawn from a “tame” distribution such as
a Gaussian (later, we also consider model misspecification).
Throughout this paper, we consider the classical notion of Mean Squared Error (MSE) to measure
the performance (risk) of an estimator. As expected, the minimax risk depends on the family F that
f comes from. The natural case that f is linear is fully understood: It is well-known that the leastsquares estimator is statistically efficient and runs in sample-linear time. The more general case that
f is non-linear, but satisfies some well-defined structural constraint has been extensively studied in a
variety of contexts (see, e.g., [GA73, Fed75, Fri91, BP98, YP13, KRS15, ASW13, Mey08, CGS15]).
In contrast to the linear case, this more general setting is not well-understood from an informationtheoretic and/or computational aspect.
In this paper, we focus on the case that the function f is promised to be piecewise linear with a
given number k of unknown pieces (segments). This is known as fixed design segmented regression,
and has received considerable attention in the statistics community [GA73, Fed75, BP98, YP13].
The special case of piecewise polynomial functions (splines) has been extensively used in the context
of inference, including density estimation and regression, see, e.g., [WW83, Fri91, Sto94, SHKT97,
Mey08].
Information-theoretic aspects of the segmented regression problem are well-understood: Roughly
speaking, the minimax risk is inversely proportional to the number of samples. In contrast, the
computational complexity of the problem is poorly understood: Prior to our work, known algorithms
for this problem with provable guarantees were quite slow. Our main contribution is a set of new
provably fast algorithms that outperform previous approaches both in theory and in practice. Our
algorithms run in time that is nearly-linear in the number of data points n and the number of
intervals k. Their computational efficiency is established both theoretically and experimentally. We
also emphasize that our algorithms are robust to model misspecification, i.e., they perform well even
if the function f is only approximately piecewise linear.
Note that if the segments of f were known a priori, the segmented regression problem could
be immediately reduced to k independent linear regression problems. Roughly speaking, in the
general case (where the location of the segment boundaries is unknown) one needs to “discover” the
right segments using information provided by the samples. To address this algorithmic problem,
previous works [BP98, YP13] relied on dynamic programming that, while being statistically efficient,
is computationally quite slow: its running time scales at least quadratically with the size n of the
data, hence it is rather impractical for large datasets.
Our main motivation comes from the availability of large datasets that has made computational
efficiency the main bottleneck in many cases. In the words of [Jor13]: “As data grows, it may be
beneficial to consider faster inferential algorithms, because the increasing statistical strength of the
data can compensate for the poor algorithmic quality.” Hence, it is sometimes advantageous to

sacrifice statistical efficiency in order to achieve faster running times because we can then achieve
the desired error guarantee faster (provided more samples). In our context, instead of using a slow
dynamic program, we employ a subtle iterative greedy approach that runs in sample-linear time.
Our iterative greedy approach builds on the work of [ADH+ 15, ADLS15], but the details of
our algorithms here and their analysis are substantially different. In particular, as we explain in
the body of the paper, the natural adaptation of their analysis to our setting fails to provide any
meaningful statistical guarantees. To obtain our results, we introduce novel algorithmic ideas and
carefully combine them with additional probabilistic arguments.

2

Preliminaries

In this paper, we study the problem of fixed design segmented regression. We are given samples
xi ∈ Rd for i ∈ [n] ( = {1, . . . , n}), and we consider the following classical regression model:
yi = f (xi ) + i .

(1)

Here, the i are i.i.d. sub-Gaussian noise variables with variance proxy σ 2 , mean E[i ] = 0, and
variance s2 = E[2i ] for all i.1 We will let  = (1 , . . . , n ) denote the vector of noise variables. We
also assume that f : Rd → R is a k-piecewise linear function. Formally, this means:
Definition 1. The function f : Rd → R is a k-piecewise linear function if there exists a partition
of the real line into k disjoint intervals I1 , . . . , Ik , k corresponding parameters θ1 , . . . , θk ∈ Rd , and
a fixed, known j such that for all x = (x1 , . . . , xd ) ∈ Rd we have that f (x) = hθi , xi if xj ∈ Ii . Let
Lk,j denote the space of k-piecewise linear functions with partition defined on coordinate j.
Moreover, we say f is flat on an interval I ⊆ R if I ⊆ Ii for some i = 1, . . . , k, otherwise, we
say that f has a jump on the interval I.
Later in the paper (see Section 6), we also discuss the setting where the ground truth f is not
piecewise linear itself but only well-approximated by a k-piecewise linear function. For simplicity
of exposition, we assume that the partition coordinate j is 1 in the rest of the paper.
Following this generative model, a regression algorithm receives the n pairs (xi , yi ) as input.
The goal of the algorithm is then to produce an estimate fb that is close to the true, unknown f
with high probability over the noise terms i and any randomness in the algorithm. We measure
the distance between our estimate fb and the unknown function f with the classical mean-squared
error:
n
1X
b
MSE(f ) =
(f (xi ) − fb(xi ))2 .
n
i=1

Rn×d

Throughout this paper, we let X ∈
be the data matrix, i.e., the matrix whose j-th row
T
is xj for every j, and we let r denote the rank of X. We also assume that no xi is the all-zeros
vector, since such points are trivially fit by any linear function.
The following notation will also be useful. For any function f : Rd → R, we let f ∈ Rn denote
the vector with components fi = f (xi ) for i ∈ [n]. For any interval I, we let X I denote the data
matrix consisting of all data points xi for i ∈ I, and for any vector v ∈ Rn , we let v I ∈ R|I| be the
vector of vi for i ∈ I.
1

We observe that s2 is guaranteed to be finite since the i are sub-Gaussian. The variance s2 is in general not
equal to the variance proxy σ 2 , however, it is well-known that s2 ≤ σ 2 .

3

We remark that this model also contains the problem of (fixed design) piecewise polynomial
regression as an important subcase. Indeed, imagine we are given x1 , . . . , xn ∈ R and y1 , . . . , yn ∈ R
so that there is some k-piece degree-d polynomial p so that for all i = 1, . . . , n we have yi = p(xi )+i ,
and our goal is to recover p in the same mean-squared sense as above. We can write this problem as a
k-piecewise linear fit in Rd+1 , by associating the data vector vi = (1, xi , x2i , . . . , xdi ) to each xi . If the
piecewise polynomial p has breakpoints at z1 , . . . zk−1 , the associated partition for the k-piecewise
linear function is I =P{(−∞, z1 ), [z1 , z2 ), . . . , (zk−2 , zk−1 ], (zk−1 , ∞)}. If the piecewise polynomial
is of the form p(x) = d`=0 aI` x` for any interval I ∈ I, then for any vector v = (v1 , v2 , . . . , vd+1 ) ∈
P
I
Rd+1 , the ground truth k-piecewise linear function is simply the linear function f (v) = d+1
`=1 a`−1 v`
for v2 ∈ I. Moreover, the data matrix associated with the data points v is a Vandermonde matrix.
For n ≥ d + 1 it is well-known that this associated Vandermonde matrix has rank exactly d + 1
(assuming xi 6= xj for any i 6= j).

2.1

Our Contributions

Our main contributions are new, fast algorithms for the aforementioned segmented regression problem. We now informally state our main results and refer to later sections for more precise theorems.
Theorem 2 (informal statement of Theorems 18 and 19). There is an algorithm GreedyMerge,
which, given X (of rank r), y, a target number of pieces k, and the variance of the noise s2 , runs in
time O(nd2 log n) and outputs an O(k)-piecewise linear function fb so that with probability at least
0.99, we have
!
r
kr
k
2
MSE(fb) ≤ O σ
+σ
log n .
n
n
That is, our algorithm runs in time which is nearly linear in n and still achieves a reasonable rate
of convergence. While this rate is asymptotically slower than the rate of the dynamic programming
estimator, our algorithm is significantly faster than the DP so that in order to achieve a comparable
MSE, our algorithm takes less total time given access to a sufficient number of samples. For more
details on this comparision and an experimental evaluation, see Sections 2.2 and 7.
At a high level, our algorithm proceeds as follows: it first forms a fine partition of [n] and then
iteratively merges pieces of this partitions until only O(k) pieces are left. In each iteration, the
algorithm reduces the number of pieces in the following manner: we group consecutive intervals
into pairs which we call “candidates”. For each candidate interval, the algorithm computes an error
term that is the error of a least squares fit combined with a regularizer depending on the variance
of the noise s2 . The algorithm then finds the O(k) candidates with the largest errors. We do not
merge these candidates, but do merge all other candidate intervals. We repeat this process until
only O(k) pieces are left.
A drawback of this algorithm is that we need to know the variance of the noise s2 , or at least
have a good estimate of it. In practice, we might be able to obtain such an estimate, but ideally
our algorithm would work without knowing s2 . By extending our greedy algorithm, we obtain the
following result:
Theorem 3 (informal statement of Theorems 21 and 22). There is an algorithm BucketGreedyMerge,
which, given X (of rank r), y, and a target number of pieces k, runs in time O(nd2 log n) and outputs

4

an O(k log n)-piecewise linear function fb so that with probability at least 0.99, we have
!
r
kr
log
n
k
MSE(fb) ≤ O σ 2
+σ
log n .
n
n
At a high level, there are two fundamental changes to the algorithm: first, instead of merging
with respect to the sum squared error of the least squares fit regularized by a term depending on
s2 , we instead merge with respect to the average error the least squares fit incurs within the current
interval. The second change is more substantial: instead of finding the top O(k) candidates with
largest error and merging the rest, we now split the candidates into log n buckets based on the
lengths of the candidate intervals. In this bucketing scheme, bucket α contains all candidates with
length between 2α and 2α+1 , for α = 0, . . . , log n − 1. Then we find the k + 1 candidates with largest
error within each bucket and merge the remaining candidate intervals. We continue this process
until we are left with O(k log n) buckets. Intuitively, this bucketing allows us to control the variance
of the noise without knowing s2 because all candidate intervals have roughly the same length.
A potential disadvantage of our algorithms above is that they produce O(k) or O(k log n) pieces,
respectively. In order to address this issue, we also provide a postprocessing algorithm that converts
the output of any of our algorithms and decreases the number of pieces down to 2k + 1 while
preserving the statistical guarantees above. The guarantee of this postprocessing algorithm is as
follows.
Theorem 4 (informal statement of Theorems 23 and 24). There is an algorithm Postprocessing
that takes as input the output of either GreedyMerge or BucketGreedyMerge together with
a target number of pieces k, runs in time O k 3 d2 log n , and outputs a (2k + 1)-piecewise linear
function fbp so that with probability at least 0.99, we have
!
r
kr
k
p
2
MSE(fb ) ≤ O σ
+σ
log n .
n
n
Qualitatively, an important aspect this algorithm is that its running time depends only logarithmically on n. In practice, we expect k to be much smaller than n, and hence the running time
of this postprocessing step will usually be dominated by the running time of the piecewise linear
fitting algorithm run before it.

2.2

Comparison to prior work

Dynamic programming. Previous work on segmented regression with statistical guarantees
[BP98, YP13] relies heavily on dynamic programming-style algorithms to find the k-piecewise linear
least-squares estimator. Somewhat surprisingly, we are not aware of any work which explicitly gives
the best possible running time and statistical guarantee for this algorithm. For completeness, we
prove the following result (Theorem 5), which we believe to be folklore. The techniques used in its
proof will also be useful for us later.
Theorem 5 (informal statement of Theorems 15 and 16). The exact dynamic program runs in time
O(n2 (d2 + k)) and outputs an k-piecewise linear estimator fb so that with probability at least 0.99 we
have


kr
MSE(fb) ≤ O σ 2
.
n
5

We now compare our guarantees with those of the DP. The main advantage of our approaches
is computational efficiency – our algorithm runs in linear time, while the running time of the DP
has a quadratic dependence on n. While our statistical rate of convergence is slower, we are able
to achieve the same MSE as the DP in asymptotically less time (and also in practice) if enough
samples are available.
For instance, suppose we wish to obtain a MSE of η with a k-piecewise linear function, and
suppose for simplicity that d = O(1) so that r = O(1) as well. Then Theorem 5 tells us that the
DP requires n = k/η samples and runs in time O(k 3 /η 2 ). On the other hand, Theorem 2 tells us
2 ) samples (ignoring log factors) and thus runs in time
e
that GreedyMerging requires n = O(k/η
2 ). For non-trivial values of k, this is a significant improvement in time complexity.
e
O(k/η
This gap also manifests itself strongly in our experiments (see Section 7). When given the same
number of samples, our MSE is a factor of 2-4 times worse than that achieved by the DP, but our
algorithm is about 1, 000 times faster already for 104 data points. When more samples are available
for our algorithm, it achieves the same MSE as the DP about 100 times faster.
Histogram Approximation Our results build upon the techniques of [ADH+ 15], who consider
the problem of histogram approximation. In this setting, one is given a function f : [n] → R, and
the goal is to find the best k-flat approximation to f in sum-squared error. This problem bears
nontrivial resemblance to the problem of piecewise constant regression, and indeed, the results in this
paper establish a connection between the two problems. The histogram approximation problem has
received significant attention in the database community (e.g. [JKM+ 98, GKS06, ADH+ 15]), and
it is natural to ask whether these results also imply algorithms for segmented regression. Indeed, it
is possible to convert the algorithms of [GKS06] and [ADH+ 15] to the segmented regression setting,
but the corresponding guarantees are too weak to obtain good statistical results. At a high level,
the algorithms
in these papers
be adapted to output an O(k)-piecewise linear function fb so
Pcan
Pn
n
2
b
that i=1 (yi − f (xi )) ≤ C · i=1 (yi − f LS (xi ))2 , where C > 1 is a fixed constant and f LS is the
best k-piecewise linear least-squares fit to the data. However, this guarantee by itself does not give
meaningful results for segmented regression.
As a toy example, consider the k = 1 case. Let x1 , . . . , xn ∈ Rn be arbitrary, let f (x) = 0 for
all x, and i ∼ N (0,
P 1), for i = 1, . . . , n. Then it is not too hard to see that the least squares fit
has squared error ni=1 (yi − f LS (xi ))2 = Θ(n) with high probability. Hence the function g which
is constantly C/2 for all x achieves the approximation guarantee
n
X

(yi − g(xi ))2 ≤ C

i=1

n
X
(yi − f LS (xi ))2
i=1

described above with non-negligible probability. However, this function clearly does not converge
to f as n → ∞, and indeed its MSE is always Θ(1).
To get around this difficulty, we must extend the algorithms presented in histogram approximation literature in order to adapt them to the segmented regression setting. In the process, we
introduce new algorithmic ideas and use more sophisticated proof techniques to obtain a meaningful
rate of convergence for our algorithms.

2.3

Agnostic guarantees

We also consider the agnostic model, also known as learning with model misspecification. So far,
we assumed that the ground truth is exactly a piecewise linear function. In reality, such a notion
6

is probably only an approximation. While the ground truth may be close to a piecewise linear
function, generally we do not believe that it exactly follows a piecewise linear function. In this case,
our goal should be to recover a piecewise linear function that is competitive with the best piecewise
linear approximation to the ground truth.
Formally, we consider the following problem. We assume the same generative model as in (1),
however, we no longer assume that f is a piecewise linear function. Instead, the function f can now
be arbitrary. We define
OPTk = min MSE(g)
g∈Lk

to be the error of the best fit k-piecewise linear function to f , and we let f ∗ be any k-piecewise
linear function that achieves this minimum. Then the goal of our algorithm is to achieve an MSE
as close to OPTk as possible. We remark that the qualitatively interesting regime is when OPTk is
small and comparable to the statistical error, and we encourage readers to think of OPTk as being
in that regime.
For clarity of exposition, we first prove non-agnostic guarantees. In Section 6, we then show how
to modify these proofs to obtain agnostic guarantees with roughly the same rate of convergence.

2.4

Mathematical Preliminaries

In this section, we state some preliminaries that our analysis builds on.
Tail bounds. We require the following tail bound on sub-exponential random variables. Recall
that for any sub-exponential random variable X, the sub-exponential norm of X, denoted kXkψ1 ,
is defined to be the smallest parameter K so that (E[|X|p ])1/p ≤ Kp for all p ≥ 1 (see [Ver10]).
Moreover, if Y is a sub-Gaussian random variable with variance σ 2 , then Y 2 − σ 2 is a centered
sub-exponential random with kY 2 − σ 2 kψ1 = O(σ 2 ).
Fact 6 (Bernstein-type inequality, c.f. Proposition 5.16 in [Ver10]). Let X1 , . . . , XN be centered
sub-exponential random variables, and let K = maxi kXi kψ1 . Then for all t > 0, we have
"
Pr

#
N
1 X
√
Xi ≥ t ≤ 2 exp −c min
N i=1

√ !!
t2 t N
,
.
K2 K

This yields the following straightforward corollary:
Corollary 7. Fix δ > 0 and let 1 , . . . , n be as in (1). Recall that s = E[2i ]. Then, with probability
1 − δ, we have that simultaneously for all intervals I ⊆ [n] the following inequality holds:
X

p
2i − s2 |I| ≤ O(σ 2 · log(n/δ)) |I| .

(2)

i∈I

Proof. Let δ 0 = O(δ/n2 ). For any interval I ⊆ [n], by Fact 6 applied to the random variables 2i − s2
for i ∈ I, we have that (2) holds with probability 1 − δ 0 . Thus, by a union bound over all O(n2 )
subintervals of [n], we conclude that Equation 2 holds with probability 1 − δ.

7

Linear regression. Our analysis builds on the classical results for fixed design linear regression.
In linear regression, the generative model is exactly of the form described in (1), except that f
is restricted to be a 1-piecewise linear function (as opposed to a k-piecewise linear function), i.e.,
f (x) = hθ ∗ , xi for some unknown θ ∗ .
The problem of linear regression is very well understood, and the asymptotically best estimator
is known to be the least-squares estimator.
Definition 8. Given x1 , . . . , xn P
and y1 , . . . , yn , the least squares estimator f LS is defined to be the
linear function which minimizes ni=1 (yi − f (xi ))2 over all linear functions f . For any interval I,
we let fILS denote the least squares estimator for the data points restricted to I, i.e. for the data
pairs {(xi , yi )}i∈I . We also let LeastSquares (X, y, I) denote an algorithm which solves linear
least squares for these data points, i.e., which outputs the coefficients of the linear least squares fit
for these points.
Following our previous definitions, we let f LS ∈ Rn denote the vector whose i-th coordinate is
f LS (xi ), and similarly for any I ⊆ [n] we let fILS ∈ R|I| denote the vector whose i-th coordinate for
i ∈ I is fILS (xi ).
The following prediction error rate is known for the least-squares estimator:
Fact 9. Let x1 , . . . , xn and y1 , . . . , yn be generated as in (1), where f is a linear function. Let
f LS (x) be the least squares estimator. Then,
 r


E MSE(f LS ) = O σ 2
.
n
Moreover, with probability 1 − δ, we have
MSE(f

LS



r + log 1/δ
)=O σ
n
2


.

Fact 9 can be proved with the following lemma, which we also use in our analysis. The lemma
bounds the correlation of a random vector with any fixed r-dimensional subspace:
Lemma 10 (c.f. proof of Theorem 2.2 in [Rig15]). Fix δ > 0. Let 1 , . . . , n be i.i.d. sub-Gaussian
random variables with variance proxy σ 2 . Let  = (1 , . . . , n ), and let S be a fixed r-dimensional
affine subspace of Rn . Then, with probability 1 − δ, we have
 p

|v T |
≤ O σ r + log(1/δ) .
v∈S\{0} kvk
sup

This lemma also yields the two following consequences. The first corollary bounds the correlation
between sub-Gaussian random noise and any linear function on any interval:
Corollary 11. Fix δ > 0 and v ∈ Rn . Let  = (1 , . . . , n ) be as in Lemma 10. Then with probability
at least 1 − δ, we have that for all intervals I, and for all non-zero linear functions f : Rd → R,
p
|hI , fI + vI i|
≤ O(σ r + log(n/δ)) .
kfI + vI k

8

Proof. Fix any interval I ⊆ [n]. Observe that any linear function on I can only take values in
the range {X I θ : θ ∈ Rd }, and hence the range of functions of the form f (xi ) + v is at most an
r-dimensional affine subspace. Thus, by Lemma 10, we know that for any linear function f ,
p
|hI , fI + vI i|
≤ O(σ r + log(n/δ)) .
kfI + vI k
with probability 1 − O(δ/n2 ). By union bounding over all O(n2 ) intervals, we achieve the desired
result.
Corollary 12. Fix δ > 0 and v ∈ Rn . Let  = (1 , . . . , n ) be as in Lemma 10. Then with probability
at least 1 − δ, we have

 r
|h, fI + vI i|
n
sup
.
≤ O σ kr + k log
δ
f ∈Lk kfI + vI k
Proof. Fix any partition of [n] into k intervals I, and let SI be the set of k-piecewise linear functions
which are flat on each I ∈ I. It is not hard to see that SI is a kr-dimensional linear subspace, and
hence the set of all k-piecewise linear functions which are flat on each I ∈ I when translated by v
is a kr-dimensional affine subspace. Hence Lemma 10 implies that
!
r
|h, fI + vI i|
1
sup
≤ O σ kr + log 0
δ
f ∈SI kfI + vI k

= nO(k)
with probability at least 1−δ 0 . A basic combinatorial argument shows that there are n+k−1
k−1

such different partitions I. Let δ 0 = δ/ n−k
k . Then the result follows by union bounding over all
the different possible partitions.

2.5

Runtime of linear least squares

The appeal of linear least squares is not only its statistical properties, but also that the estimator
can be computed efficiently. In our algorithm, we invoke linear least squares multiple times as a
black-box subrountine. Primarily, we use the following theorem:
Fact 13. Let A ∈ Rn×d be an arbitrary data matrix, and let y be the set of responses. Then there
is an algorithm LeastSquares(A, y) that runs in time O(nd2 ) and computes the least squares fit
to this data.
There has been work on faster, approximate algorithms that would suffice for our purposes in
theory. These algorithms offer a better dependence on the dimension d in exchange for slightly more
complicated approximation guarantees and somewhat more complicated algorithms (for instance,
see [CW13]). However, the classical algorithms for least squares with time complexity O(nd2 ) are
more commonly used in practice. For this reason and to simplify our exposition, we thus present our
results using the running time of the classical algorithms for least squares regression. Specifically,
we write LeastSquares(X, y, I) to denote an algorithm that computes a least squares fit on a
given interval I ⊆ [n] and assume that LeastSquares(X, y, I) runs in time O(|I| · d2 ) for all
I ⊆ [n].

9

3

Finding the least squares estimator via dynamic programming

In this section, we first present a dynamic programming approach (DP) to piecewise linear regression.
We do not believe these results to be novel, but to the best of our knowledge, these results appear
primarily as folklore in the literature. For completeness, we demonstrate the fastest DP we are
aware of, and we also prove its statistical guarantees. Not only will this serve as a good warm-up
for the later, more complex proofs, but we will also need a variant of this result in the later analyses.

3.1

The exact DP

We first describe the exact dynamic program. It will simply find the k-piecewise linear function
which minimizes the sum-squared error to the data. In other words, it outputs
arg min
f ∈Lk

n
X

(yi − f (xi ))2

=

arg minky − f k2 ,
f ∈Lk

i=1

which is simply the least-squares fit to the data amongst all k-piecewise linear functions. The
dynamic program computes the estimator f as follows. For i = 1, . . . , n and j = 1, . . . , k, let A(i, j)
denote the best error achievable by a j-piecewise linear function on the data pairs {(x` , y` )}i`=1 ,
that is, the best error achievable by j pieces for the first i data points. Then it is not hard to see
that

A(i, j) = min
err(X, y, {i0 + 1, . . . , i}) + A(i0 , j − 1) ,
0
i <i

where for any interval I, we let err(X, y, I) denote the sum-squared error to the data of the best
least squares fit to the data points {(x` , y` )}`∈I . That is, if fILS is the least squares fit to the data
in I, then err(X, y, I) = kyI − fILS k2 .
The algorithm then uses dynamic programming to fill out this n × k sized table of A values,
starting at i = 1 and j = 1. After having done so, the algorithm does one additional pass backwards
over the table to actually find the path through the table which achieves the best error. We can
optimize this by first computing the error quantities err(X, y, {i0 + 1, . . . , i}) for all i0 < i, and then
using a lookup table to find their values while actually executing the DP. Given such a lookup table,
the DP runs in time O(n2 k). Naively, the construction of this look-up table would take time which
is O(n3 d2 ) since there are O(n2 ) linear regression problems of size O(n), each of which takes O(nd2 )
time to solve.
However, we can speed this up. Consider a fixed interval I ⊂ [n]. Then the least squares fit on
this interval is of the form XIT XI θI = XIT yI . Assuming the matrix XIT XI is invertible, we have
−1 T
that θI = XIT XI
XI yI . Moreover, the error of the fit is given by
kXI θI − yI k2 = (XI θI − yI )T (XI θI − yI )
= θIT XIT XI θI − 2y T XI θI + yIT yI

−1 T
= yIT XI XIT XI
XI yI − 2yIT XI θI + yIT yI .
−1
The main point is the following: given XIT XI
and XIT yI , we can compute all the remaining
−1
quantities in time which is O(d2 ). To compute XIT XI
and XIT yI , we use the fact that our
regression instances are not arbitrary but very closely related. As a result, we can re-use computation

10

from previous regression instances that we have already solved in order to speed up later regression
instances.
Concretely, the calculations above indicate that we need to compute XIT yI for all intervals
I ⊆ [n]. But these are all just sub-vectors of the vector X T y, which we compute once in time
O(nd) and then use later on in all of the remaining computations. More non-trivially, we also
need to compute (XI XI )−1 for all I ⊆ [n]. However, observe that if we let I = {`, . . . , p − 1} for
T
` < p, then XI∪{p}
XI∪{p} = XIT XI + xp xTp , so that adding a single data point to the data matrix
corresponds to a rank-one update of the data matrix. Moreover, the effect of a rank-one update to
the inverse of the matrix is well-known:
Theorem 14 (Sherman-Morrison formula). Suppose M ∈ Rd×d is an invertible square matrix, and
let v ∈ Rd be such that 1 + v T M −1 v 6= 0. Then
M + vv T

−1

= M −1 −

M −1 vv T M −1
.
1 + v T M −1 v

Thus, given M −1 and v satisfying these conditions, we may compute M + vv T

−1

in time O(d2 ).
−1
Therefore, the dynamic program can do as follows: for each ` = 1, . . . , n, first compute XIT XI
for the interval I of length d starting at `, and then use the Sherman-Morrison update formula2 to
−1
compute XIT XI
for every interval I starting at ` in total time O(nd2 ). Thus the algorithm
−1
takes O(n2 d2 ) time total to compute all of the XIT XI
. As demonstrated earlier, this implies
the following:
Theorem 15. The exact dynamic program runs in time O(n2 (d2 + k)).

3.2

Error analysis for the exact DP

We now turn our attention to the learning rate of the exact DP. We show:
Theorem 16. Let δ > 0, and let f LS be the k-piecewise linear estimator returned by the exact DP.
Then, with probability 1 − δ, we we have that

n
LS
2 kr + k log δ
.
MSE(f ) ≤ O σ
n
Proof. We follow the proof technique for convergence of linear least squares as presented in [Rig15].
Recall f is the true k-piecewise linear function. Then, by the definition of the least squares fit, we
have that
ky − f LS k2 ≤ ky − f k2 = kk2 .
By expanding out ky − f LS k2 , we have that
ky − f LS k2 = kf +  − f LS k2
= kf LS − f k2 + 2 h, f − f LS i + kk2 .
2

(3)

In general, our matrices may not be invertible, or the update vector may not satisfy the condition in the ShermanMorrison formula, but in practice it seems these issues don’t come up and thus we do not worry about them here.
In general there are more complicated formulas for rank one updates for pseudo-inverses here but we will not cover
them for simplicity.

11

From these two calculations, we gather that
kf LS − f k2 ≤ 2 h, f − f LS i

 r
n
· kf LS − f k ,
≤ O σ kr + k log
δ
with probability 1−δ, where the last line follows from Corollary 12. A simple algebraic manipulation
from this last inequality yields the desired statement.
The same proof technique can also be easily adapted to yield the following slight extension of
Theorem 16, which can be proven via a union bound over all sets of k disjoint intervals. We omit
a proof here for conciseness.
Lemma 17. Fix δ > 0. Then with probability 1 − δ we have the following: for all disjoint sets of k
intervals I1 , . . . , Ik of [n] so that f is flat on each I` , the following inequality holds:
k
X

kfILS
− fI` k2 ≤ O(σ 2 k(r + log(n/δ))) .
`

`=1

4

A simple greedy merging algorithm

In this section, we give a novel greedy algorithm which runs much faster than the DP, but which
achieves a somewhat worse learning rate. However, we show both theoretically and experimentally
that the tradeoff between speed and statistical accuracy for this algorithm is markedly better than
it is for the exact DP.

4.1

The greedy merging algorithm

The overall structure of the algorithm is quite similar to [ADH+ 15], however, the merging criterion
is different, and as explained above, the guarantees proved in that paper are insufficient to give
non-trivial learning guarantees for regression.
Our algorithm here also requires an additional input s2 , which is defined to be the variance of
the i variables, i.e., s2 = E[2i ]. Requiring that we know s2 is a drawback, and in Section 5 we give
a slightly more complicated algorithm which does not require knowledge of s.
We give the formal pseudocode for the procedure in Algorithm 1. In the pseudocode we provide
two additional tuning parameters τ, γ. This is because in general our algorithm cannot provide a
k-histogram, but instead provides an O(k) histogram, which for most practical applications suffices.
The tuning parameters allow us to trade off running time for fewer pieces. In the typical use case
we will have τ, γ = Θ(1), in which case our algorithm will output an O(k)-piecewise linear function
in time O(nd2 log n) time.

4.2

Runtime of GreedyMerging

In this section we prove that our algorithm has the following, nearly-linear running time. The
analysis is similar to the analysis presented in [ADH+ 15].

Theorem 18. Let X and y be as in (1). Then GreedyMerging(τ, γ, s, X, y) outputs a (2 + τ2 )k + γ piecewise linear function and runs in time O(nd2 log(n/γ)).
12

Algorithm 1 Piecewise linear regression by greedy merging.
1: function GreedyMerging(τ, γ, s, X, y)
2:
. Initial partition of [n] into intervals of length 1.
3:
I 0 ← {{1}, {2}, . . . , {n}}
4:
. Iterative greedy merging (we start with j ← 0).
5:
while |I j | > (2 + τ2 )k + γ do
6:
Let sj be the current number of intervals.
7:
. Compute the least squares fit and its error for merging neighboring pairs of intervals.
s
8:
for u ∈ {1, 2, . . . , 2j } do
9:
θu ← LeastSquares(X, y, I2u−1 ∪ I2u )
10:
eu = kyI − XI θu k22 − s2 |I2u−1 ∪ I2u |
11:
end for
12:
Let L be the set of indices u with the (1 + τ1 )k largest errors eu ,
breaking ties arbitrarily.
13:
Let M be the set of the remaining indices.
14:
. Keep the
S intervals with large merging errors.
j+1
{I2u−1 , I2u }
15:
I
←
u∈L

16:
17:
18:
19:
20:
21:

. Merge the remaining intervals.
I j+1 ← I j+1 ∪ {I2u−1 ∪ I2u | u ∈ M }
j ←j+1
end while
return the least squares fit to the data on every interval in I j
end function

Before we prove this theorem, we compare this with the running time for the exact DP as given in
Theorem 15. Our main advantage is that our runtime scales linearly with n instead of quadratically.
This manifests itself as a substantially win theoretically in most reasonable regimes, and also as a
big win in practice—see our experiments for more details there.
Proof of Theorem 18. We first bound the time it takes to run any single iteration of the algorithm.
In any iteration j, we do a linear least squares regression problem on each interval I ∈ I j ; each such
problem takes O(|I|d2 ) time; hence solving them all takes O(nd2 ) time. Computing
 the eu given
the least squares fit takes no additional time asymptotically, and finding the 1 + τ1 largest errors
takes linear time [CLRS09]. Afterwards the remaining computations in this iteration can easily be
seen to be done in linear time. Hence each iteration takes O(nd2 ) time to complete.
We now bound the number of iterations of the algorithm. By the same analysis as that done in
the proof of Theorem 3.4 in [ADH+ 15] one can show that the algorithm terminates after at most
log(n/γ) iterations. Thus the whole algorithm runs in time O(nd2 log(n/γ)) time, as claimed.

13

4.3

Analysis of GreedyMerging

Theorem 19. Let δ > 0, and let fb be the estimator returned by GreedyMerging. Let m =
(2 + τ2 )k + γ be the number of pieces in fb. Then, with probability 1 − δ, we have that
!
√


n
τ
+
m(r
+
log(n/δ))
k
2
+σ √
.
MSE(fb) ≤ O σ
log
n
δ
n
Proof. We first condition on the event that Corollaries 7, 11 and 12, and Lemma 17 all hold with error
parameter O(δ), so that together they all hold with probability at least 1 − δ. Let I = {I1 , . . . , Im }
be the final partition of [n] that our algorithm produces. Recall f is the ground truth k-piecewise
linear function. We partition the intervals in I into two sets:
F = {I ∈ I : f is flat on I} ,
J = {I ∈ I : f has a jump on I} .
We first bound the error over the intervals in F. By Lemma 17, we have
X
kfI − fbI k2 ≤ O(σ 2 |F|(r + log(n/δ))) ,

(4)

I∈F

with probability at least 1 − O(δ).
We now turn our attention to the intervals in J and distinguish two further cases. We let J1
be the set of intervals in J which were never merged, and we let J2 be the remaining intervals. If
the interval I ∈ J1 was never merged, the interval contains one point, call it i. Because we may
assume that xi 6= 0, we know that for this one point, our estimator satisfies fb(xi ) = yi , since this is
clearly the least squares fit for a linear estimator on one nonzero point. Hence Corollary 7 implies
that the following inequality holds with probability at least 1 − O(δ):
X
X
kfI − fbI k2 =
kI k2
I∈J1

I∈J1


≤ σ2 

X



|I| + O log

I∈J1

≤ σ2



s
n X
δ


|I|

I∈J1

n √ 
m + O log
m .
δ


(5)

We now finally turn our attention to the intervals in J2 . Fix an interval I ∈ J2 . By definition,
the interval I was merged in some iteration of the algorithm. This implies that in that iteration,
there were (1 + 1/τ )k intervals M1 , . . . , M(1+1/τ )k so that for each interval M` , we have
LS 2
kyI − fbI k2 − s2 |I| ≤ kyM` − fM
k − s2 |M` | .
`

(6)

Since the true, underlying k-piecewise linear function f has at most k jumps, this means that there
are at least k/τ intervals of the M` on which f is flat. WLOG assume that these intervals are
M1 , . . . , Mk/τ .

14

Fix any ` = 1, . . . , k/τ . Expanding out the RHS of (6) using the definition of yi gives
LS
yM` − fM
`

2

LS
− s2 |M` | = fM` − fM
`
LS
= fM` − fM
`

2

LS
i + kM` k2 − s2 |M` |
+ 2hM` , fM` − fM
`
X
2
LS
(2i − s2 ) .
i
+
+ 2hM` , fM` − fM
`
i∈M`

Thus, we have that in aggregate,
k/τ
X

LS
yM` − fM
`

2

− s2 |M` | =

`=1

k/τ
X

LS
fM` − fM
`

2

`=1

+2

k/τ
k/τ
X
X
X
LS
i
+
(2i − s2 ) .
hM` , fM` − fM
`
`=1

`=1 i∈M`

(7)
We will upper bound each term on the RHS in turn. First, since the function f is flat on each M`
for ` = 1, . . . , k/τ , Lemma 17 implies
k/τ
X

fM` −

`=1

LS 2
fM
`




n
2k
≤O σ
,
r + log
τ
δ

(8)

with probability at least 1 − O(δ).
LS is a linear function on M of the form f LS (x) = xT β̂,
Moreover, note that the function fM
`
M`
`
LS
where β̂ ∈ Rd is the least-squares fit on M` . Because f is just a fixed vector, the vector fM` − fM
`
lives in the affine subspace of vectors of the form fM` + (XM` )η where η ∈ Rd is arbitrary. So
Corollary 12 and (8) imply that
v
u k/τ
k/τ
uX
X
LS
LS i · sup |hM` , Xηi|
hM` , fM` − fM` i ≤ t hM` , fM` − fM
`
kXηk
η
`=1
`=1


k
n
≤ O σ2
r + log
.
(9)
τ
δ
with probability 1 − O(δ).
By Corollary 7, we get that with probability 1 − O(δ),


k/τ

X
X
n √

2i − s2 |M` | ≤ O σ log
n.
δ
`=1

i∈M`

Putting it all together, we get that
k/τ 
X
i=1

yM` −

LS 2
fM
`





n √
k 2
n
σ r + log
+ O σ log
n
− s |M` | ≤ O
τ
δ
δ
2

(10)

with probability 1 − O(δ). Since the LHS of (6) is bounded by each individual summand above,
this implies that the LHS is also bounded by their average, which implies that
k/τ

τ X
2
2
LS 2
2
b
k
−
s
|M
|
kyI − fI k − s |I| ≤
kyM` − fM
`
`
k
i=1
√
 
 n 

n τ n
2
≤ O σ r + log
+ O σ log
.
δ
δ
k

15

(11)

We now similarly expand out the LHS of (6):
kyI − fbI k2 − s2 |I| = kfI + I − fbI k2 − s2 |I|
= kfI − fbI k2 + 2hI , fI − fbI i + kI k2 − s2 |I| .

(12)

P
From this we are interested in obtaining an upper bound on i∈I (f (xi ) − fb(xi ))2 , hence we seek
to lower bound the second and third terms of (12). The calculations here will very closely mirror
those done above.
By Corollary 11, we have that
 r
 n 
b
2hI , fI − fI i ≥ −O σ r + log
kfI − fbI k ,
δ
and by Corollary 7 we have

n p
kI k2 − s2 |I| ≥ −O σ log
|I| ,
δ
and so
 r
 n 

n p
2
2
2
b
b
kyI − fI k − s |I| ≥ kfI − fI k − O σ r + log
kfI − fbI k − O σ log
|I| .
δ
δ

(13)

Putting (11) and (13) together yields that with probability 1 − O(δ),
 
 n 
kfI − fbI k2 ≤ O σ 2 r + log
δ
 r
 √

 n 

n τ n p
+ O σ r + log
kfI − fbI k + O σ log
+ |I| .
δ
δ
k
Letting z 2 = kfbI − fI k2 , then this inequality is of the form z 2 ≤ bz + c where b, c ≥ 0. In this
specific case, we have that

 r
n
, and
b = O σ r + log
δ
 √

 n 
 

n τ n p
c = O σ 2 r + log
+ O σ log
+ |I|
δ
δ
k
We now prove the following lemma about the behavior of such quadratic inequalities:
Lemma 20. Suppose z 2 ≤ bz + c where b, c ≥ 0. Then z 2 ≤ O(b2 + c).
Proof. From the quadratic formula, the inequality implies that z ≤
straightforward to demonstrate the desired claim.

√
b+ b2 +4c
.
2

From this, it is

Thus, from the lemma, we have
 √

 
 n 

n τ n p
kfI − fbI k2 ≤ O σ 2 r + log
+ O σ log
+ |I| .
δ
δ
k

16

Hence the total error over all intervals in J2 can be bounded by:
!

Xp
√
n
+ O σ log
τ n+
|I|
kfI − fbI k ≤ O kσ r + log
δ
δ
I∈J
I∈J2


 n 

√ 
n  √
2
(14)
≤ O kσ r + log 0
+ O σ log
τ n + kn .
δ
δ
In the last line we use that the intervals I ∈ J2 are disjoint (and hence their cardinalities sum up to
at most n), and that there are at most k intervals in J2 because the function f is k-piecewise linear.
Finally, applying a union bound and summing (4), (5), and (14) yields the desired conclusion.
X

5

2



2



 n 



A variance-free merging algorithm

In this section, we give a variant of the above algorithm that does not require knowledge of the
noise variance s2 . The formal pseudo code is given in Algorithm 2.
Algorithm 2 Variance-free greedy merging with bucketing.
1: function BucketGreedyMerging(γ, X, y)
2:
. Initial histogram.
3:
Let I 0 ← {{1}, {2}, . . . , {n}} be the initial partition of [n] into intervals of length 1.
4:
. Iterative greedy merging (we start with j = 0).
5:
while |I j | > (2(k + 1) + γ) log n do
6:
Let sj be the current number of intervals.
7:
. Compute the least squares fit and its error for merging neighboring pairs of intervals.
s
8:
for u ∈ {1, 2, . . . , 2j } do
9:
Iu0 ← I2u−1 ∪ I2u
10:
θu ← LeastSquares(X, y, Iu0 )
11:
eu = |I10 | kyI − XI θu k22
u
12:
end for
13:
for α = 0, . . . , log(n) − 1 do
14:
Let Bα be the set of indices u so that 2α ≤ |Iu0 | ≤ 2α+1
15:
Let Lα be the set of indices u ∈ Bα with the k + 1 largest eu amongst u ∈ Bα ,
breaking ties arbitrarily.
16:
Let Mα be the set of the remaining indices u in Bα .
17:
. Keep theSintervals in each Bα with large merging errors.
18:
I j+1 ←
{I2u−1 , I2u }
u∈Lα

19:
20:
21:
22:
23:
24:
25:

. Merge the remaining intervals in each Bα .
I j+1 ← I j+1 ∪ {Iu0 | u ∈ Mα }
end for
j ←j+1
end while
return the least squares fit to the data on every interval in I j
end function

We now state the running time for BucketGreedyMerge. The running time analysis for
BucketGreedyMerge is almost identical to that of GreedyMerge; hence we omit its proof.
17

Theorem 21. Let X and y be as in (1). Then BucketGreedyMerging(γ, X, y) outputs a
(2(k + 1) + γ) log n-piecewise linear function and runs in time O(nd2 log(n/γ)).

5.1

Analysis of BucketGreedyMerge

This section is dedicated to the proof of the following theorem:
Theorem 22. Let fb be the m-piecewise linear function that is returned by BucketGreedyMerge,
where m = (2(k + 1) + γ) log n. Then, with probability 1 − δ, we have
!
r


n
k
2 m(r + log(n/δ))
b
+σ
MSE(f ) ≤ O σ
log
.
n
n
δ
Proof. As in the proof of Theorem 19, we let I be the final partition of [n] that our algorithm
produces. We also similarly condition on the event that Corollaries 7 and 12 and Lemma 17 all hold
with parameter O(δ). We again partition I into two sets F and J as before, and further subdivide
J into J1 and J2 . The error on F and J1 is the same as the error given in (4) and (5). The only
difference is the error on intervals in J2 .
Let I ∈ J2 be fixed. By definition, this means there was some iteration and a collection of some
k + 1 disjoint intervals M1 , . . . , Mk+1 so that |M` |/2 ≤ |I| ≤ 2|M` | and
1
yI − fbI
|I|

2

≤

1
LS
yM` − fM
`
|M` |

2

for all ` = 1, . . . , k + 1. Since f has at most k jumps, this means that there is at least one interval
M` so that f is flat on M` . WLOG assume that f is flat on M1 . By the same kinds of calculations
done in the proof of Theorem 19, we have that with probability 1 − δ,

p
n 
LS 2
2
yM1 − fM
) + σ 2 |M1 | + O(σ log(n/δ)) |M1 |
≤
O
σ
(r
+
log
1
δ
and

 p
p
2
2
yI − fbI ≥ fI − fbI − O σ r + log(n/δ) yI − fbI + σ 2 |I| − O(σ log(n/δ)) |I| ,
and putting these two equations together and rearranging, we have
fI − fbI

2


≤ O σ 2 (r + log(n/δ)) + σ 2 · O(log(n/δ))

!
p
|I|
p
+ |I|
|M` |

 p

+ O σ r + log(n/δ) yI − fbI
 p


p 

≤ O σ 2 (r + log(n/δ)) + O σ 2 log(n/δ) |I| + O σ r + log(n/δ) yI − fbI

,

where in the last inequality we used that |M` | ≥ |I|/2. By Lemma 20, this implies that

p 
2

fI − fbI ≤ O σ 2 (r + log(n/δ)) + O σ 2 log(n/δ) |I| .
Since there are at most k elements in J2 , and they are all disjoint, this yields that



n √ 
X
2
n 
fI − fbI ≤ O kσ 2 r + log
+ O σ log
kn
δ
δ
I∈J2

and putting this equation together with (4) and (5) yields the desired conclusion.
18

(15)

5.2

Postprocessing

One unsatisfying aspect of BucketGreedyMerge is that it outputs O(k log n) pieces. Not only
does this increase the size of the representation nontrivially when n is large, but it also increases the
error rate: it is the reason why the first term in the error rate given in Theorem 22 has an additional
log n factor as opposed the rate in Theorem 19. In this section, we give an efficient postprocessing
procedure for BucketGreedyMerge which, when run on the output of BucketGreedyMerge,
outputs an O(k) histogram with the same rate as before. In fact, the rate is slightly improved as
we are able to remove the log n factor mentioned above.
The postprocessing algorithm Postprocessing takes as input a partition I of [n] of size
O(k log n) and a target number of pieces k. It then performs the following steps: starting from
the O(k log n) partition I, run the dynamic program (DP) on intervals with breakpoints amongst
the breakpoints of I to find the 2k + 1 partition Ip (whose endpoints are also endpoints of I) that
minimizes the sum squared error to the data. The running time analysis is identical to that of the
DP with two exceptions: first, we only need to fill out a O(k log n) × (2k + 1) size table (as compared
to an n × k sized table). Second, we are no longer performing rank one updates because we instead
compute updates in large chunks, we cannot use the Sherman-Morrison formula to speed up the
computation of the least-squares fits. Hence, we obtain the following theorem:
Theorem 23. Given a partition I of [n] into O(k log n) pieces, Postprocessing(I, k) runs in
e 3 d2 ), and outputs a (2k+1)-piecewise linear function, where the O
e hides poly log(n) factors.
time O(k
We now show that the algorithm still provides the same (in fact, slightly better) statistical
guarantees as the original partition:
Theorem 24. Fix δ > 0. Let fbp be the estimator output by PostprocessedBucketGreedyMerge.
Then, with probability 1 − δ, we have
!

n
n rk
k
r
+
log
p
2
δ
MSE(fb ) ≤ O σ
+ σ log
.
n
δ
n
Proof. Let I and J be as in the proof of Theorem 22. Let Ip = {J1 , . . . , J2k+1 } be the intervals in
the partition that we return. We again condition on the event that Corollaries 7 and 12 and Lemma
17 all hold with parameter O(δ). The following will then all hold with probability 1 − δ.
Define the partition K to be the partition that contains every interval in J and exactly one
interval between any two non-consecutive intervals in J (i.e., K merges all flat intervals). Moreover,
let g be the (2k + 1)-piecewise linear function which is the least squares fit to the data on each
interval in I ∈ K. This is clearly a possible solution for the dynamic program, and therefore we
have
kfb − yk2 ≤ kg − yk2
X
X
=
kgI − yI k2 +
kgI − yI k2 .

(16)

I∈J

I∈K\J

We will expand and then upper bound the RHS of (16). First, by calculations similar to those
employed in the proof of Theorem 19, we have that


X
X
n 
kgI − yI k2 ≤ O kσ 2 r + log
+
kI k2 ,
δ
I∈K\J

I∈K\J

19

and from (15) and Corollary 12 we have
X
X
kgI − yI k2 =
kfI + I − gI k2
I∈J

I∈J

=

X

X

kfI − gI k2 + 2

I∈J



≤ O kσ

I∈J

2



hI , fI − gI i +

X

kI k2

I∈J


n √  X
n 
kn +
kI k2 ,
r + log
+ O σ log
δ
δ
I∈J

so all together now we have



n √ 
X
n 
kgI − yI k2 ≤ O kσ 2 r + log
+ O σ log
kn + kk2 .
δ
δ
I∈K

Moreover, by the same kinds of calculations, we can expand out the LHS of (16):
kfb − yk2 ≥ kfb − f k2 + h, fb − f i + kk2
r


√
n
2
b
k · σ r + log
kfb − f k + kk2
≥ kf − f k − O
δ
and hence, combining, cancelling, and moving terms around, we get that
r




n √ 
√
n 
2
b
kf − f k ≤ O
1 + kfb − f k + O σ log
k · σ r + log
kn .
δ
δ
By Lemma 20, this implies that

n √ 

n
+ σ log
kn ,
kfb − f k2 ≤ O kσ 2 r + log
δ
δ
with probability 1 − δ, as claimed.
We remark that Postprocessing can also be run on the output of GreedyMerging to
decrease the number of pieces from O(k) to 2k + 1 if so desired. The proof that it maintains similar
statistical guarantees is almost identical to the one presented above.

6

Obtaining agnostic guarantees

In this section, we demonstrate how to show agnostic guarantees for the algorithms in the previous
sections. Recall now f is arbitrary, and f ∗ is a k-piecewise linear function which obtains the
best approximation in mean-squared error to f amongst all k-piecewise linear functions. For all
i = 1, . . . , n, define ζi = f (xi ) − f ∗ (xi ) to be the error at data point i of the approximation, so that
for all i, we have
yi = f ∗ (xi ) + ζi + i .
(17)
By definition, we have that kζk2 = n · OPTk .
As a warm-up, we first show the following agnostic guarantee for the exact DP:

20

Theorem 25. Fix δ > 0. Let fLS be the k-piecewise linear function returned by the exact DP.
Then, with probability 1 − δ, we have
!
 r
n
OPTk
1
2 kr + k log δ
MSE(fLS ) ≤ O σ
+ σ log
+ OPTk
n
δ
n
In the regimes that are most interesting, i.e., when OPTk is small, the middle term does not
contribute significantly to the error.
Proof. The overall proof structure stays the same. From the definition of the least-squares fit, we
have
ky − f LS k2 ≤ ky − f ∗ k2
= k + ζk2
= kk2 + 2h, ζi2 + n · OPTk


1
2
(n · OPTk )1/2 + n · OPTk ,
≤ kk + O σ log
δ
with probability 1 − O(δ), where the last inequality follows from the r = 1 case of Lemma 10. The
LHS can be expanded out exactly as before in (3), and putting these two sides together we obtain
that


1
LS
2
LS
kf − f k ≤ 2h, f − f i + O σ log
(n · OPTk )1/2 + n · OPTk
δ
 r



n
1
LS
≤ O σ kr + k log
· kf − f k + O σ log
(n · OPTk )1/2 + n · OPTk
δ
δ
with probability at least 1 − O(δ) and so by Lemma 20 we obtain that
 

 
n
1
1/2
LS
2
2
+ σ log
(n · OPTk ) + n · OPTk .
kf − f k ≤ O σ kr + k log
δ
δ
Dividing both sides by n yields the desired conclusion.
Reviewing the proof, we observe the following general pattern: in all upper bounds we wish
to obtain (generally for formulas which occur on the RHS of the equations above) we obtain new
OPTk terms, whereas in the lower bounds (generally for formulas on the LHS of the equations
above) nothing changes. In fact, all folllowing proofs exhibit this pattern. We first establish some
useful concentration bounds. By the same union-bound technique used throughout this paper, one
can show the following. We omit the proof for conciseness.
Lemma 26. Fix δ > 0. Let i and ζi be as in (17), for i = 1, . . . , n. Then, with probability 1 − δ,
we have that for all collections of k disjoint intervals J1 , . . . , Jk , the following holds:
!1/2
k
k

X
X
n
2
hJ` , ζJ` i ≤ O kσ log
·
kζJ` k
δ
`=1
`=1

n
≤ O kσ log
· (n · OPTk )1/2
δ
21

Observe that with this lemma, we can easily adapt the proof of Theorem 25 to show the following
agnostic version of Lemma 17:
Lemma 27. Fix δ > 0. Then with probability 1 − δ, we have that for all disjoint sets of k intervals
J1 , . . . , Jk of [n] so that f ∗ is flat on each J` , the following inequality holds:
k
X
`=1

 
n

n
1/2
2
2
k
≤
O
σ
kr
+
k
log
kfJLS
−
f
+
kσ
log
(n
·
OPT
)
+
n
·
OPT
.
J`
k
k
`
δ
δ

With this lemma, we can now prove the following theorem for our greedy algorithm, which shows
that in the presence of model misspecification, our algorithm still performs at a good rate:
Theorem 28. Let δ > 0, and let fb be the estimator returned by GreedyMerging. Let m =
(2 + τ2 )k + γ be the number of pieces in fb. Then, with probability 1 − δ, we have that
!
√


n
 n  r OPT
k
m(r
+
log(n/δ))
τ
+
k
+σ √
log
MSE(fb) ≤ O σ 2
+ kσ log
+ τ · OPTk .
n
δ
δ
n
n
Proof. As before, we condition on the event that Corollaries 7, 11, and 12, and Lemmas 26 and 27
all hold with error parameter O(δ), so that together they all hold with probability at least 1 − δ.
We also let I, F, J1 , and J2 denote the same quantities as before, except we define the partition
with respect to the jumps of f ∗ instead of f , since the latter does not have well-defined jumps. For
instance, F is the set of intervals in I on which f ∗ has no jumps.
We again bound the error on these three sets separately. First, by Lemma 27 we have that

n
 
X
n
(18)
+ mσ log
(n · OPTk )1/2 + n · OPTk .
kfI − fbI k2 ≤ O σ 2 mr + m log
δ
δ
I∈F

Second, we bound the error on J1 . By modifying the calculations as those preceding (5) in the
same way as we did above in the proof of Theorem 25, we may show that




X
n √ 
n
kfI − fbI k2 ≤ σ 2 k + O log
m + O kσ log (n · OPTk )1/2 + n · OPTk .
(19)
δ
δ
I∈J1

Finally, we bound the error on J2 . Fix an I ∈ J2 , and let M1 , . . . , Mk/τ be as before. Recall that
this proof had two components: an upper bound for the M1 , . . . , M` and a lower bound for I. We
first compute the upper bound. By following the calculations for (11) and using Lemma 27, we have
k/τ
X
`=1

kyM` −

LS 2
fM
k
`

2



 n 
n √
k 2
σ r + log
+ σ log
n
τ
δ
δ

n
k
1/2
+ σ log
(n · OPTk ) + n · OPTk
τ
δ


− s |M` | ≤ O

The lower bound is in fact unchanged: we may use (13) as stated. Putting these two terms together
and simplifying as we did previously, we obtain that
 
 √

 n 
n τ n p
2
2
b
kfI − fI k ≤ O σ r + log
+ σ log
+ |I|
δ
δ
k
n

τ
+σ log
(n · OPTk )1/2 + n · OPTk .
δ
k
As before, summing up all the bounds we have achieved proves the desired claim.
22

Piecewise constant

MSE

10−1
10−2
103
n

104

4
2
103
n

104

103

100
10−2

101

10−4
102

Piecewise linear

102

103
n

104

102

103
n

104

103
n

104

10−2
103
n

104

8
6
4
2
102

103
n

Merging k

104

102
Speed-up

MSE

10−1

Running time (s)

10

100

102

6

102

Relative MSE ratio

102

8

Speed-up

Relative MSE ratio

100

Running time (s)

10

100
10−2

102

10−4
102

Merging 2k

103

103
n

Merging 4k

104

102

Exact DP

Figure 1: Experiments with synthetic data: results for piecewise constant models with k = 10
segments (top row) and piecewise linear models with k = 5 segments (bottom row, dimension
d = 10). Compared to the exact dynamic program, the MSE achieved by the merging algorithm
is worse but stays within a factor of 2 to 4 for a sufficient number of output segments. The
merging algorithm is significantly faster and achieves a speed-up of about 103 × compared to the
dynamic program for n = 104 . This leads to a significantly better trade-off between statistical and
computational performance (see also Figure 2).
Through virtually identical methods we also obtain the same upper bound for BucketGreedyMerge
and PostprocessedBucketGreedyMerge. We state the results but omit their proofs for this
reason.
Theorem 29. Let fb be the m-piecewise linear function that is returned by BucketGreedyMerge,
where m = (2 + τ2 )k log n + γ. Then, with probability 1 − δ, we have
!
r


n
 n  r OPT
m(r
+
log(n/δ))
k
k
2
MSE(fb) ≤ O σ
+σ
log
+ kσ log
+ τ · OPTk .
n
n
δ
δ
n
Theorem 30. Fix δ > 0. Let fbp be the estimator output by PostprocessedBucketGreedyMerge.
Then, with probability 1 − δ, we have
!

n
n rk
 n  r OPT
k
2 k r + log δ
b
MSE(fp ) ≤ O σ
+ σ log
+ kσ log
+ τ · OPTk .
n
δ
n
δ
n

23

Piecewise constant

Piecewise linear

100
MSE

MSE

100

10−2

10−4
10−4

10−2

10−2
Time (s)

10−4
10−4

100

10−1
Time (s)

102

Figure 2: Computational vs. statistical efficiency in the synthetic data experiments. The solid lines
correspond to the data in Figure 1, the dashed lines show the results from additional runs of the
merging algorithms for larger values of n. The merging algorithm achieves the same MSE as the
dynamic program about 100× faster if a sufficient number of samples is available.
Dow Jones data

Piecewise linear approximation

Figure 3: Results of fitting a 5-piecewise linear function (d = 2) to a Dow Jones time series. The
merging algorithm produces a fit that is comparable to the dynamic program and is about 200×
faster (0.013 vs. 3.2 seconds).

7

Experiments

In addition to our theoretical analysis above, we also study the empirical performance of our new
estimator for segmented regression on both real and synthetic data. As baseline, we compare our
estimator (GreedyMerging) to the dynamic programming approach. Since our algorithm combines both combinatorial and linear-algebraic operations, we use the Julia programming language3
(version 0.4.2) for our experiments because Julia achieves performance close to C on both types of
operations. All experiments were conducted on a laptop computer with a 2.8 GHz Intel Core i7
CPU and 16 GB of RAM.
Synthetic data. Experiments with synthetic data allow us to study the statistical and computational performance of our estimator as a function of the problem size n. Our theoretical
p bounds
indicate that the worst-case performance of the merging algorithm scales as O( kd
+
k/n log n)
n
kd
for constant error variance. Compared to the O( n ) rate of the dynamic program, this indicates
that the relative performance of our algorithm can depend on the number of features d. Hence we
use two types of synthetic data: a piecewise-constant function f (effectively d = 1) and a piecewise
linear function f with d = 10 features.
We generate the piecewise constant function f by randomly choosing k = 10 integers from the
3

http://julialang.org/

24

set {1, . . . , 10} as function value in each of the k segments.4 Then we draw n/k samples from each
segment by adding an i.i.d. Gaussian noise term with variance 1 to each sample.
For the piecewise linear case, we generate a n × d data matrix X with i.i.d. Gaussian entries
(d = 10). In each segment I, we choose the parameter values βI independently and uniformly
at random from the interval [−1, 1]. So the true function values in this segment are given by
fI = XI βI . As before, we then add an i.i.d. Gaussian noise term with variance 1 to each function
value.
Figure 1 shows the results of the merging algorithm and the exact dynamic program for sample
size n ranging from 102 to 104 . Since the merging algorithm can produce a variable number of output
segments, we run the merging algorithm with three different parameter settings corresponding to k,
2k, and 4k output segments, respectively. As predicted by our theory, the plots show that the exact
dynamic program has a better statistical performance. However, the MSE of the merging algorithm
with 2k pieces is only worse by a factor of 2 to 4, and this ratio empirically increases only slowly
with n (if at all). The experiments also show that forcing the merging algorithm to return at most
k pieces can lead to a significantly worse MSE.
In terms of computational performance, the merging algorithm has a significantly faster running
time, with speed-ups of more than 1, 000× for n = 104 samples. As can be seen in Figure 2, this
combination of statistical and computational performance leads to a significantly improved trade-off
between the two quantities. When we have a sufficient number of samples, the merging algorithm
achieves a given MSE roughly 100× faster than the dynamic program.
Real data. We also investigate whether the merging algorithm can empirically be used to find
linear trends in a real dataset. We use a time series of the Dow Jones index as input, and fit a
piecewise linear function (d = 2) with 5 segments using both the dynamic program and our merging
algorithm with k = 5 output pieces. As can be seen from Figure 3, the dynamic program produces
a slightly better fit for the rightmost part of the curve, but the merging algorithm identifies roughly
the same five main segments. As before, the merging algorithm is significantly faster and achieves
a 200× speed-up compared to the dynamic program (0.013 vs 3.2 seconds).

Acknowledgements
Part of this research was conducted while Ilias Diakonikolas was at the University of Edinburgh,
Jerry Li was an intern at Microsoft Research Cambridge (UK), and Ludwig Schmidt was visiting
the EECS department at UC Berkeley.
Jayadev Acharya was supported by a grant from the MIT-Shell Energy Initiative. Ilias Diakonikolas was supported in part by EPSRC grant EP/L021749/1, a Marie Curie Career Integration Grant, and a SICSA grant. Jerry Li was supported by NSF grant CCF-1217921, DOE grant
DE-SC0008923, NSF CAREER Award CCF-145326, and a NSF Graduate Research Fellowship.
Ludwig Schmidt was supported by grants from the MIT-Shell Energy Initiative, MADALGO, and
the Simons Foundation.
4

We also repeated the experiment for other values of k. Since the results are not qualitatively different, we only
report the k = 10 case here.

25

References
[ADH+ 15] J. Acharya, I. Diakonikolas, C. Hegde, J. Z. Li, and L. Schmidt. Fast and near-optimal
algorithms for approximating distributions by histograms. In PODS, pages 249–263,
2015.
[ADLS15] J. Acharya, I. Diakonikolas, J. Zheng Li, and L Schmidt. Sample-optimal density estimation in nearly-linear time. CoRR, abs/1506.00671, 2015.
[ASW13]

H. Avron, V. Sindhwani, and D. Woodruff. Sketching structured matrices for faster
nonlinear regression. In NIPS, pages 2994–3002. 2013.

[BP98]

J. Bai and P. Perron. Estimating and testing linear models with multiple structural
changes. Econometrica, 66(1):47–78, 1998.

[CGS15]

S. Chatterjee, A. Guntuboyina, and B. Sen. On risk bounds in isotonic and other shape
restricted regression problems. Annals of Statistics, 43(4):1774–1800, 08 2015.

[CLRS09] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein. Introduction to Algorithms.
3rd edition, 2009.
[CW13]

K. L. Clarkson and D. P. Woodruff. Low rank approximation and regression in input
sparsity time. In STOC, 2013.

[Fed75]

P. I. Feder. On asymptotic distribution theory in segmented regression problems– identified case. Annals of Statistics, 3(1):49–83, 01 1975.

[Fri91]

J. H. Friedman. Multivariate adaptive regression splines. Annals of Statistics, 19(1):1–
67, 03 1991.

[GA73]

A. R. Gallant and Fuller W. A. Fitting segmented polynomial regression models whose
join points have to be estimated. Journal of the American Statistical Association,
68(341):144–147, 1973.

[GKS06]

S. Guha, N. Koudas, and K. Shim. Approximation and streaming algorithms for histogram construction problems. ACM Trans. Database Syst., 2006.

[JKM+ 98] H. V. Jagadish, Nick Koudas, S. Muthukrishnan, Viswanath Poosala, Kenneth C. Sevcik,
and Torsten Suel. Optimal histograms with quality guarantees. In VLDB ’98, 1998.
[Jor13]

M. I. Jordan. On statistics, computation and scalability. Bernoulli, 19(4):1378–1390, 09
2013.

[KRS15]

R. Kyng, A. Rao, and S. Sachdeva. Fast, provable algorithms for isotonic regression in
all lp -norms. In NIPS, pages 2701–2709, 2015.

[Mey08]

M. C. Meyer. Inference using shape-restricted regression splines. Annals of Applied
Statistics, 2(3):1013–1033, 09 2008.

[MT77]

F. Mosteller and J. W. Tukey. Data analysis and regression: a second course in statistics.
Addison-Wesley, Reading (Mass.), Menlo Park (Calif.), London, 1977.
26

[Rig15]

P. Rigollet. High dimensional statistics. 2015.

[SHKT97] C. J. Stone, M. H. Hansen, C. Kooperberg, and Y. K. Truong. Polynomial splines and
their tensor products in extended linear modeling: 1994 wald memorial lecture. Annals
of Statistics, 25(4):1371–1470, 1997.
[Sto94]

C. J. Stone. The use of polynomial splines and their tensor products in multivariate
function estimation. Annals of Statistics, 22(1):pp. 118–171, 1994.

[Ver10]

Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices.
Chapter 5 of: Compressed Sensing, Theory and Applications. Edited by Y. Eldar and
G. Kutyniok. Cambridge University Press, 2012, 2010.

[WW83]

E. J. Wegman and I. W. Wright. Splines in statistics. Journal of the American Statistical
Association, 78(382):pp. 351–365, 1983.

[YP13]

Y. Yamamoto and P. Perron. Estimating and testing multiple structural changes in
linear models using band spectral regressions. Econometrics Journal, 16(3):400–429,
2013.

27

