# Exact deconvolution of two matrices using numerical techniques

Suppose that I am given two $n \times m$ matrices $\bf{A}$ and $\bf{C}$, and let $\bf{B}$ be a matrix that is convolved with $\bf{A}$, such that:

$\bf{A} * B = C$

In the above, $*$ is the convolution operator (link).

I am searching for a numerical procedure (and ideally a reference book or paper on the mathematics of the operation) that will allow me to determine $\bf{B}$, given $\bf{A}$ and $\bf{C}$.

I know both $\bf{A}$ and $\bf{C}$ exactly. The $\bf{A}$ is the system input, whereas the $\bf{C}$ is the system output.

Although the convolution is a product in the Fourier domain, I might be able to determine $\bf{B}$ by division. However, using the 2D Fast Fourier Transform (FFT), do I have to pad the matrices in the spatial domain or deal with division by zero? How do I do this in an elegent manner?

Is there a way to do this in Matlab using the fft2 and ifft2 functions?

-

In your case, there is no noise in the output, then the Fourier method is a safe way for deconvolution: first zero-pad $\mathbf{A}$ to the size of $\mathbf{C}$ (without zero-padding will leads to result corresponding to circulant convolution); then compute the spectra using fft2; divide the spectrum of $\mathbf{C}$ by the one of $\mathbf{A}$, then you can get the deconvolution result by using ifft2.

However, if the output contains noise, then the Fourier method may not a proper way for deconvolution, since division in frequency domain could amplify the noise and contaminate the results. The better way is, as Seyhmus Güngören mentioned, to convert the deconvolution problem into a system of linear equations or optimization problem with certain regularizations, then use iterative algorithms to solve the equations or optimization problem.

-
Thank you. (1) If $\bf{A}$ is the same size as $\bf{C}$, then is zero-padding required? (2) Do $\bf{A}$ and $\bf{C}$ have to be both zero-padded to size $(2n - 1) \times (2m - 1)$ to prevent circular convolution? –  Nicholas Kinar Sep 2 '12 at 23:56
If the convolution is linear, then size of $\mathbf{C}$ should be $(n_A+n_B-1) \times (m_A+m_B-1)$ provided $\mathbf{A}$ is $n_A \times m_A$ and $\mathbf{B}$ is $n_B \times m_B$. Also, if it's linear convolution, $\mathbf{C}$ does not to be zero padded. –  chaohuang Sep 3 '12 at 0:25
Thanks again! So I need to pad all matrices $\bf{A},\bf{B},\bf{C}$ in the spatial domain to be $(n_A+n_B-1) \times (m_A+m_B-1)$, the same size as $\bf{C}$? –  Nicholas Kinar Sep 3 '12 at 0:31
You just need to zero pad $\mathbf{A}$ to the size of $\mathbf{C}$, and because you want to recover $\mathbf B$, the result of ifft2 needs to be trimmed down to the size of $\mathbf B$ –  chaohuang Sep 3 '12 at 1:30
Thanks; this makes sense. I think that the cropping is similar to any convolution of two matrices. An example of this cropping (for Matlab) and a numerical algorithm for 2D FFT convolution can be found on pg. 178 of R. C. Gonzalez et al., Digital image processing using MATLAB, 2nd ed. S.l.: Gatesmark Publishing, 2009. –  Nicholas Kinar Sep 3 '12 at 2:30

Yes you need to make zero padding if you intend to do everything in the DFT domain. However there is an obvious risk of division with zeros. There are iterative algorithms to get rid of division with zero. These algorithms are not specific for the applications in the DFT domain.

I am not completely sure but I think it might be possible to convert the problem to a $1D$ convolution, reshaping the matrices as $1D$ vectors. Then one can use MATLAB with deconv command.

Another way is to write codes for your implementation. Basically you have a linear system of equations with some unknowns. If you can write them in a suitable form what remains is to solve the system. It is simply the multiplication of the inverse of the matrix at the left side with the matrix at the right side. But one should construct these matrices before.

-
The linear system is a good idea. –  Nicholas Kinar Sep 2 '12 at 23:57
(1) Do I write the linear system of equations in the frequency domain or the spatial domain? (2) Is there an equation showing the numerical procedure that I have to do? –  Nicholas Kinar Sep 2 '12 at 23:59
Usually it's in spatial domain. –  chaohuang Sep 3 '12 at 0:27
OK, I think that I understand now. –  Nicholas Kinar Sep 3 '12 at 0:32
I will write you more today. DCT had to be DFT. Linear system of equations in the spatial domain. –  Seyhmus Güngören Sep 3 '12 at 7:31

An extremely-readable paper on the mathematics of the deconvolution problem [1] discusses deconvolution for 1D vectors and 2D matrices. The paper is written from the perspective of an applied mathematician, and recasts the deconvolution problem as a special case of Fredholm integral equations. It discusses periodic and nonperiodic deconvolution problems, outlines a framework for deconvolution without using the FFT, and is probably one of the most exciting papers that I have ever read.

[1] P. C. Hansen, “Deconvolution and Regularization with Toeplitz Matrices,” Numerical Algorithms, vol. 29, no. 4, pp. 323–378, 2002.

-
This paper also introduces blind deconvolution as well, and contains much of the theory (and a number of references) involving the mathematics of deconvolution. It is a good theoretical reference. –  Nicholas Kinar Sep 3 '12 at 2:32
Blind deconvolution is useful when the kernel $\bf{B}$ is not known exactly, and $\bf{A}$ needs to be recovered from $\bf{C}$. –  Nicholas Kinar Sep 3 '12 at 2:34