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.