How to rebuild a $N\times N$ matrix from the sums of all of its submatrices of size $M\times M$ ($M\ll N$)

What I have is sums of smaller submatrices of size $M\times M$ ($M$ is much smaller than $N$, say $N/6$ or less). The sums of all possible submatrices at all positions are known. I am trying to rebuild the entire $N\times N$ matrix.

For example, if there is a $4\times 4$ matrix $X$ with the values

$$\begin{array}{cccc} 1 & 2 & 3 & 4\\ 5 & 6 & 7 & 8\\ 9 & 1 & 2 & 3\\ 4 & 5 & 6 & 7\\ \end{array}$$

And I have the sums of submatrices of size $2\times 2$ (for all possible positions - assuming that the matrix is surrounded by infinite number of zero values).

The first known submatrix sum would thus be $0+0+0+1=1$: $$\hskip-0.3in \begin{array}{ccccc} \fbox{\begin{matrix} 0 & 0 \\ 0 & 1\end{matrix}}\!\!\!\! & {\atop 2} & {\atop 3} & {\atop 4}\\ \begin{matrix} \hphantom{0} & 5\end{matrix}\!\!\!\! & 6 & 7 & 8\\ \begin{matrix} \hphantom{0} & 9\end{matrix}\!\!\!\! & 1 & 2 & 3\\ \begin{matrix} \hphantom{0} & 4\end{matrix}\!\!\!\! & 5 & 6 & 7\\ \end{array}$$

...and the second submatrix sum, $0+0+1+2=3$:

$$\begin{array}{cccc} \fbox{\begin{matrix} 0 & 0 \\ 1 & 2\end{matrix}}\!\!\!\!\!\! & {\atop 3} & {\atop 4} \\ \begin{matrix} 5 & 6\end{matrix}\!\!\!\!\!\! & 7 & 8\\ \begin{matrix} 9 & 1\end{matrix}\!\!\!\!\!\! & 2 & 3\\ \begin{matrix} 4 & 5\end{matrix}\!\!\!\!\!\! & 6 & 7\\ \end{array}$$

The third sum would be $5$, the fourth $7$, and so on for each row and column in $X$. There are $(N+1)\times(N+1)$ of these sums, and they can be written as the matrix $Y$ (containing all available submatrix sums):

$$\begin{matrix} 1 & 3 & 5 & 7 & 4\\ 6 & 14& 18&22 & 12\\ 14 & 21 & 16&20 & 11\\ 13 & 19& 14& 18&10 \\ 4 & 9& 11&13 & 7 \end{matrix}$$

The question is, how to rebuild / calculate the original $N\times N$ matrix $X$ from the $(N+1)\times(N+1)$ matrix $Y$, from the available submatrix sums? If it is not possible to calculate exactly, how well can it be approximated?

Any hints are appreciated!

-

Wait a minute, if your matrix is padded with infinite zeros, then this is really easy.

You have $$Y(i,j) = \sum_{\substack{i'=0\text{ to }M-1 \\ j'=0\text{ to }M-1}} X(i-i',j-j'),$$ so $$X(i,j) = Y(i,j) - \sum_{\substack{i'=0\text{ to }M-1 \\ j'=0\text{ to }M-1 \\ (i,j) \neq (0,0)}} X(i-i',j-j').$$ That just says that for any $M\times M$ submatrix, you can find the value at the bottom right corner if you know all the other entries.

So start with the upper left corner: clearly, $X(1,1) = Y(1,1)$ because all the rest of the entries are zero as they are outside the $N\times N$ matrix. Now that you know $X(1,1)$, you can find $X(1,2) = Y(1,2) - X(1,1)$, and then $X(1,3) = Y(1,3) - \ldots$, and thus fill in the whole first row of $X$. Then you do the same thing, marching down the rest of the rows, and you're done.

-
Well, what can I say :-) Tamas, I suggest you accept this answer instead of mine; I was using a sledgehammer to crack a nut... –  joriki Aug 7 '11 at 19:01
Thanks Rahul! Thats very useful. Would you also know a solution (or an approximation) in case there are unknown values around the matrix X, instead of zeroes? –  Tamas Madl Aug 7 '11 at 19:57

You have effectively applied a box blur filter. This is diagonal in Fourier space, so you can undo it by performing a Fourier transform, dividing by the filter's frequency response, and transforming back. For that to work, you need to choose a transform size such that none of the Fourier modes is annihilated by the filter (else you can't reconstruct them). This will be the case e.g. if the transform size is a sufficiently large prime. Since I don't know how much you already know about these things, I'll stop here; feel free to ask if any of that isn't clear.

-
Thank you very much for the quick reply, joriki! A few questions: What is a filter frequency response? Is the size of the submatrices M the transform size? And finally, do you have a website or paper where I can read up on the subject? (Unfortunately I know little about this stuff...) –  Tamas Madl Aug 7 '11 at 16:56
@Tamas: This looks like it contains a lot of what you might want to read up on: cs.toronto.edu/~jepson/csc2503/linearFilters.pdf. The frequency response is the (complex) factor by which the filter multiplies each of the Fourier components. No, the transform size must be at least $N+1$, plus some padding to avoid artefacts from the periodicity (I think $2$ or $3$ should be enough in this case; just increase it until the result doesn't change), plus potentially some padding to avoid one of the Fourier modes being annihilated by the filter. –  joriki Aug 7 '11 at 18:55
Thanks for the pdf, joriki, it certainly helped - however, its still not clear to me how i can do this in my case, the case of oversampling (as far as i have understood the aliasing and area sampling reconstruction in the pdf, and in other sources i have looked at, they only work with undersampling, since that is the case in an overwhelming majority of image processing and similar applications). Could you maybe provide a short approach about how i could apply fourier analysis to this problem with oversampling? (Sorry if there was something in the pdf addressing that which I have overlooked...) –  Tamas Madl Aug 7 '11 at 20:00
@Tamas: I think all of this has been superseded by Rahul's much simpler solution, but if you're still interested anyway: There's no question of over- or undersampling; you're not sampling. You have some data, and you transform them to Fourier space and then back -- the Nyquist theorem theorem doesn't enter into it; all you have to make sure is that you have enough padding so the two boundaries of the matrix don't interact with each other, and that you choose a transform size that doesn't cause one of the Fourier modes to be annihilated; in the $2\times2$ case, that means the transform size... –  joriki Aug 7 '11 at 20:50
... should be odd. So the recipe is basically: Choose appropriate transform size (in the example above, $7$ should be good), DFT in $x$ direction, DFT in $y$ direction, divide by the filter's frequency response, inverse DFT in $y$ direction, inverse DFT in $x$ direction. I linked to a PDF that concentrated on image manipulation because I thought it would have the most appropriate stuff for your 2D calculations (instead of 1D in music or 3D in solid state physics, for example), but some of that image-related stuff isn't relevant here. This is just about transforming to eigenmodes of the filter. –  joriki Aug 7 '11 at 20:52