Vtyjkyuk

I'm working in image restoration so this is a inverse problem that usually it is ill-posed but I don't understand why inverse problems are usually ill-posed.

Direct problem: given $x$ find $y$: $Kx = y$

Inverse problem: given $y$ find $x$: $ Kx=y$

I'm working in deconvolution of images that is, I have a blurred image $y$ that it is the result of the original imagen $x$ with the convolution with a kernel $h$ so $y = h * x = H x$ where $H$ is a matrix (so we can see the convolution as a matrix product).

For obtaining the original image $x$ we have to minimize on $x$: $(y-Hx)^2$ so the minimum must have be obtained at $H^TH x = H^T y$. All that I find in the internet is that $H^T H$ is not one-to-one usually or in the best cases its eigenvalues are very small which makes the problem unstable.

Why the eigenvalues are small? why if the eigenvalues are small makes it unstable?

Let $K:Xto Y$ be a compact linear operator between two Hilbert spaces. If we have an equation of the form

$$Kx=y,$$

then we say that the problem is "ill posed" if there is not a "continuous dependence between the solution and the data". In this case, the data is normally some perturbed version of the RHS, i.e. $y_delta=y+e$, for some $deltage|e|$; and the solution is $x^dagger=K^dagger y$, where $K^dagger:Yto X$ is the Moore-Penrose pseudo-inverse. Thus, the part in bold is equivalent to saying that $K^dagger$ is continuous. Note that this condition is violated when $dimoperatornamerange K=infty$ and $overlineoperatornamerangeKne Y$

In particular, $x^dagger$ is equal to the least squares solution of minimum norm, i.e. it solves the normal equation $K^ast Kx^dagger=K^ast y$. Then you can easily compute your solution using spectral theory, i.e.

$$
beginaligned
x^dagger&=(K^star K)^-1K^ast y
\
&=int_0^^2+frac1lambda,mathrmdE_lambda K^ast y,
endaligned
$$
where $E_lambda$ denotes the spectral family of $K^ast K$. Now, you can clearly see that for small eigenvalues (i.e. $lambda$), one encounters problems. In particular, if $K$ is compact and $dim operatornamerange K=infty$, then the eigenvalues accumulate at zero!

Therefore, one needs to use regularisation. A regularisation essentially consists of replacing $1/lambda$ with a parametric spectral filter function $g_alpha(lambda)$ such that $g_alpha(lambda)to 1/lambda$ as $alphato 0$. Note that this is equivalent to replacing $K^dagger$ by a family of parametric continuous operators $R_alpha:Yto X$ called a "regularisation operator".

Regularisation

Truncated singular value decomposition:

If we set $g_alpha(lambda)=1_[alpha,infty)(lambda)lambda^-1$, then we can rewrite the integral above in the form of a regularised solution which we denote as

$$x_alpha=int_alpha^^2+frac1lambda,mathrmdE_lambda K^ast ylongrightarrow x^daggermbox textas alphato 0.$$

Notice that this corresponds to "cutting off the bad eigenvalues".

Tikhonov regularisation:

Recall the normal equation $K^ast K x=K^ast y$. We can regularise it by shifting the LHS, i.e. $(K^ast K+alpha I)x=K^ast y$. Then the solution which satisfies this equation is denoted by

$$
beginaligned
x_alpha&=(K^ast K+alpha I)^-1K^ast y
\
&=int_0^^2+frac1lambda+alpha,mathrmdE_lambda K^ast y
\
&to x^daggermbox as alphato 0.
endaligned
$$
In this case $g_alpha(lambda)=1/(lambda+alpha)$.

Of course, in practice, the $y$ is replaced by $y_delta$.

There are also many other regularisation methods which I will not list.

Parameter Choice Rules

The question still remains of how to choose the so-called regularisation parameter $alpha$? Recall that the data is usually perturbed. Thus one must choose $alpha_ast:=alpha(delta,y_delta)$ such that if $delta_n$ is a sequence tending to zero, then $x_alpha_ast,delta:= R_alpha_asty^deltato K^dagger y=:x^dagger$ as $delta_nto 0$.

A-priori parameter choice rules:

If one has knowledge of the noise level $delta$, then one can simply choose $alpha_astsimdelta$. In particular, if $x^daggerinoperatornamerange(K^astK)^mu$, for $mu>0$, which is known in the literature as a source condition (i.e. a certain condition on the smoothness of your solution), then one can even choose

$$alpha_astsimdelta^frac22mu+1,$$
which then yields optimal convergence rate:

$$|x_alpha_ast,delta-x^dagger|le Cdelta^frac2mu2mu+1,$$

for Tikhonov regularisation.

A-posteriori parameter choice rules:

In practice, one does not normally know the $mu$. An a-posteriori parameter choice rule select the parameter depending on both the noise and the noisy data. One method, which was suggested by Morozov, is the so-called discrepancy principle. That is, choose

$$alpha_ast=supKx_alpha,delta-y_delta,$$

for a constant $tau$. If the aforementioned source condition is satisfied, then this also yields optimal convergence rates for Tikhonov regularisation.

Heuristic parameter choice rules:

Suppose that you have no knowledge of the noise level (which is usually the case in practical problems), then one may opt to use a heuristic paramter choice rule (i.e. selecting a parameter which depends only on the noisy data). The drawback, due to Bakushinskii's veto, is that a heuristic regularisation method cannot converge in the worst case (i.e. if we take the supremum over all noise levels). However, it has been proven that in many situations (i.e. when the noise satisfies certain conditions, e.g. it is sufficiently irregular), then heuristic regularisation methods are indeed convergent.

A simple example is the so-called heuristic discrepancy principle:

$$alpha_ast=operatornameargmin_alphafracsqrtalpha$$

Other regularisation methods:

There is also a whole family of iterative regularisation methods, e.g. Landweber iteration which are semi-convergent in the sense that the iteration converges to an optimal solution and then diverges. Therefore, the "parameter" in this case is the stopping index and the aforementioned parameter choice rules may be slightly modified to become "stopping rules". I can expand in more detail upon request.