# Decomposing a discrete signal into a sum of rectangle functions

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.

• You haven't been specific about what you mean by $|w-w'|$. Is this $| w - w' | = \sum_{i=1}^{n} | w_{i}-w'_{i} |$ or $| w - w' | = \max_{i=1, 2, \ldots, n} | w_{i}-w'_{i} |$ or something else? – Brian Borchers Apr 27 '11 at 13:54
• Indeed, thanks for the comment. I have clarified the question now. – rodrigob Apr 27 '11 at 15:45
• You've stated that your signal is discrete, which is hugely important. How long are your signals- hundreds of samples? thousands of samples? millions of samples? – Brian Borchers Apr 27 '11 at 17:33
• the length of $w$ will be in the order of $10^4$ – rodrigob Apr 27 '11 at 19:52

I think that essentially what you want to do is to find the projection of your function onto the vector space generated by the first vectors in the Haar wavelet. More specifically, the Haar wavelet with mother wavelet function $\psi(t)$ can be described as

$$\psi(t) \begin{cases} 1 & 0 \leq t < 1/2,\\\ -1 & 1/2 \leq t < 1,\\\ 0 &\text{otherwise.} \end{cases}$$

Its scaling function $\phi(t)$ can be described as $$\phi(t) = \begin{cases} 1 \quad & 0 \leq t < 1,\\\ 0 &\mbox{otherwise.} \end{cases}$$

The Haar systems denotes the set of Haar wavelets $$t \mapsto \psi_{n,k}(t)=\psi(2^n t-k)$$ with $n \in \mathbb{N}$ and $0 \leq k < 2^n$. In Hilbert space terms, this constitutes a complete orthogonal system for the functions on the unit interval.

My understading is that the $w'(t)$ that you are looking for, is the projection of your signal $w(t)$ onte the subspace generated by the first $\psi_{n,k}$.

• No, this doesn't really do it. The problem here is that your Haar wavelets $\Psi_{n,k}$ have fixed scaling and shifts. The original poster wants to use a limited number of functions but adjust the scaling factors and shifts. – Brian Borchers Apr 27 '11 at 13:48
• Brian is right, Haar wavelet does not answer the question, because it is decomposing the function in a different basis than the one I am looking for. – rodrigob Apr 27 '11 at 14:27
• This was the path 1 I mentioned in the question. – rodrigob Apr 27 '11 at 18:26

You can formulate this problem as a 0-1 mixed integer linear programming problem. Whether this formulation would be of any practical use depends a lot on the length of the signal $w$. You can also relax the integer linear programming formulation to get a convex optimization problem that might be much more easily solvable and that might yield a practically useful solution.

Suppose that you're given a signal $w$ of length $m$. That is, $w$ is a vector in $R^{m}$.

Your goal is to minimize $\| w - w' \|_1$.

Construct a 0-1 matrix $A$ with each column consisting of a run of 0's, followed by a run of 1's, followed by a run of 0's. These columns are the rectangular basis signals out of which your signal will be built. There are $n=C(m,2)$ such columns (simply choose the start and end of the run of 1's.) We'll express $w'$ in terms of these basis vectors by using the constraint $w'=Ax$.

Now, your problem can be written as

$\min \| w-w' \|_{1}$

subject to

$w'=Ax$

$\| x \|_{0} \leq N.$

Here $\| x \|_{0}$ denotes "the number of nonzero entries in the vector $x$." It's not really a p-norm, but this notation is commonly used by folks in compressive sensing.

At this point, there are a couple of options.

An exact formulation of the sparsity constraint can be done using $n$ 0-1 integer variables $y_{i}$, $i=1, 2, \ldots n$. Pick a huge constant $M$, add the constraints

$| x_{i} | \leq My_{i} \;\; i=1, 2, \ldots n,$

and write $\| x \|_{0} \leq N$ as

$\sum_{i=1}^{n} y_{i} \leq N.$

In practice this won't be a very useful formulation unless $m$ is quite small, since $n=C(m,2)$ grows as the square of $m$. Integer programming problems with thousands of variables aren't practically solvable in most cases, although if your $N$ is very small it might help in making the problem easier for a branch and bound solver.

An alternative is to relax the sparsity constraint as is commonly done in compressive sensing. You'd replace the 0-norm constraint on $x$ with a constraint like $\| x \|_{1} \leq \delta$ where $\delta$ would be adjusted until you got a suitably sparse solution. If you go this route you should probably normalize the columns of the $A$ matrix.

Theoretical results on compressive sensing by $L_{1}$ regularization require matrices with special properties (e.g. RIP) that this $A$ matrix simply won't have. Thus you shouldn't expect to be able to get any of the nice theoretical results to apply to your problem.

• Although I didn't include this in my answer, it's easy to translate the $\| w -w' \|_{1}$ objective and inequalities involving $| x_{i} |$ into linear programming constraints- this is standard textbook material. – Brian Borchers Apr 27 '11 at 17:20
• This proposal corresponds to path 2 in my question. I though of this, however in my setup $A$ in too big (w is actually 2D, so the entries of $A$ grow as $n^4$, in practice the columns of $A$ will be $\sim10^{6}$). Is it there a way to solve this without explicitly computing $A$ ? – rodrigob Apr 27 '11 at 18:31
• If your problems are that large, than the integer programming formulation is out of the question. However, in compressive sensing practice, solving problems with $~ 10^6$ variables is not at all uncommon. First note that you don't necessarily need to be able to store the $A$ matrix explicitly to solve this problem- there are iterative methods that only need the ability to do matrix-vector multiplies with $A$ and $A^{T}$. Even if you did have to store $A$ explicitly, general purpose LP solvers are often used to solve LP's with millions of variables. So, this may not be out of reach. – Brian Borchers Apr 27 '11 at 20:48
• Simple iterative schemes such as "orthogonal matching pursuit" are often used on very large problems in compressive sensing. Here, you would start by finding the one rectangle that does the most to reduce $\| w - w' \|_{1}$, then recursively repeat this process on what's left until you've got as many terms as you're willing to use. This won't give you a solution that is necessarily optimal, but it might work well in practice. – Brian Borchers Apr 27 '11 at 20:51
• It also occurs to me that since your signals are actually 2D, you need to be using 2D basis functions- this would greatly expand the dimensions of the problem. That being the case, I've got nothing better to offer than the iterative scheme mentioned in my previous comment. – Brian Borchers Apr 27 '11 at 20:54

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 and here.

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

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) Example result 1 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 selected 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.