# Deconvolution (Linear Convolution) with an Under Determined System of Equations?

If I have a measured signal $\mathbf{y}$ which is the result of the true signal $\mathbf{x}$ convolved with another function (linear and not circular convolution), I always seem to get an underdetermined system. Convolution expressed as a matrix equation:

$$\mathbf{y} = \mathbf{A} \mathbf{x}$$

I am interested in finding $\mathbf{x}$, given $M=5$ observations and knowing that the convolution matrix is defined by a boxcar function of window length $n=3$. The corresponding $M\times N$ convolution matrix is defined as the following, where $N=M+n-1$. $$\mathbf{A} = \begin{pmatrix} 1 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 1 \\ \end{pmatrix}$$

In the literature, the convolution matrix is often $N\times N$ (or $M\times{}M$, however you want to define it). But it seems that with a finite number of observations, the convolution matrix corresponding to the available observations would have to be $M \times (M+n-1)$, which is always underdetermined (except when $n \equiv 1$, in which you don't have convolution of your signal).

In most convolution/deconvolution examples, the matrix is often represented as being "square," in which there exist rows in the matrix where the number of non-zero elements are less than $n$. This does not correspond to any possible measurement. Padding is common for the input signal ($\mathbf{x}$) but not for the output ($\mathbf{y}$), so how can the system not be underdetermined?

Edit for clarification:

I understand that convolution matrix is traditionally written as $N+m+1$:

$$\mathbf{A} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & 1 & 0 & 0 & 0 & 0 & 0\\ 1 & 1 & 1 & 0 & 0 & 0 & 0\\ 0 & 1 & 1 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 1 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 1 & 1 & 0\\ 0 & 0 & 0 & 0 & 1 & 1 & 1\\ 0 & 0 & 0 & 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}$$

In practice, you never measure the signals corresponding to the first two and last two rows, so it would seem that you cannot use this matrix for creating the deconvolution filter.

• You can still use the Least Squares solution. Yet not using the Normal Equations as it can't be inverted. You should use the SVD Least Squares solution which always exists (It brings the Least Squares and the Least Norm solution). – Royi Sep 25 '17 at 12:56
• Yes, the deconvolution problem is underdetermined in the general case, which is why regularized approaches are usually applied. The unregularized inverse convolution (division in the Fourier domain by the kernel) leads to, in general, useless results. – Cris Luengo Jul 14 '19 at 14:16

On the contrary, if you compute the convolution of a signal $x$ having a span of $M$ non-zero entries with a filter with a span of $N$ non-zero entries, then the resulting sequence will have $M+N-1$ non-zero entries, so that your system should be over-determined.

Computing the minimizer $x$ of $\|[0,1,0]*y-[1,1,1]*x\|_2$ leads by standard variational methods to the system of equations

$$[0,1,1,1,0]*y=[1,2,3,2,1]*x$$

with the first two and the last two removed, and for some reason mathematics / linear algebra claims that this system is solvable even if $[0,1,0]*y=[1,1,1]*x$ is not.

• I understand that this is how convolution is defined, but in practice, a series of measurements will not yield the additional $N-n$ measurements. – hatmatrix Jan 2 '14 at 2:35
• The observed blurred samples are taken from a window within a larger (infinite) domain. In the real world a convolution (e.g. blurring caused by a lens) is not applied to a set of samples within a finite domain and padded with zeros. It is applied to a longer signal and the result is cropped. This is what OP is asking about. – Cris Luengo Jul 14 '19 at 14:13

Your best bet is to use the least squares solution to inverting the matrix A. I won't go into detail of deriving it, but the wikipedia article is pretty good. The solution you would get is

$$\hat{\mathbf{x}} = \left(A^TA\right)^{-1}A^T\mathbf{y}$$

What does this solution mean? First, notice that we're solving for $\hat{\mathbf{x}}$, not $\mathbf{x}$, because there may not be a solution that solves your system perfectly. However, if you look at the difference vector $\mathbf{d} = \mathbf{x} - \hat{\mathbf{x}}$, then this solution guarantees that $\mathbf{d}$ has the smallest possible length out of all possible solutions in the usual Euclidean sense on length.

• Is not the least-norm solution the one preferred for underdetermined systems? – hatmatrix Jan 2 '14 at 15:46
• But this would also not take advantage of the Toeplitz structure of the convolution matrix in the final solution? – hatmatrix Jan 2 '14 at 16:11
• @hatmatrix It could if you wanted to. The Toeplitz structure of $A$ will impose certain nice properties on $A^TA$, so it would be easy to invert. Many numerical tricks should still hold, since there are many ways to compute LMS inversion problems. – Phonon Jan 3 '14 at 3:55
• sorry what is 'LMS'? And would it not be $\hat{\mathbf{x}}=A^T(AA^T)^{-1}\mathbf{y}$ for the least norm, rather than least squares solution? – hatmatrix Jan 3 '14 at 18:57
• Why least squares rather than least norm solution in this case? I have an undetermined system. – hatmatrix Jan 4 '14 at 13:33

You can always use the Least Squares Solution:

$$\arg \min_{x} \frac{1}{2} \left\| A x - y \right\|_{2}^{2}$$

Now, in cases $A$ isn't full rank you won't be able to solve it using the Normal Equations.
What you should do (And always works for Least Squares Solution) is use the Pseudo Inverse of $A$.

In MATLAB it will go like this:

vX = pinv(mA) * vY;


Nice property of this solution is the returned answer is both the Least Squares Solution and the Least Norm Solution (As the system has infinite number of solutions).