Solving inverse problem using black box implementation of the kernel

Question

My question is related to Solving regularized least squares problem using black-box computation of $\mathbf{A}\mathbf{x}$ and $\mathbf{A}^T\mathbf{x}$.

In case, the problem is formulated as:

\begin{align} \min_{\mathbf{x}} \frac{1}{2} \| \mathbf{y} - \mathbf{h} \ast \mathbf{x}\|^2 + \lambda \phi(\mathbf{x}) \end{align}

I would like to solve it using fast implementation of the convolution with the kernel $\mathbf{h}$.

In my case I am using regularization of the form:

• Ridge: $\phi(\mathbf{x}) = \| \mathbf{L} \mathbf{x} \|_{2}^{2}$.
• TV: $\phi(\mathbf{x}) = \| \mathbf{L} \mathbf{x} \|_{1}$.

Is there an efficient solver in that case?
When would you chose which regularization?

Answer 1

I will try to illustrate a solution.

The general form of the problem is:

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

There 2 variants to tackle this question using the ADMM framework:

1. Modeling the problem as (Model 001):

$$ \arg \min_{\boldsymbol{x}} f \left( \boldsymbol{x} \right) + g \left( \boldsymbol{x} \right) $$

2. Modeling the problem as (Model 002):

$$ \arg \min_{\boldsymbol{x}} f \left( \boldsymbol{x} \right) + g \left( \boldsymbol{Lx} \right) $$

Each model will generate a different sub problems.
Also, the ADMM requires calculating the Proximal operator of the functions.
For instance, it is hard to solve the proximal of $ {\left\| L \boldsymbol{x} \right\|}_{1} $ (Requires basically running a sub algorithm as there is no closed form) while solving $ {\left\| \boldsymbol{x} \right\|}_{1} $ is trivial (Hence model 2 is better for this case).

In general the Scaled ADMM optimization will generate something along:

$$ \begin{align} \arg \min_{\boldsymbol{x}, \boldsymbol{z}} \quad & f \left( \boldsymbol{x} \right) + g \left( \boldsymbol{z} \right) \\ \text{subject to} \quad & \begin{aligned} P \boldsymbol{x} + Q \boldsymbol{z} - \boldsymbol{c} = \boldsymbol{0} \\ \end{aligned} \end{align} $$

Then the solution is given by iterating:

1. $ \boldsymbol{x}^{k + 1} = \arg \min_{\boldsymbol{x}} f \left( \boldsymbol{x} \right) + \frac{\rho}{2} {\left\| P \boldsymbol{x} + Q \boldsymbol{z}^{k} - \boldsymbol{c} + \boldsymbol{u}^{k} \right\|}_{2}^{2} $.
2. $ \boldsymbol{z}^{k + 1} = \arg \min_{\boldsymbol{z}} g \left( \boldsymbol{z} \right) + \frac{\rho}{2} {\left\| P \boldsymbol{x}^{k} + Q \boldsymbol{z} - \boldsymbol{c} + \boldsymbol{u}^{k} \right\|}_{2}^{2} $.
3. $ \boldsymbol{u}^{k + 1} = \boldsymbol{u}^{k} + P \boldsymbol{x}^{k + 1} + Q \boldsymbol{z}^{k + 1} - \boldsymbol{c} $.

In the models we'll use the matrix form of the convolution, namely $ \boldsymbol{h} \ast \boldsymbol{x} \Rightarrow H \boldsymbol{x} $.
This basically makes the question equivalent to the linked one. We'll deal with using the convolution later on.

Model 001

In this model we set $ \boldsymbol{x} = \boldsymbol{z} \implies \boldsymbol{x} - \boldsymbol{z} = 0 $.
Basically we just want solve function on its own.
Hence the iterations becomes:

1. $ \boldsymbol{x}^{k + 1} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{\rho}{2} {\left\| \boldsymbol{x} - \boldsymbol{z}^{k} + \boldsymbol{u}^{k} \right\|}_{2}^{2} = {\left( {H}^{T} H + \rho I \right)}^{-1} \left( {H}^{T} \boldsymbol{y} + \rho \left( \boldsymbol{z}^{k} - \boldsymbol{u}^{k} \right) \right) $.
2. $ \boldsymbol{z}^{k + 1} = \arg \min_{\boldsymbol{z}} \lambda \phi \left( \boldsymbol{z} \right) + \frac{\rho}{2} {\left\| \boldsymbol{x}^{k} - \boldsymbol{z} + \boldsymbol{u}^{k} \right\|}_{2}^{2} $.
3. $ \boldsymbol{u}^{k + 1} = \boldsymbol{u}^{k} + \boldsymbol{x}^{k + 1} - \boldsymbol{z}^{k + 1} $.

The tricky part in this model is step (2). We have 2 variants:

• $ \phi \left( \boldsymbol{x} \right) = \frac{1}{2} {\left\| L \boldsymbol{x} \right\|}_{2}^{2} \implies g \left( \boldsymbol{z} \right) = \frac{\lambda}{2} {\left\| L \boldsymbol{z} \right\|}_{2}^{2} $
  Then $ \boldsymbol{z}^{k + 1} = \frac{\rho}{\lambda} {\left( {L}^{T} L + \frac{\rho}{\lambda} I \right)}^{-1} \left( \boldsymbol{x}^{k} + \boldsymbol{u}^{k} \right) $.
