Solving linear equation of discrete convolution kernels using black box model for the convolution

In Solving inverse problem using black box implementation of the kernel the solution depends on solving the equations of the form:

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

Where $$H$$ and $$G$$ represent convolution operation of the $$h$$ and $$g$$ kernels.

Given we have a a black box implementation of the convolution by $$h$$ and $$g$$, how can we solve those equations efficiently?

Those forms are very common in many image processing algorithms and solving them by the matrix representation is not efficient.

• The inverse operation is called deconvolution. Classical method are Richardson-Lucy and Iterative Constrained Tikhonov-Miller. There are some newer solutions, such as Fast Iterative Shrinkage Thresholding (FISTA). Commented Apr 14, 2023 at 23:56
• @CrisLuengo, If you look at the context he used, the idea is as following, if you have a kernel $\boldsymbol{h}$ and $\boldsymbol{g}$, how do you apply above fast using convolution operations.
– Royi
Commented Apr 15, 2023 at 5:52
• @Royi If you look up the above methods, they're always implemented using the FFT. There are lots of implementations out there to use (or just to look at if OP insists in implementing it themselves). Commented Apr 15, 2023 at 7:31
• @CrisLuengo, If you look on the context of the answer it is something like the FISTA model (TV regularization) yet with ADMM which is even faster than the accelerated proximal gradient descent in FISTA.
– Royi
Commented Apr 15, 2023 at 16:57

This is not an easy question, not because ideas are complex but since the details matter a lot.

As I wrote in my answer in the question Solving Inverse Problem Using Black Box Implementation of the Convolution Kernel, the issue is defining how to handle the boundaries (See The Different Solutions for Filter Coefficients Estimation for Periodic Convolution and Full Convolution). This is super important for images.

Matrix Case

In all cases, for this equation:

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

Can be solved using sparse solvers taking advantage of the sparsity pattern of the matrices and the matrix being symmetric positive definite.

This will work for any variant of the convolution.
Though we have 2 options, direct solvers and iterative ones.

The iterative ones allow us taking advantage of the convolution black box as they usually require applying the operator $$A$$. Many iterative solvers require $${A}^{T}$$ but if they do, it means they don't model it as a symmetric operator. Hence it means they are not optimal. So stick with those assuming symmetric operator, hence require $$A$$ only.

Now the question, how to implement $$A$$ efficiently?

Periodic / Cyclic Boundary Conditions

If we assume that the $$\ast$$ (Convolution) is applied using cyclic boundary convolution we can apply everything in frequency domain.

The Convolution of Type full / same / valid

For the other convolution types (I Use the lingo of MATLAB's convolution functions, such as conv2()) things are trickier.

I created the matrix $${H}^{T} H$$ for the different convolution shapes from the same 3 coefficients kernel applied on 5 samples signal:

As one can see, while all cases generates a 5x5 matrix, yet only the full case generates a Toeplitz Matrix.
Hence, only this case can be represented by a convolution kernel with some convolution type.

General Solution

We can always apply each matrix by itself. Something like:

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

So $$H \boldsymbol{x}$$ / $$G \boldsymbol{x}$$ is just applying the convolution, what about $${H}^{T} \boldsymbol{z}$$ / $${G}^{T} \boldsymbol{z}$$? Well, we can have a recipe per case, where the adjoint is done by the flipped kernel:

• Convolution Type full -> Correlation type valid.
• Convolution Type same -> Correlation type same.
• Convolution Type valid -> Correlation type full.

Now we can apply all of the above one by one.

Optimized Solution for full

As we can see, in case our convolution is of type full the normal matrix $${H}^{T} H$$ is a Toeplitz Matrix, so we can generate the composite operation using a single kernel convolution.

Results

I implemented all cases above and using Conjugate Gradient I got the following results for 1D:

As can be seen, we can get the same results using the iterative method.
This method only applies convolution operations and doesn't require the matrix form.

Some Other Resources

• You can get around the periodic boundary condition of the FFT-based convolution by simply padding the image a bit. This is a fairly standard procedure, I think? Commented Apr 15, 2023 at 7:32
• $H^TH$ translates to $H^*H = |H|^2$ if $H$ is the Fourier transform of the kernel, rather than the convolution matrix. Commented Apr 15, 2023 at 7:35
• I wrote that I don't take care of the DFT in this answer as I wrote about it in many other places (I linked, some include the extension trick). I am not sure it will be faster from using the full method. I am working on extending the answer.
– Royi
Commented Apr 15, 2023 at 16:55
• @CrisLuengo, Does you library have something like that pre defined as a function? I am talking about generating this in Fourier domain using the proper padding and afterwards proper extraction of the signals?
– Royi
Commented Apr 16, 2023 at 7:21
• @CrisLuengo, OK. Let's see if the OP is interested in it. Then I might compare high level results to your code (I read it has Python / MATLAB wrapper, I don't "speak" C++ :)).
– Royi
Commented Apr 16, 2023 at 14:36