As discussed in the comments, here is an idea for a greedy algorithm. Basically, the goal is to iteratively find the rectangle that will minimize |w-w'| in that particular situation. A combination of a maximum subarray and minimum subarray algorithm finds a rectangle that comes pretty close, if you use the average over the resulting subarray as $a_i$. The only problem is that values below that average (or above, for negative rectangles) should already count as negative (positive); I don't know if there is a simple way to solve this (and it's a bit difficult to explain, too). One way to reduce this problem is to trim the rectangle if that improves the result. Another is to limit the rectangle length, as longer rectangles tend to have a higher percentage of values below the average. There may be better options, possibly depending on the shape of the signal.
Anyway, I needed to experiment with this to figure out whether it actually works, so here is some code (in C, but should be easy to port):
/* Type of index into signal. */ typedef unsigned int Index; /* This should be a floating-point type. If you decide to use integers, change the division to round away from zero instead of towards zero. */ typedef double Value; /* Represents a rectangle at [start,end). */ typedef struct { Index start, end; Value height; } Rectangle; /* Trims the rectangle by removing samples whose absolute value will increase when preliminaryHeight is subtracted. */ static Index optimizeRectangle(const Value signal[], Rectangle *rectangle, Value preliminaryHeight) { Index result = 0; /* Divided into two cases for better efficiency. */ if (preliminaryHeight >= 0) { while (rectangle->start + 1 < rectangle->end) { Value value = signal[rectangle->start]; if (preliminaryHeight - value < value) { break; } rectangle->start++; rectangle->height -= value; result++; } while (rectangle->start + 1 < rectangle->end) { Value value = signal[rectangle->end - 1]; if (preliminaryHeight - value < value) { break; } rectangle->end--; rectangle->height -= value; result++; } } else { while (rectangle->start + 1 < rectangle->end) { Value value = signal[rectangle->start]; if (preliminaryHeight - value > value) { break; } rectangle->start++; rectangle->height -= value; result++; } while (rectangle->start + 1 < rectangle->end) { Value value = signal[rectangle->end - 1]; if (preliminaryHeight - value > value) { break; } rectangle->end--; rectangle->height -= value; result++; } } return result; } /* Finds the rectangle that seems most appropriate at the moment, and whose length is at most maxResultLength. Returns the original length of the rectangle, before optimization. This value should be used for maxResultLength in the next iteration. If it is zero, the entire signal was zero. */ static Index findRectangle(const Value signal[], Index signalLength, Rectangle *result, Index maxResultLength) { result->start = result->end = 0; result->height = 0; Value resultHeightAbs = 0; /* Start index and height of current maximum and minimum subarrays. */ Index posStart = 0, negStart = 0; Value posHeight = 0, negHeight = 0; for (Index index = 0; index < signalLength; index++) { /* If the length of a subarray exceeds maxResultLength, backtrack to find the next best start index. */ Index nextIndex = index + 1; if (nextIndex - posStart > maxResultLength && index > 0) { Index oldStart = posStart; posStart = index; posHeight = 0; Value newHeight = 0; for (Index newStart = index - 1; newStart > oldStart; newStart--) { newHeight += signal[newStart]; if (newHeight > posHeight) { posStart = newStart; posHeight = newHeight; } } } if (nextIndex - negStart > maxResultLength && index > 0) { Index oldStart = negStart; negStart = index; negHeight = 0; Value newHeight = 0; for (Index newStart = index - 1; newStart > oldStart; newStart--) { newHeight += signal[newStart]; if (newHeight < negHeight) { negStart = newStart; negHeight = newHeight; } } } /* Extend the subarrays. */ Value value = signal[index]; posHeight += value; negHeight += value; /* Throw away the entire maximum subarray if it is negative or zero. */ if (posHeight <= 0) { posStart = nextIndex; posHeight = 0; } else { /* Update result if current maximum subarray is better. */ if (posHeight > resultHeightAbs) { result->start = posStart; result->end = nextIndex; result->height = posHeight; resultHeightAbs = posHeight; } } /* Throw away the entire minimum subarray if it is positive or zero. */ if (negHeight >= 0) { negStart = nextIndex; negHeight = 0; } else { /* Update result if current minimum subarray is better. */ Value negHeightAbs = -negHeight; if (negHeightAbs > resultHeightAbs) { result->start = negStart; result->end = nextIndex; result->height = negHeight; resultHeightAbs = negHeightAbs; } } } Index resultLength = result->end - result->start; if (!resultLength) { /* Return now to avoid division by zero. */ return 0; } /* Trim rectangle. */ Value normalizedHeight; Index trimLength; do { normalizedHeight = result->height / (result->end - result->start); trimLength = optimizeRectangle(signal, result, normalizedHeight); } while (trimLength); /* Normalize height. */ result->height = normalizedHeight; return resultLength; } /* Subtracts a rectangle from a signal. */ static void subtractRectangle(Value signal[], const Rectangle *rectangle) { for (Index index = rectangle->start; index < rectangle->end; index++) { signal[index] -= rectangle->height; } } /* Decomposes a signal into rectangles. Stores at most resultLength rectangles in result, and returns the actual number of rectangles written. All rectangles are subtracted from the signal. */ unsigned int extractRectangles(Value signal[], Index signalLength, Rectangle result[], unsigned int resultLength) { Index rectangleLength = signalLength; for (unsigned int resultIndex = 0; resultIndex < resultLength; resultIndex++) { Rectangle *rectangle = &(result[resultIndex]); rectangleLength = findRectangle(signal, signalLength, rectangle, rectangleLength); if (!rectangleLength) { return resultIndex; } subtractRectangle(signal, rectangle); } return resultLength; }
The main function is extractRectangles
. It returns the rectangles as an array, but this can be adapted easily. Have fun.