# How to Solve a Composition of Convolutions from Regularized Least Squares Model in Frequency Domain

Assume we need to solve the model:

$$\arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| \boldsymbol{h} \ast \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| \boldsymbol{g} \ast \boldsymbol{x} \right\|}_{2}^{2}$$

Where $$\ast$$ is the convolution operator.

How can we solve this in the frequency domain?
Assuming different convolution models (full, same and valid).

The question is derived from Solving Linear Equation of Discrete Convolution Kernels Using Black Box Model for the Convolution.

In order to solve the problem let's formulate it in its matrix form:

$$\arg \min_{\boldsymbol{x}} {\left( {H}^{T} H + \lambda {G}^{T} G \right)} \boldsymbol{x} = {H}^{T} \boldsymbol{y}$$

Where $$H$$ and $$G$$ are the matrix form of the convolutions.
Since in the frequency domain all we can apply is a circular / cyclic / periodic convolution we need to create a form of the convolutions using a circulant matrix.

We can do as following:

$$H = E {C}_{\boldsymbol{h}} P$$

Where:

• $$P$$ is the padding matrix, it adds elements to the argument.
• $${C}_{\boldsymbol{h}}$$ is a circulant matrix build on the convolution kernel of $$H$$.
• $$E$$ is the extracting matrix. It extract what's needed from the cyclic convolution.

For instance, assume input $$\boldsymbol{x} \in \mathbb{R}^{5}$$ and a kernel $$\boldsymbol{h} \in \mathbb{R}^{3} = {\left[1, 2, 3 \right]}^{T}$$.
If we want to apply a valid convolution using the above, then:

$$E = \begin{bmatrix} 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix}, \; {C}_{\boldsymbol{h}} = \begin{bmatrix} 1 & 0 & 0 & 3 & 2 \\ 2 & 1 & 0 & 0 & 3 \\ 3 & 2 & 1 & 0 & 0 \\ 0 & 3 & 2 & 1 & 0 \\ 0 & 0 & 3 & 2 & 1 \end{bmatrix}, \; P = {I}_{5}$$

So for any $$\boldsymbol{x} \in \mathbb{R}^{5}$$ the matrix operation $$E {C}_{\boldsymbol{h}} P \boldsymbol{x}$$ will be equivalent to $$\boldsymbol{h} \ast \boldsymbol{x}$$ with valid type of convolution.

• One can this in my implementation for Replicate MATLAB's conv2() in Frequency Domain.
• The valid / same types of convolution is not commutative. In all the above we assume the case where the kernel $$\boldsymbol{h}$$ is rolling on the signal $$\boldsymbol{x}$$. In order to make sense for vlaid we also assume the length of the signal is bigger than the length of the kernel.

Since $${C}_{\boldsymbol{h}}$$ is a circulant matrix is can be diagonalized using the DFT Matrix:

$${C}_{\boldsymbol{h}} = {F}^{H} {D}_{\boldsymbol{h}} F$$

Where:

• $${D}_{\boldsymbol{h}}$$ (Diagonal Matrix) has the DFT transform of $$\boldsymbol{h}$$ along its diagonal (Padded in zeros to the number of rows of $$P$$).
• $${F}$$ (Unitary Matrix) is the DFT Matrix.

Let's analyze $${H}^{T} H$$ from above:

$${H}^{T} H = {\left( E {C}_{\boldsymbol{h}} P \right)}^{T} \left( E {C}_{\boldsymbol{h}} P \right) = {\left( E {F}^{H} {D}_{\boldsymbol{h}} F P \right)}^{H} \left( E {F}^{H} {D}_{\boldsymbol{h}} F P \right) = {P}^{T} {F}^{H} {D}_{\boldsymbol{h}}^{H} F {E}^{T} E {F}^{H} {D}_{\boldsymbol{h}} F P$$

As can be seen above, the only way to have a closed form solution in the Frequency Domain is when $${E}^{T} E = I$$. This happens when the convolution type is full only.

Once that happens, we indeed get:

$${H}^{T} H = {P}^{T} {F}^{H} {D}_{\boldsymbol{h}}^{H} {D}_{\boldsymbol{h}} F P$$

Which basically means:

• Pad the vector to length length(vH) + length(vX) - 1 (With zeros).
• Transform vX into DFT. Namely the 2 steps are equal to DFT with zero padding.
• Multiply by the squared magnitude of the padded transform of the kernel.
• Go back into the original domain (Time / Spatial, etc..).

For the composition, we get:

$${H}^{T} H + \lambda {G}^{T} G = {P}^{T} {F}^{H} \left( {D}_{\boldsymbol{h}}^{H} {D}_{\boldsymbol{h}} + \lambda {D}_{\boldsymbol{g}}^{H} {D}_{\boldsymbol{g}} \right) F P$$

So the whole process is the same.
From here, solving the equation, in the case of full or periodic is trivial.

Unless I am missing a trick to by pass $${E}^{T} E$$, it can not be done for the other types of convolution.
This is unfortunate because classic deconvolution (Not based on Deep Learning) is best when we assume valid convolution.
Namely where vX is larger than the data vY. Then when we estimate vX we basically only look on its center and not edges (Estimated by far less data).

• Interesting thoughts. Thank you for sharing! As you know I think a bit differently about the deconvolution problem, and am not too keen on MATLAB’s conv implementation with “same” or “full” modes: padding the image with zeros is one of the worst possible ways to handle the boundary. And with “valid” mode you always end up with too small a result — apply a couple of filters in a row and there’s nothing left of the image! :) But this is interesting work nonetheless. Apr 18, 2023 at 14:11
• @CrisLuengo, Usually I don't use convolution (Toeplitz) matrices. I build matrices which are similar to replicate or something like that. But since those are not Toeplitz, you can use the the frequency domain to calculate them in the frequency domain. The matrix $P$ in the above can be any padding model (Mirror, Replicate, etc...). Unless the $E$ matrix is identity, the composite in the frequency won't match.
– Royi
Apr 18, 2023 at 14:15
• Sure. But my point of view is that it doesn’t matter if it matches or not. You need to find a suitable way to handle the boundary so that it has least influence in the output. That handling might be different depending on how you compute the convolutions. In the end, you need to estimate what the infinite-length signal before the convolution looked like, not replicate one particular discrete solution with a different convolution implementation. The optimization problem in your question at the top is posed in the continuous domain. Apr 18, 2023 at 14:32
• @CrisLuengo, Could you show such case? Again, I'd argue that just doing the brute force DFT division / multiplication with padding has a different effect than you imagine. It applies basically a kernel which is different from what you assume.
– Royi
Apr 18, 2023 at 14:33
• @CrisLuengo, To make it more concrete, let's select a small image and solve it with any padding you'd like and compare results. I'd argue that the result from the DFT division won't match the LS solution of the padding of your choice.
– Royi
Apr 18, 2023 at 14:34