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?

  • 1
    $\begingroup$ What is an "assignment of resolution"? $\endgroup$ – Federico Poloni Dec 17 '20 at 14:06
  • 2
    $\begingroup$ Anyhow, it seems to me that there is no (useful) analogue of a point spread function because with a non-L2 regularization the problem is not linear (i.e., the solution map $s$ is such that $s(b_1+b_2) \neq s(b_1) + s(b_2)$). $\endgroup$ – Federico Poloni Dec 17 '20 at 14:07
  • 3
    $\begingroup$ "Spike tests" are often used with nonlinear inversions. Start with a $\delta$ function like model, run it through the forward model, run that through the inversion process, and then see what comes out. This is not really equivalent to finding the resolution matrix because the result is specific to that particular spike model that you started with. $\endgroup$ – Brian Borchers Dec 17 '20 at 20:04

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.

  • $\begingroup$ Shouldn't the change $e_i$ be given by feeding a unit change $c_i$ into the forward operator, i.e. $e_i=Ac_i$? Because $B$ denotes only $A^\dagger$ in this case. $\endgroup$ – Ron Dec 20 '20 at 16:43
  • 1
    $\begingroup$ @Ron That depends on what you want to test and what the application is. If you have a deblurring problem, then you want to know how far a one-pixel noise in the original image $b$ spreads to the deblurred image $x$. On the other hand, if your $A$ is modeling a forward process, then you will want to know how sharp a reconstruction of a point source would be, and in that case the perturbation $\delta b$ should be $Ae_i$. $\endgroup$ – Wolfgang Bangerth Dec 21 '20 at 20:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.