6.2. Deconvolution as inverse problem

Definition 6.1 (Inverse Problem)

Let \(X,Y\) be Banach spaces, let \(u \in X , f \in Y\)and let \(K: X \longrightarrow Y\) be an operator. We are interested to solve the equation:

\[Ku = (Ku)(x) = f(x) = f\]

based on the observation \(f\) for an unknown \(u\). This problem is called an inverse problem

Definition 6.2 (Hadamard conditions of inverse problems)

1. there exists a \(u \in X. Ku = f\) (existence)\

2. there exists exactly one solution (uniqueness)\

3. the solution \(u\) depends continuously on the measurement \(f\) (stability)

We observe that during our naive reconstruction via the convolution theorem that small perturbations of \(f\) lead to huge deviaions in \(u\). Thus the third Hadamard condition is violated and the deconvolution is an ill-posed, linear inverse problem (IP). We want to understand why deconvolution is an ill-posed inverse problem and investigate methods to approximate solutions.

We investigate \(Ku=f\) for which the operator \(K\) is given as an integral operator on a smooth bounded domain \(\Omega\subset\R^n\) as:

\[f(x) \ = \ (Ku)(x) \ = \ \int_{\Omega} k(x,y) \cdot u(y) dy \quad \text{for all } x\in\Omega \,.\]

For the kernel \(k(x,y) := g(x-y)\) we get the special case of a convolution. It gets clear that \(Ku=f\) is a linear inverse problem.

Definition 6.3 (Compact Operator)

Every linear map \(K: X \longrightarrow Y\) is called compact if each bounded sequence \((u_n)_{n \in \N} \subset X\) has a subsequence \((U_n)_{n \in \N} \subset X\), such that \((KU_n)_{n \in \N} \subset Y\) converges to an element \(y \in Y\).

For the moment we assume we know the blur kernel \(k(x,y)\) up to neglectable small deviations. The case for \(k\) unknown is denoted as “blind deconvolution” and will be covered later.

Lemma 6.1 (Boundedness)

Let \(K:L^2(\Omega)\to L^2(\Omega)\) be the integral operator defined above and let \(k\in L^2(\Omega\times\Omega)\). Then, \(K\) is a linear, bounded operator.

Proof. We simply calculate with the Cauchy-Schwartz inequality

\[\begin{aligned} \|Ku\|_{L^2(\Omega)}^2 &= \int_\Omega(Ku)^2(x)\, dx = \int_{\Omega}\left(\int_\Omega k(x,y)u(y)\, dy\right)^2\, dx \\ &\leq \int_{\Omega}\left(\int_\Omega k^2(x,y)\, dy\, \int_\Omega u^2(y)\, dy\right)^2\, dx = \|k\|^2_{L^2(\Omega\times\Omega)} \, \|u\|^2_{L^2(\Omega)} \end{aligned}\]

From which we follow

\[\|K\|_{L^2(\Omega)}\leq \|k\|_{L^2(\Omega\times\Omega)}\]

◻

Lemma 6.2 (Compactness of the integral operator)

The integral operator \(K\) as defined above is compact.

Proof. This statement can be shown via the Arzenà-Ascoli theorem. ◻

Any inverse problem of the form \(Ku = f\) with a compact operator \(K\) is ill-posed with respect to the Hadamard conditions in Definition 6.2. This can be seen from the following argumentation.

Let \(K: L^2(\Omega) \longrightarrow L^2(\Omega)\) be a compact operator, and let \((u_n)_{n \in \mathbb{N}} \subset L^2(\Omega)\) be an orthonormal sequence. Since \(\|u_n\|_{L^2(\Omega)} = 1\) for all \(n \in \mathbb{N}\), the sequence \((u_n)_{n \in \mathbb{N}}\) is bounded. By Definition 6.3 the image sequence \((Ku_n)_{n \in \mathbb{N}}\) contains a subsequence \((Ku_{n_k})_{k \in \mathbb{N}}\) that converges to some element \(v \in L^2(\Omega)\) in the norm topology.

