Unstable linear inverse problem: which “dampening” Tikhonov matrix should I use?

A linear inverse problem is given by:

$\ \mathbf{d}=\mathbf{A}\mathbf{m}+\mathbf{e}$

where d: observed data, A: theory operator, m: unknown model and e: error.

The Least Square Error (LSE) model estimate is given by:

$\ \mathbf{\tilde{m}}=(\mathbf{A^\top A})^{-1}\mathbf{A^\top d}$

Typically the linear equation is ill-conditioned and as a consequence the model estimate tends to oscillate.

The LSE solution to the linear equation is equivalent to minimizing the following cost function:

$\ C=\|\mathbf{d-Am}\|^2$

In order to stabilize the solution we may introduce a "dampening" term:

$\ C=\|\mathbf{d-Am}\|^2 + \|\mathbf{\Gamma m}\|^2$

where $\ \Gamma$ is called the Tikhonov matrix.

The Tikhonov matrix should penalize unwanted aspect of a solution.

The stabilized LSE model estimate then becomes:

$\ \mathbf{\tilde{m}}=(\mathbf{A^\top A+\Gamma^\top \Gamma})^{-1}\mathbf{A^\top d}$

As an example problem say I have the model $\ m(x)=x^2, 0 \leq x \leq 1$.

Instead of observing the model directly I am observing the derivative of the model + some noise (n):

$\ d(x)=\frac{dm}{dx}+e$

From measurements of d(x) I want to estimate m(x). Say I know that m(0)=0.

My linear inverse problem is now:

$\ \mathbf{d}=\mathbf{A}\mathbf{m}+\mathbf{e}, \mathbf{A=D_x}$

where $\ \mathbf{D_x}$ is the differential operator (matrix) for which I use central differences.

The LSE solution to this problem oscillates a lot:

Therefore I want to dampen it by using an appropriate Tikhonov matrix.

Which Tikhonov matrix should I use?

-
Have you tried penalizing large variations in your solution by using a differentiation matrix as $\Gamma$? – begeistzwerst Jul 5 '12 at 9:57
Thanks begeistzwerst! Would this not "penalize" not flat estimations? I tried it and it seems to drag the solution towards a horizontal line. – Andy Jul 5 '12 at 11:00

Try $\Gamma = \alpha \Delta$, where $\Delta$ is a discretized version of the Laplace operator (ie, a second difference operator), and $\alpha$ is a scalar regularization parameter. The action of the Laplacean is to amplify high frequency modes, so such regularization will smooth your solution.
The laplacean's eigenfunctions are like sinusoids of higher and higher frequency, with linearly increasing eigenvalues. If you take powers of the laplacean, the eigenfunctions stay the same but the eigenvalues get raised to that power, so now the eigenvalues are increasing quadratically/cubically/whatever instead of linearly. If you want to penalize higher frequency modes more, then raising the laplacean to a power will do that. However you might also try increasing the scalar regularization parameter $\alpha$ out front instead. – Nick Alger Jul 5 '12 at 11:12