1
$\begingroup$

I have a one-dimensional set of points, i.e. $(n,y_n), 1\leq n \leq N$. I want to fit them with a linear combination of $k$ rectangular functions in a least-squared-error sense. Each rectangle is parametrised by a left edge, right edge, and height. It is similar to the Lebesgue sum but I am looking at horizontal strips of flexible width and height to achieve minimum error.

It seems to be a really straightforward and very linear problem, but I don't see this facility in MATLAB curve fitting toolbox, and the internet is not helping at all.

I am sure (I think) I can formulate the equations and write the code, it might take me a couple of hours. But should I do it? Is there some reason it's not widely used? Is it actually ill-defined (in the sense there is no optimum)?

Improve this question
$\endgroup$
14
  • $\begingroup$ What do you mean by a 1D set of points? A dotted or discrete curve? $\endgroup$
    – Laurent Duval
    Mar 20, 2016 at 12:57
  • $\begingroup$ Discrete set of points. X-axis is just natural numbers 1 to N. $\endgroup$
    – Milind R
    Mar 20, 2016 at 13:39
  • $\begingroup$ So you have a set of coordinates $(n, Y_n)$, $1\le n \le N$. What do you mean by "fit" ? Have each point in at least one rectangle? Are you somehow paving your discrete curve with rectangles? $\endgroup$
    – Laurent Duval
    Mar 20, 2016 at 13:57
  • $\begingroup$ @LaurentDuval Thanks for the help, question edited for clarity. $\endgroup$
    – Milind R
    Mar 20, 2016 at 14:16
  • 2
    $\begingroup$ Well, for the time being I've given a solution to a problem that is at least related to what you want to do -- please have a look. $\endgroup$
    – j_random_hacker
    Apr 5, 2016 at 9:59

1 Answer 1

Reset to default
1
$\begingroup$

Assuming that

  • the rectangles do not overlap at all on the $x$ axis, and that
  • the sum of the rectangles should never exceed the value of the original function at any point,

the following simple dynamic programming algorithm will calculate the optimal set of at most $k$ rectangles w.r.t. least squared error in $O(n^2k)$ time and $O(nk)$ space in the worst case. Although this might seem prohibitive for large $n$, there are some refinements that mean that many inputs will take much less time -- and if necessary, some shortcuts that improve the runtime at the cost of sacrificing optimality.

Define $f(i, j)$ to be the least error that can be achieved on the subsequence $1, \dots, i$ ($i \le n$) using at most $j \le k$ rectangles. To avoid confusion about which data points are supported by which rectangles, assume all rectangles begin and end at positions having fractional part 0.5. Also define $g(i, j)$ to be the least error that can be achieved for the subsequence $i, ..., j$ using a single rectangle that ends at $x$ co-ord $j + 0.5$ and begins at any position $\ge i + 0.5$ and $\le j - 0.5$ and having fractional part 0.5. Finally, let $adjErr(i, j, h)$ be the sum of squared residuals over the range $i, ..., j$ under the assumption that a single rectangle of height $h$ spans this region. Then:

$f(i, j) = \min_{0 \le m < j}{(f(m, j-1) + g(m+1, i))}$

$g(i, j) = \min_{i \le m < j}{(adjErr(i, m, 0) + adjErr(m+1, j, min_{m<r\le j}x_r))}$

$adjErr(i, j, h) = \sum_{m=i}^j(x_m - h)^2$

The overall lowest possible error is given by $f(n, k)$. An actual solution (i.e. a set of rectangles) having this best-possible error level can be found by tracing back through the DP matrix, looking for the $m$ value that allowed $f(i, j)$ to take its minimum value.

If this is useful, comment and I'll come back to flesh out the details a bit more.

Improve this answer
$\endgroup$
6
  • $\begingroup$ That looks great! I came as far as $adjErr()$, but then only thought of brute-forcing it instead of DP, which is what I'm trying to understand now. Will update. $\endgroup$
    – Milind R
    Apr 6, 2016 at 6:09
  • $\begingroup$ Glad to hear it :) Don't have much time right now, but I'm now wondering whether computing $adjErr()$ causes the time complexity to exceed $O(n^2k)$... Anyway, if you're prepared to (or would prefer to!) use rectangles that can be higher than some input points, then the optimal height $h$ to use with $adjErr()$ is just the mean of all samples in the range $i, ..., j$, and $adjErr(i, j)$ is then just the variance times $j - i + 1$, which can be computed efficiently online. $\endgroup$
    – j_random_hacker
    Apr 6, 2016 at 10:11
  • $\begingroup$ I had assumed exactly that (rectangle height can be greater than input values), and got exactly what you mentioned :). I have a feeling though that keeping heights strictly lesser than input values is more realistic. But with more than one rectangle, the local optimum is no longer an optimum. Have to verify though, after understanding your DP. $\endgroup$
    – Milind R
    Apr 6, 2016 at 10:16
  • $\begingroup$ I'm sure the formulation in my answer is correct for the always-stay-at-or-below-the-data-points version of the problem (unless I made an off-by-one error with the indices, as I often do...), but it's probably slower to compute. We don't have any "local" optima to worry about: every possible solution to a given subproblem $(i, j)$ either has a rightmost rectangle ending at $x = i + 0.5$, or its rightmost rectangle ends at some $i' + 0.5 < i + 0.5$ (imagine an empty solution has a rectange ending at -0.5). ... $\endgroup$
    – j_random_hacker
    Apr 6, 2016 at 18:14
  • $\begingroup$ ... In the latter case, the best possible solution is given by $f(i', j)$. In the former case, we try all possible rectangles ending at $i + 0.5$ and combine each of them with the best overall solution that ends anywhere to the left of that particular rectangle and uses $\le j-1$ rectangles. We don't need to try combining a rightmost rectangle with all possible such leftward sub-solutions because the two parts (rightmost rectangle, and leftward sub-solution) don't overlap, so their costs are simply summed. $\endgroup$
    – j_random_hacker
    Apr 6, 2016 at 18:39

Your Answer

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

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