Since this subsequence converges in a normed space, it is also a Cauchy sequence. Thus, for any \(\varepsilon > 0\), there exists an index \(k_0 \in \mathbb{N}\) such that for all \(i, j \geq k_0\) it holds:

(6.1)\[\|Ku_{n_i} - Ku_{n_j}\|_{L^2(\Omega)} \, < \, \varepsilon.\]

Now, suppose that a bounded (continuous) inverse operator \(K^{-1}: \mathcal{R}(K) \longrightarrow L^2(\Omega)\) of \(K\) exists. Then, there must exist a constant \(C > 0\) such that for all \(i, j \geq k_0\):

\[\|u_{n_i} - u_{n_j}\|_{L^2(\Omega)} \ = \ \|K^{-1}(Ku_{n_i} - Ku_{n_j})\|_{L^2(\Omega)} \ \leq \ C \|Ku_{n_i} - Ku_{n_j}\|_{L^2(\Omega)}\]

Applying the estimate from (6.1) yields:

(6.2)\[\|u_{n_i} - u_{n_j}\|_{L^2(\Omega)} < C \varepsilon\]

However, since \((u_n)_{n \in \mathbb{N}}\) is an orthonormal sequence, the distance between any two distinct elements for \(i \neq j\) is strictly bounded away from zero by the Pythagorean theorem:

\[\|u_{n_i} - u_{n_j}\|_{L^2(\Omega)} \ = \ \sqrt{\|u_{n_i}\|^2_{L^2(\Omega)} -2 \langle u_{n_i}, u_{n_j} \rangle + \|u_{n_j}\|^2_{L^2(\Omega)}} = \sqrt{1 - 0 + 1} = \sqrt{2}\]

Choosing \(\varepsilon \le \frac{\sqrt{2}}{C}\) in (6.2) yields a direct contradiction, as \(\sqrt{2} < \sqrt{2}\). Consequently, the inverse operator \(K^{-1}\) cannot be continuous, which violates the third Hadamard condition.

6.2.1. The discrete deconvolution model

In the following we will coarsely follow the argumentation in:

[1] Jennifer L. Müller, Samuli Siltanen: Linear and Nonlinear Inverse Problems with Practical Applications, SIAM (2012).

For the sake of simplicity we will restrict ourselves to the case \(n = 1\). Let \(u \colon [0,1] \to \mathbb{R}\) and \(\Phi \colon [-a,a] \to \mathbb{R}\) with compact support in \((-a,a)\) and \(0 < a < \tfrac{1}{2}\). We model our given data as discrete measurements \(f \coloneqq \bigl[(\Phi * u)(x_0),\, \ldots,\, (\Phi * u)(x_{k-1})\bigr] \in \mathbb{R}^k.\) We choose \(r \coloneqq \min\{k > 0 : (k+1) \cdot h > a\}\) for the step size \(h = \tfrac{1}{n}\). Then we get a discretisation of the point spread function \(\Phi\) on \(\Omega_h = \{x_0, \ldots, x_{n-1}\}\) by

\[p(x_i) = \frac{\Phi(x_i)}{h \cdot \displaystyle\sum_{j=-r}^{r} \Phi(x_j)}, \qquad i = 0, \ldots, n-1.\]

By this we know that \(p\) is normalised, i.e. \(h \cdot \sum_{j=-r}^{r} p(x_j) = 1.\) We now get a discretisation of the (continuous) convolution as:

\[(p * u)(x_i) = \sum_{j=-r}^{r} p(x_j) \cdot u(x_{i-j \bmod n}), \qquad i = 0, \ldots, n-1.\]

It gets clear that \(h \cdot (p * u) \approx (\Phi * u)\) on \([0,1]\) and thus:

\[f(x_i) = h \cdot (p * u)(x_i), \qquad i = 0, \ldots, k-1.\]