• $ \phi \left( \boldsymbol{x} \right) = {\left\| L \boldsymbol{x} \right\|}_{1} \implies g \left( \boldsymbol{z} \right) = \lambda {\left\| L \boldsymbol{z} \right\|}_{1} $
  Then there's no closed form solution but basically we have to solve something similar to the main problem just without the convolution with the kernel $ \boldsymbol{h} $.

So this model fits the case the regularization is the ridge regularization.

Model 002

In this model we set $ L \boldsymbol{x} = \boldsymbol{z} \implies L \boldsymbol{x} - \boldsymbol{z} = 0 $.
Basically this will help use to reduce the regularization to work directly on the input vector and not a linear mapping of it.
Hence the iterations becomes:

1. $ \boldsymbol{x}^{k + 1} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| H \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{\rho}{2} {\left\| L \boldsymbol{x} - \boldsymbol{z}^{k} + \boldsymbol{u}^{k} \right\|}_{2}^{2} = {\left( {H}^{T} H + \rho {L}^{T} L \right)}^{-1} \left( {H}^{T} \boldsymbol{y} + \rho {L}^{T} \left( \boldsymbol{z}^{k} - \boldsymbol{u}^{k} \right) \right) $.
2. $ \boldsymbol{z}^{k + 1} = \arg \min_{\boldsymbol{z}} \lambda \phi \left( \boldsymbol{z} \right) + \frac{\rho}{2} {\left\| L \boldsymbol{x}^{k} - \boldsymbol{z} + \boldsymbol{u}^{k} \right\|}_{2}^{2} $.
3. $ \boldsymbol{u}^{k + 1} = \boldsymbol{u}^{k} + L \boldsymbol{x}^{k + 1} - \boldsymbol{z}^{k + 1} $.

The (1) step became more complex but still doable using convolutions (2).
Yet step (2) becomes a bit simpler:

• $ \phi \left( \boldsymbol{x} \right) = \frac{1}{2} {\left\| L \boldsymbol{x} \right\|}_{2}^{2} \implies g \left( \boldsymbol{z} \right) = \frac{\lambda}{2} {\left\| \boldsymbol{z} \right\|}_{2}^{2} $
  Then $ \boldsymbol{z}^{k + 1} = \frac{\rho}{1 + \lambda} \left( L \boldsymbol{x}^{k} + \boldsymbol{u}^{k} \right) $.
• $ \phi \left( \boldsymbol{x} \right) = {\left\| L \boldsymbol{x} \right\|}_{1} \implies g \left( \boldsymbol{z} \right) = \lambda {\left\| \boldsymbol{z} \right\|}_{1} $
  There is a simple closed form solution by the Soft Threshold operator: $ \boldsymbol{z}^{k + 1} = \mathbf{S}_{ \frac{\lambda}{\rho} } \left( L \boldsymbol{x}^{k + 1} + \boldsymbol{u}^{k} \right) $ where $ \mathbf{S}_{\lambda} \left( \cdot \right) $ is the Soft Threshold operator with parameter $ \lambda $.

So in this model we can solve pretty efficiently for both regularization variants.

Solving Using Convolution Operator

I assume the motivation to the question is avoiding working with the explicit matrices $ H $ and $ L $.
I assume we have them in the form of convolution kernels.

Once we work with convolution of discreet signals we should choose the type of convolution (Mostly how to handle the boundaries). In case of cyclic convolution (Periodic boundary condition) all is easy to implement using the DFT Matrix.

For the general matrix we can use another approach.
If $ \boldsymbol{h} \ast \boldsymbol{x} = H \boldsymbol{x} $ then the operator $ {H}^{T} $ is basically applying the same kernel using correlation instead of convolution (See my answer to Meaning of the Transpose of Convolution).
So we basically have a fast implementation of both $ H $ and $ {H}^{T} $.

We have 3 equations to solve (In the 2 models):

1. $ {\left( {H}^{T} H + \rho I \right)} \boldsymbol{x} = \left( {H}^{T} \boldsymbol{y} + \rho \left( \boldsymbol{z}^{k} - \boldsymbol{u}^{k} \right) \right) $.
2. $ {\left( {L}^{T} L + \frac{\rho}{\lambda} I \right)} \boldsymbol{z} = \frac{\lambda}{\rho} \left( \boldsymbol{x}^{k} + \boldsymbol{u}^{k} \right) $.
3. $ {\left( {H}^{T} H + \rho {L}^{T} L \right)} \boldsymbol{x} = \left( {H}^{T} \boldsymbol{y} + \rho {L}^{T} \left( \boldsymbol{z}^{k} - \boldsymbol{u}^{k} \right) \right) $.

In all cases the left hand side matrix is symmetric and positive definite.
In case we drop the option of working with matrices, One, relatively efficient, iterative solver is the Conjugate Gradient method which takes advantage of the properties mentioned above.
It also only requires the forward method, namely all needed are the following operators (The right hand side can be calculated by correlation):

• $ {\left( {H}^{T} H + \rho I \right)} $.
• $ {\left( {L}^{T} L + \frac{\rho}{\lambda} I \right)} $.
• $ {\left( {H}^{T} H + \rho {L}^{T} L \right)} $

Remark: @mlbj, in his other question (Matrix Vector Multiplication Representation of Total Variation Function) found an interesting set of algorithms (Marcelo V. W. Zibetti, Chuan Lin, Gabor T. Herman - Total Variation Superiorized Conjugate Gradient Method for Image Reconstruction). They seem to be fast, yet the quality of their result seems lower than FISTA (See Figure 4 in the paper). The ADMM above is usually even faster and better than FISTA.