(deconvolution-as-inverse-problem)= # Deconvolution as inverse problem ````{prf:definition} 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: ```{math} Ku = (Ku)(x) = f(x) = f ``` based on the observation $f$ for an unknown $u$. This problem is called an *inverse problem* ```` ````{prf:definition} Hadamard conditions of inverse problems :label: def:hadamard_conditions 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: ```{math} 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. ````{prf:definition} Compact Operator :label: def: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. ````{prf:lemma} 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. ```` ````{prf:proof} We simply calculate with the Cauchy-Schwartz inequality ```{math} \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 ```{math} \|K\|_{L^2(\Omega)}\leq \|k\|_{L^2(\Omega\times\Omega)} ```  ◻ ```` ````{prf:lemma} Compactness of the integral operator The integral operator $K$ as defined above is compact. ```` ````{prf:proof} This statement can be shown via the [Arzenà-Ascoli theorem](https://en.wikipedia.org/wiki/Arzel%C3%A0%E2%80%93Ascoli_theorem). ◻ ```` Any inverse problem of the form $Ku = f$ with a compact operator $K$ is ill-posed with respect to the Hadamard conditions in {prf:ref}`def:hadamard_conditions`. 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 {prf:ref}`def:compact_operator` 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: ```{math} :label: eq:cauchy_image \|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$: ```{math} \|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 {eq}`eq:cauchy_image` yields: ```{math} :label: eq:continuity_contradiction \|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: ```{math} \|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 {eq}`eq:continuity_contradiction` yields a direct contradiction, as $\sqrt{2} < \sqrt{2}$. Consequently, the inverse operator $K^{-1}$ cannot be continuous, which violates the third Hadamard condition. (the-discrete-deconvolution-model)= ## 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 ```{math} 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: ```{math} (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: ```{math} 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 ```{math} 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}. ``` (ill-posedness-of-the-discrete-problem)= ## Ill-posedness of the discrete problem In the simple case $k = n$ one could try to perform deconvolution by inversion of $A$, i.e., ```{math} u = A^{-1} f. ``` However, if we consider the condition of $A$ with respect to the spectral norm, we observe the following: ```{math} \mathrm{cond}(A) = \frac{|\lambda_{\max}(A)|}{|\lambda_{\min}(A)|} \gg 1. ``` The following inequality gives us some idea what that tells us. ````{prf:lemma} 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: ```{math} \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: ```{math} \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}$. (i-k-n-underdetermined-system)= #### i) $k > n$ (underdetermined system) ```{math} \Rightarrow \quad \dim(\operatorname{range}(A)) \leq k ``` ```{math} \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. (ii-n-k-overdetermined-system)= #### ii) $n > k$ (overdetermined system) ```{math} \Rightarrow \quad \dim(\ker(A)) \geq 0 ``` ```{math} \Rightarrow \quad \text{we can choose an arbitrary } v_0 \in \ker(A) \text{ even for } \varepsilon = 0 ``` ```{math} \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. (iii-n-k)= #### iii) $n = k$ For our case of deconvolution, $A$ is invertible for proper convolution kernels. We have: ```{math} A^{-1} f = A^{-1}(Au + \varepsilon) = u + A^{-1}\varepsilon. ``` We can bound the error $A^{-1}\varepsilon$ by: ```{math} \|A^{-1}\varepsilon\| \leq \|A^{-1}\|\,\|\varepsilon\|. ``` Thus, continuity depends on $\|A^{-1}\|$.