We know that deconvolution is a linear inverse problem and hence we can write our discrete problem as \(f=Au\) with

\[A = h \cdot \begin{pmatrix} p_0 & p_{-1} & \cdots & p_{-r} & 0 & \cdots & 0 & p_r & \cdots & p_1 \\ p_1 & p_0 & p_{-1} & \cdots & p_{-r} & \cdots & 0 & p_r & \cdots & p_2 \\ \vdots & & & & & & & & & \vdots \\ p_r & & & & & & & & & \vdots \\ 0 & & & & & & & & & \vdots \\ \vdots & & & & & & & & & \vdots \\ p_{-1} & \cdots & p_{-r} & 0 & \cdots & 0 & p_r & \cdots & p_1 & p_0 \end{pmatrix} \in \mathbb{R}^{k \times n}.\]

6.2.2. Ill-posedness of the discrete problem

In the simple case \(k = n\) one could try to perform deconvolution by inversion of \(A\), i.e.,

\[u = A^{-1} f.\]

However, if we consider the condition of \(A\) with respect to the spectral norm, we observe the following:

\[\mathrm{cond}(A) = \frac{|\lambda_{\max}(A)|}{|\lambda_{\min}(A)|} \gg 1.\]

The following inequality gives us some idea what that tells us.

Lemma 6.3 (Condition number bound)

Let \(x\) be the solution of \(Ax = b\) and let \(\tilde{x}\) be the solution of \(A\tilde{x} = \tilde{b} = b + \varepsilon\). Then the following holds:

\[\frac{\|\tilde{x} - x\|}{\|x\|} \leq \mathrm{cond}(A)\, \frac{\|\tilde{b} - b\|}{\|b\|}.\]

For the case \(k = n\) the convolution operator \(A\) is a square matrix with rapidly decaying eigenvalues:

\[\begin{aligned} &\Rightarrow \quad \mathrm{cond}(A) \gg 1 \\ &\Rightarrow \quad \text{any inversion of } A \text{ is highly sensitive to noise} \\ &\Rightarrow \quad \text{deconvolution is an ill-posed problem.} \end{aligned}\]

We introduce regularisation to perform robust deconvolution. We try to stabilise by changing the singular values accordingly.

First, we will investigate the violation of Hadamard’s well-posedness conditions for the noisy model \(Au + \varepsilon = f\), e.g. for \(\varepsilon \sim \mathcal{N}(0,\varepsilon)\), based on the shape of the convolution matrix \(A \in \mathbb{R}^{k \times n}\).

6.2.2.1. i) \(k > n\) (underdetermined system)

\[\Rightarrow \quad \dim(\operatorname{range}(A)) \leq k\]

\[\Rightarrow \quad \text{even for } \varepsilon = 0 \text{ there exists an } m_0 \in \mathbb{R}^k \text{ such that no } v \in \mathbb{R}^n \text{ exists with } Av = m_0.\]

\(\Rightarrow\) no solutions exist.

6.2.2.2. ii) \(n > k\) (overdetermined system)

\[\Rightarrow \quad \dim(\ker(A)) \geq 0\]

\[\Rightarrow \quad \text{we can choose an arbitrary } v_0 \in \ker(A) \text{ even for } \varepsilon = 0\]

\[\text{such that } A\bigl(\underbrace{A^{-}f}_{=0}\bigr) = f = A\bigl(\underbrace{A^{-}f + v_0}_{=:u}\bigr)\]

\(\Rightarrow\) no uniqueness of solutions.

6.2.2.3. iii) \(n = k\)

For our case of deconvolution, \(A\) is invertible for proper convolution kernels. We have:

\[A^{-1} f = A^{-1}(Au + \varepsilon) = u + A^{-1}\varepsilon.\]

We can bound the error \(A^{-1}\varepsilon\) by:

\[\|A^{-1}\varepsilon\| \leq \|A^{-1}\|\,\|\varepsilon\|.\]

Thus, continuity depends on \(\|A^{-1}\|\).