A notion of resolution in inverse problems

Question

Suppose I have a linear inverse problem of the form:

\begin{align} Ax=b \end{align}

I would like to reconstruct $x$ from the measurement $b$ via the objective

$$\min_x\{\vert\vert Ax-b\vert\vert^2_2+\lambda\mathcal{R}(x)\}$$

With $\mathcal{R}(x)$ being the regularizer and $\lambda$ being the regularization parameter.

If we chose no regularization, i.e. $\mathcal{R}(x)=0$, the optimum solution will be achieved when applying the (generalized) inverse operator $A^\dagger$ to the data $b$. At the same time, it is given that:

$$A^\dagger A = I$$

In reality, the full inverse operator is not used, but only an approximate inverse, in order to not fit noise present in $b$. This is either achieved by truncation at an iteration before noise is fitted (Gradient Descent or Truncated Singular Value Decomposition) or if the regularizer is differentiable, e.g. $\mathcal{R}(x) = \vert\vert x\vert\vert^2_2$ (Tikhonov regularization), while setting $\lambda$ appropriately.

In these cases, we are able to find analytical expressions for the inverse $A^\dagger$ and hence the matrix product $A^\dagger A$. We further denote the approximate inverse via a subscript $k$, i.e. $A_k^\dagger$.

The matrix product

$$A_k^\dagger A = R_k$$

is called the resolution matrix (also called Backus-Gilbert Resolution Kernel). This matrix can tell us about the spread of a single parameter in $x$ in relation to the resulting values in $b$ (in imaging this would tell us about the Point-Spread-Function). Plotting a column of the matrix then yields a localized peak, which can be used for a heuristic assignment of resolution via e.g. its full width half maximum (FWHM) for the corresponding parameter.

My Question: Suppose the problem involves a non-differentiable regularizer, such as an $\ell^1$ norm regularization. Can one then also come up with a similar assignment of resolution in such a case?

Answer 1

You need to go back to why it is that a column of $A_k^\dagger A$ gives you the point-spread function. In essence, the question you're trying to answer is this: If I changed a single entry of $b$ by a unit change, how would that change propagate to $x$?

For the quadratic optimization problem you outline, you have $$ x = B b $$ with some matrix $B$ that comes out of $A$ and the regularization matrix in the appropriate way, and so if you had a change $\tilde b = b + e_i$ then you'd get a solution $$ \tilde x = x + \delta x = B \tilde b = Bb + B e_i. $$ In other words, the change is $$ \delta x = B e_i, $$ so looking at the $i$th column of matrix $B$ makes sense -- because that's just the right hand side $Be_i$.

For non-quadratic regularizers, there is no such simple matrix $B$ that produces the map from $b$ to $x$. Instead, it will be a nonlinear operator: $$ x = B(b). $$ But you can compute $B(b)$ for a given $b$, and so you can also compute $$ \delta x = \tilde x - x = B(b+e_i) - B(b). $$ Looking at how fast the entries of $\delta x$ decay away from entry $i$ gives you exactly the same information as the $i$th column of the matrix in the quadratic case. Of course, you now have to say what your "background" $b$ is because this kind of resolution information depends on your linearization point $b$, but the idea is still valid: $b$ is likely going to be something like a "reference" or "noiseless" point.