Choosing the order of Tikhonov regularization of an inverse problem

Question

This question is migrated from math.stackexchange.

Let me first describe the problem I am trying to solve and then the question I have. I greatly appreciate anyone who can shine some light on it.

There are two parameters with the relationship as follows:

$a(x) = \int_{k_1}^{k_2} b(k)\exp(-kx)dk$

where $a(x)$ is the experimental data I have (quantized in x axis), and $b(k)$ is the information I want to get.

I am trying to solve the problem with such procedure:

1. Discretize $b(k)$ in k axis, so the above integral can be formatted as matrix multiplication like $a=Mb$, where the matrix $M$ contains the information of $\exp(-kx)dk$;

2. Apply zero-order (L0) Tikhonov regularization with regularization parameter $\lambda$, and the problem can be formatted as $a=(M+\lambda I) b$ where $I$ is an identity matrix. The term $\lambda I$ has the meaning of solution smoothness. It is for damping the effect of experimental data noise on the final result $b(k)$ to avoid overfitting.

3. Use some criterion (such as L-curve criterion http://www.sintef.no/upload/IKT/9011/SimOslo/vskoler/2005/notes/Lcurve.pdf) to determine the regularization parameter $\lambda$, and solve for $b(k)$ from the problem described by Step 2.

However, I find if I change the order of Tikhonov regularization to first-order (L1) or second order (L2), meaning the smoothness term in Step 2 becomes $\lambda I_1$ or $\lambda I_2$ where $I_1$ is the first-order derivative matrix and $I_2$ is the second-order derivative matrix, I can still get solution but it is very different to that from L0.

The difference is justified because L1 and L2 put more weight on the solution smoothness rather than fitting residual. But I do not have prior information about solution $b(k)$, so I do not know which solution I should use.

Can anyone provide any guide on choosing the order of such regularization?

Answer 1

To start with: What you call Tikhonov regularization, is usually called Lavrentiev regularization (in the case of self-adjoint, non-negative definite $M$). The idea there is to shift the spectrum of $M$ away from zero. In the case of non-self-adjoint $M$, Tikhonov regularization is $$ b_\lambda = (M^*M + \lambda I)^{-1} M^* a. $$ The idea is basically the same: shift the spectrum of $M^*M$ (i.e. the singular values of $M$) away from zero. Another way to look at this is to see that $b_\lambda$ is the minimizer of $$ \|Mb-a\|^2 + \lambda\|b\|^2. $$

For choosing $\lambda$ there are several methods: Morozov's discrepancy principle (mentioned by JJ Green in the comments) amounts to choose $\lambda$ such that $$ \|Mb_\lambda - a\|\approx C\delta $$ where $\delta$ is an estimate on the noise which is present in the measurements and $C$ is a constant which is a bit larger than 1. Note that this $b_\lambda$ is also the solution of $$ \min \|b\|^2\quad\text{s.t.}\quad\|Mb - a\|\leq C\delta. $$ For this method to be applicable you need some (upper) estimate of the noise in the signal. If you think that this quantity (and also estimates on that) is not available, you need to resort to so-called heuristic parameter choice methods. They are called heuristic because you can not guarantee any error estimate for these methods (a fact known as Bakushinksii veto). Among these methods are:

• The L-curve (I suggest the work of Hansen) as you mentioned.
• The balancing principle (also called Lepskii principle - see here or here).
• The Hanke-Raus-rule (see here or (shameless plug) here).
• Generalized cross validation (see here).

Answer 2

Try "balancing principle" for parameter choice. In this method you may find an optimal solution without any knowledge of smoothness. You just need to compare two neighboring solutions with the noise level.