Hello mathoverflow community !


I have a simple question that seems to have a non trivial answer.


Given a discrete one dimensional signal $w(x)$ defined in a finite range, and the boxcar (rectangular) function $r(x)$


$$r(x)=\begin{cases}
1 & \mbox{if }0\leq x \leq 1; \\
0 & \mbox{elsewhere}
\end{cases}$$


I would like to find the coefficients $a\_i,\ b\_i,\ c\_i $ of the sum 


$$w' = \sum\_{i=0}^{N}\ { a\_i \cdot r\left(\frac{x}{b\_i} - c\_i\right)}$$


("sum of $N$ rectangles in any range and of any height")
such as $\sum\_i\ \left| w\_i - w\_i'\right|$ is minimized (for a given $N$).


This problem seems related to:


1. Discrete wavelet transform
2. $l\_1$ regularized solution of an overdetermined linear system
3. Maximum subarray problem


However, to my understanding it does not fit any of these cases:


1. $r(x)$ is not a wavelet basis,
2. the problem cannot be solved (practically) as a linear system because the (finite) set of $a\_i,\ b\_i,\ c\_i $ values is too large to compute explicitly (length of $w$ in the order of $10^4$),
3. Since $a\_i$ is undefined, it does not fit as a maximum subarray problem.


Right now I have an approximate solution (iteratively solving the problem via maximum subarray formulation by brute force exploring a subset of possible $a\_i$ values), however the idea of "decomposing a signal as a sum of rectangles" seems general enough to think that someone has already addressed it in the past.



> 
> Do any of you have a suggestion on how to tackle this problem ?
> 
> 
> Has it already been solved in the past, by a method I am not aware of ?
> 
> 
> 


Thank you very much for your answers.



Decomposing a discrete signal into a sum of rectangle functions
Dear Raskolnikov, Brian and Sebastian,
thanks you very much for sharing your ideas. 


I cross posted my question in two forums, so you can check the complete discussion by looking [here](https://math.stackexchange.com/questions/35388) and [here](https://mathoverflow.net/questions/63144).


From the discussion I got three main elements:


1. Haar wavelet provide a solution to my problem, although a very suboptimal one
2. This problem could be solved via an approach similar matching pursuit (but the full enumeration of the dictionary of bases needs to avoided, because it is too large in my problem)
3. Sebastian Reichelt provided code for an approach based on maximum and minimum subarray search


All and all I ended up with four different solutions to the presented problem:


1. `explore_a` indicates my original solution suggested in question (quantizing the $a\_i$ values and solving multiple maximum subarray problems)
2. `brute_force` indicates the approach based on approximating $\left|w' - w\right|$ by $\left(w' - w\right)^2$. In this approach, it turns out that each rectangle possibility can be evaluated with only two memory reads, one multiplication and one division, making a brute force exploration tractable.
3. `reichelt` indicates a port of Sebastian's code to python (which I used to test the ideas)
4. `abs_max` indicates an approach where at each iteration a rectangle of a single element, placed at the maximum absolute value of the signal, is selected.


All the [code implementing these methods](http://lts.cr/add) and some example results are available at <http://lts.cr/Hzg>


At each run the code will generate a new "random" signal (with a noisy step). An example signal (labeled `a[0:n]`) and the first rectangle selected by each method can be seen below.
[Example signal http://lts.cr/adL](http://lts.cr/adL)


Then each method is run recursively to approximate the input signal, until the approximation error is very low or the maximum number of iterations have been reached.


At typical result can be seen below (and another [here](http://lts.cr/adN))
[Example result 1 http://lts.cr/adT](http://lts.cr/adT)


After running the code multiple time (for different random signals), the behavior seems consistent:


1. Unsurprisingly, `abs_max` reaches convergence in a number of iterations equal to the signal length. It does so quite fast (exponential decay).
2. `explore_a` decrease the energy fast initially and then tends to stagnate.
3. `reichelt` is consistently worse than `explore_a`, getting stuck at a higher energy level. I hope I did not do any dumb mistake in the port from C to Python. By visual inspection the [first rectangle](http://lts.cr/adL) [selected](http://lts.cr/adP) seems reasonable.
4. `brute_force` is consistently the method that decreases the energy the fastest. It is also consistently intersects the `abs_max` solution, which indicates that a better strategy would be to switch from one method to the other.


Obviously the exact behavior changes from run to run and would certainly change depending on the method used to generate the "random" signal. However I believe these initial results are a good indicator. I feel that it is reasonable now to proceed to generate my real data, and evaluate how well/fast `brute_force` and `explore_a` run there.


Feel free to add comments or play around with the code. 


Again, thank you very much for your inputs and insights ! 


