(variational-deconvolution)=
# Variational deconvolution

(classical-tikhonov-regularisation)=
## Classical Tikhonov regularisation

We have the following optimisation problem in the continuous setting:
```{math}
\min_{u \in X} \int_{\Omega} ((\Phi * u) - f)^2\, dx + \lambda \int_{\Omega} u^2\, dx, \qquad \alpha > 0.
```
Typically we choose $X = L^2(\Omega)$ and one can derive the necessary
optimality conditions via the Euler--Lagrange equations.

In the discrete setting we want to solve:
```{math}
:label: eq:variational_deconvolution_classical_tikhonov

  \min_{u \in \mathbb{R}^n} \|Au - f\|^2_2 + \lambda \|u\|^2_2.
```

(interpretation.)=
#### Interpretation.

The classical Tikhonov regularisation aims to find solutions which
balance two requirements:

1.  the solution should have a small residual $Au - f$

2.  the solution should be small in the $\ell^2$-sense.

The parameter $\lambda \geq 0$ controls this balancing. This is a
similar idea as the least-squares solution.

````{prf:lemma} Optimality conditions
:label: lem:deconvolution_classical_tikhonov_optimality
Any minimizer $u \in \R^n$ of the discrete variational deconvolution
problem {eq}`eq:variational_deconvolution_classical_tikhonov` has to
satisfy the following optimality condition (also known as normal
equations):
```{math}
u \ = \ (A^T A + \alpha I)^{-1} A^T f.
```

````

````{prf:proof} 
We begin by rewriting the variational problem
{eq}`eq:variational_deconvolution_classical_tikhonov` using the inner
product of $\R^n$.
```{math}
\begin{aligned}
E(u) &= \|Au - f\|_2^2 + \alpha \|u\|_2^2 \\
&= \langle Au - f, Au - f \rangle + \alpha \langle u, u \rangle \\
&= \langle Au, Au \rangle - \langle Au, f \rangle - \langle f, Au \rangle + \langle f, f \rangle + \alpha \langle u, u \rangle \\
&= \langle u, A^T A u \rangle - 2 \langle u, A^T f \rangle + \alpha \langle u, u \rangle + \langle f, f \rangle \\
&= \langle u, \underbrace{(A^T A + \alpha I)}_{= B} u \rangle - 2 \langle A^T f, u \rangle + \langle f, f \rangle \\
&= \langle u, Bu \rangle - 2 \langle A^T f, u \rangle + \langle f, f \rangle
\end{aligned}
```

We are interested in the directional derivative
$\langle \nabla E(u), v \rangle = \left. \frac{d}{dt} E(u + tv) \right|_{t=0}$:

```{math}
\begin{aligned}
E(u + tv) &= \langle u + tv, B(u + tv) - 2A^T f \rangle + \langle f, f \rangle \\
&= \langle u, Bu + Btv - 2A^T f \rangle + t \langle v, Bu + Btv - 2A^T f \rangle + \langle f, f \rangle \\
&= \underbrace{\langle u, Bu - 2A^T f \rangle + \langle f, f \rangle}_{= E(u)} + t \langle u, Bv \rangle + t \langle v, Bu \rangle + t^2 \langle v, Bv \rangle - t \langle v, 2A^T f \rangle \\
&\stackrel{B = B^T}{=} E(u) + 2t \langle v, Bu \rangle - 2t \langle v, A^T f \rangle + t^2 \langle v, Bv \rangle \\
&= E(u) + 2t \langle v, Bu - A^T f \rangle + t^2 \langle v, Bv \rangle
\end{aligned}
```

Thus, it follows:
```{math}
\begin{aligned}
\left. \frac{d}{dt} E(u + tv) \right|_{t=0} &= 2 \langle v, Bu - A^T f \rangle = 0 \quad \forall v \in \mathbb{R}^n \\
&\Leftrightarrow (A^T A + \alpha I)u - A^T f = 0 \\
&\Leftrightarrow u = (A^T A + \alpha I)^{-1} A^T f
\end{aligned}
```
 ◻

````

It turns out that computing a solution of the discrete variational
deconvolution problem
{eq}`eq:variational_deconvolution_classical_tikhonov` via the normal
equations in {prf:ref}`lem:deconvolution_classical_tikhonov_optimality`
is numerically too expensive as one has to invert a potentially large,
non-sparse matrix $A^TA + \alpha I$. Additionally, the condition of this
matrix is quite large and hence measurement errors in $f$ will get
amplified. For this reason we will regard the singular value
decomposition of the solution operator
$T_\lambda(f) \coloneqq (A^T A + \lambda I)^{-1} A^T f$ in the
following.

````{prf:lemma} Classical Tikhonov regularistion via SVD
Let $UDV^T = A \in \mathbb{R}^{k \times n}$ and $\lambda > 0$ fixed.
Then we can compute the Tikhonov-regularised solution of
$Au + \varepsilon = f$ as:
```{math}
u = T_\alpha(f) = VS_\alpha U^T f, \qquad \text{with } S_\lambda =
  \begin{pmatrix}
    \dfrac{\sigma_1}{\sigma_1^2 + \alpha} & & 0 \\[6pt]
    & \ddots & \\[4pt]
    0 & & \dfrac{\sigma_r}{\sigma_r^2 + \lambda}
  \end{pmatrix}
```
and $r = \min(k, n)$.

````

````{prf:proof} 
Since the SVD of the matrix $A\mathbb{R}^{k \times n}$ is given as
$A = UDV^T$ we can compute the SVD of its transposed $A^T$ via:
```{math}
A^T \: = \: (V^T)^T D^T U^T \: = \: VD^T U^T.
```

Multiplying these two identities yields:
```{math}
A^T A \: = \: (V D^T U^T)(U D V^T) \: = \: V D^T D V^T.
```

To incorporate the diagonal matrix $\lambda I$ we have to use the
identity $\underline{I} = I V V^T = V I V^T$ and thus we can write:
```{math}
A^T A + \lambda I \: = \: V D^T D V^T + V \lambda I V^T \: = \: V(D^T D + \lambda I)V^T.
```
It gets clear that the matrix $(D^T D + \lambda I)$ is a diagonal matrix
with the diagonal entries $d_{ii} = \sigma_i^2 + \lambda$, for which
$\sigma_i \geq 0$ is a singular value of $A$.

For the inversion of this operator we compute:
```{math}
(A^T A + \lambda I)^{-1} \: = \: (V(D^T D + \lambda I)V^T)^{-1} \: = \: (V^T)^{-1} (D^T D + \lambda I)^{-1} V^{-1} \: = \: V (D^T D + \lambda I)^{-1} V^T.
```
The matrix $(D^T D + \lambda I)^{-1}$ is also a diagonal matrix with the
diagonal entries $d_{ii} = \frac{1}{\sigma_i^2 + \lambda}$.

Now we can finally the SVD of the solution operator
$T_\lambda(f) = (A^T A + \lambda I)^{-1} A^T f$ as:
```{math}
\begin{aligned}
u_\lambda \: = \: V (D^T D + \lambda I)^{-1} V^T \cdot V D^T U^T f \: = \: V \underbrace{(D^T D + \lambda I)^{-1} \cdot D^T}_{=: S_\lambda} \cdot U^T f
\end{aligned}
```
The operator $S_\lambda$ is now a diagonal matrix with entries
$(S_\lambda)_{ii} = \frac{\sigma_i}{\sigma_i^2 + \lambda}$. ◻

````

**Observation:**

1.  For $\lambda = 0$ the diagonal entries of the operator $S_\lambda$
    get $(S_0)_{ii} = \frac{1}{\sigma_i}$ and hence we get the normal
    SVD.

2.  For $\lambda \to \infty$ the operator $S_\lambda$ converges to the
    zero matrix, i.e.,
    $S_\lambda \to \mathbf{0} \in \mathbb{R}^{n \times n}$.

3.  The error $T_\lambda(f)\,\varepsilon$ gets mainly amplified by
    $\dfrac{\sigma_r}{\sigma_r^2 + \lambda}$ and hence the condition of
    the operator in $T_\lambda(f)$ gets smaller for
    $\lambda \to \infty$.

4.  Since we have a quadratic problem we can give the explicit solution:
    ```{math}
    T_\lambda(f) = (A^T A + \lambda I)^{-1} A^T f.
    ```

5.  If we reinterpret the minimisation problem we see that we try to
    solve two equations simultaneously: $Au = f$ and $u = 0$. We can
    rewrite this as a linear system in a stacked form:
    ```{math}
    Bu = \begin{bmatrix} A \\ \sqrt{\lambda}\, I \end{bmatrix} u = \begin{bmatrix} f \\ 0 \end{bmatrix} =: F.
    ```
    This linear equation system can be efficiently solved.

(generalised-tikhonov-regularisation)=
## Generalised Tikhonov regularisation

We have the following optimisation problem in the continuous setting:
```{math}
\min_{u \in X} \int_{\Omega} \|(\Phi * u) - f\|^2\, dx + \alpha \int_{\Omega} \|\nabla u\|^2\, dx, \qquad \alpha > 0.
```
Typically we choose $X = H^1(\Omega) = W^{1,2}(\Omega)$ and one can
derive the necessary optimality conditions via the Euler--Lagrange
equations.

In the discrete setting we want to solve:
```{math}
\min_{u \in \mathbb{R}^n} \|Au - f\|^2 + \alpha \|Lu\|^2, \qquad \alpha > 0,
```
and the matrix $L$ is a discretisation of the differential operator,
i.e.,
```{math}
L = h \cdot
  \begin{pmatrix}
    -1 & 1 & & 0 \\
    & 0 & \ddots & \\
    & & \ddots & 1 \\
    1 & & & -1
  \end{pmatrix}
  \quad \text{for forward differences.}
```
Again we have a quadratic problem and we can compute the solution as:

1.  Normal equations (inefficient!):
    ```{math}
    T_\alpha^G(f) = (A^T A + \alpha L^T L)^{-1} A^T f.
    ```

2.  Stacked form:
    ```{math}
    Bu = \begin{bmatrix} A \\ \sqrt{\alpha}\, L \end{bmatrix} u = \begin{bmatrix} f \\ 0 \end{bmatrix} = F.
    ```

The generalised Tikhonov regularisation aims to find smooth solutions
which have a small residual $Au - f$.

(total-variation-regularisation)=
## Total Variation regularisation

We have the following optimisation problem in the continuous setting:
```{math}
\min_{u \in X} \int_{\Omega} \|(\Phi * u) - f\|^2\, dx + \alpha \int_{\Omega} |\nabla u|\, dx, \qquad \alpha > 0.
```

Compared to the general Tikhonov regularisation, this model penalises
large variations (e.g. edges) less, but small variations (e.g. noise)
stronger. As we have already seen the TV model is not differentiable in
$0$.

(primal-dual-optimization)=
### Primal-dual optimization

In the discrete setting we aim to solve the variational problem:
```{math}
\min_{u \in \mathbb{R}^n} \frac{1}{2\lambda} \|Au - f\|_2^2 + \|\nabla u\|_{1}
```

[To use]{.underline} the primal-dual optimization scheme we could set
```{math}
G(u) = \frac{1}{2\lambda} \|Au - f\|_2^2 \quad , \quad F(Ku) = \|Du\|_{1} \quad \text{mit } K := \nabla
```

[But:]{.underline} Evaluating the proximal operator
$\text{prox}_{\sigma G}(v)$ is numerically expensive
```{math}
\text{prox}_{\sigma G}(v) = \text{argmin}_{u \in \mathbb{R}^n} \frac{1}{2\sigma} \|u - v\|_2^2 + \frac{1}{2\lambda} \|Au - f\|_2^2
```
This requires solving a large, non-sparse linear system:
```{math}
\left(A^T A + \frac{\lambda}{\sigma} I\right) u = A^T f + \frac{\lambda}{\sigma} v
```

[Idea:]{.underline}\
Set $G(u) \equiv 0$ and
$F(Ku) = F\left(\begin{bmatrix} \nabla u \\ Au \end{bmatrix}\right) := F_1(\nabla u) + F_2(Au)$\
with $K := \begin{pmatrix} \nabla \\ A \end{pmatrix}$,
$F_1 := \|\cdot\|_1$, $F_2 := \|\cdot - f\|_2^2$.

Since $F$ is a sum of two functions, we need to ask ourselves what is
the Fenchel conjugate $F^*$ of $F$? For this we regard the following
simple argument. Let $v := (v_1, v_2)$ and $p := (p_1, p_2)$, then
```{math}
\begin{aligned}
F^*(p) &= \sup_{v} \langle p, v \rangle - F(v) = \sup_{v_1} \langle p_1, v_1 \rangle - F_1(v_1) + \sup_{v_2} \langle p_2, v_2 \rangle - F_2(v_2) \\
&= F_1^*(p_1) + F_2^*(p_2)
\end{aligned}
```

We already know that $F_1^*(q) = \chi_B(q)$ with
$B = \{q : \|q\|_\infty \leq 1\}$. On the other hand, we still need to
compute the Fenchel conjugate $F_2^*$ of
$F_2(v) = \frac{1}{2\lambda} \|v - f\|_2^2$. For this we regard:
```{math}
\begin{aligned}
F_2^*(q) \: = \: \sup_{v} \langle v, q \rangle - F_2(v) \: = \: \sup_{v} \langle v, q \rangle - \frac{1}{2\lambda} \|v - f\|_2^2
\end{aligned}
```
The necessary optimality conditions yield:
```{math}
\begin{aligned}
q - \frac{1}{\lambda} (\hat{v} - f) = 0 \quad &\Leftrightarrow \quad q + \frac{1}{\lambda} f = \frac{1}{\lambda} \hat{v} \\
&\Leftrightarrow \quad \hat{v} \: = \: \lambda q + f
\end{aligned}
```

Inserting this maximum $\hat{v}$ into $F_2^*(q)$ gives us:
```{math}
\begin{aligned}
F_2^*(q) \: &= \: \langle \lambda q + f, q \rangle - F_2(\lambda q + f) \: = \: \langle \lambda q + f, q \rangle - \frac{1}{2\lambda} \|\lambda q + f - f\|_2^2 \\
&= \: \lambda \|q\|_2^2 + \langle f, q \rangle - \frac{\lambda}{2} \|q\|_2^2 \: = \: \frac{\lambda}{2} \|q\|_2^2 + \langle f, q \rangle \\
& \hspace{5.5cm} \Rightarrow \partial F_2^*(q) \: = \: \lambda q + f
\end{aligned}
```

Now we investigate the proximal operators. The proximal operator
$\text{prox}_{\sigma G}(u)$ is trivial due to the choice of
$G \equiv 0$:
```{math}
\text{prox}_{\sigma G}(u) \: = \: (I + \sigma \partial G)^{-1}(u) \: = \: I^{-1}(u) \: = \: u.
```

The proximal operator $\text{prox}_{\sigma F_2^*}(p)$ is more
interesting:
```{math}
\begin{aligned}
\hat{p} = \text{prox}_{\sigma F_2^*}(p) = (I + \sigma \partial F_2^*)^{-1}(p) \quad &\Leftrightarrow \quad p = (I + \sigma \partial F_2^*)(\hat{p}) \\
&\Leftrightarrow \quad p = \hat{p} + \sigma(\lambda \hat{p} + f) \\
&\Leftrightarrow \quad p - \sigma f = (1 + \sigma \lambda) \hat{p} \\
&\Leftrightarrow \quad \hat{p} = \frac{p - \sigma f}{1 + \sigma \lambda}.
\end{aligned}
```
This point-wise operation is numerically cheap and can easily be
implemented.

Putting all necessary terms together now leads us to the following
primal-dual optimization scheme:
```{math}
\begin{aligned}
y^{k+1} &= \text{prox}_{\sigma F^*} (y^k + \sigma K \bar{x}^k) \quad \text{with } Kx := \begin{pmatrix} Dx \\ Ax \end{pmatrix} \\
x^{k+1} &= \text{prox}_{\tau G} (x^k - \tau K^* y^{k+1}) = x^k - \tau K^* y^{k+1} \quad \text{with } K^* y := -\text{div } y_1 + A^T y_2 \\
\bar{x}^{k+1} &= x^{k+1} + \theta (x^{k+1} - x^k)
\end{aligned}
```

(alternative-approach)=
### Alternative approach:

In the discrete setting we want to solve:
```{math}
\min_{u \in \mathbb{R}^n} \|Au - f\|^2 + \alpha \|Lu\|_1 \quad \Rightarrow \text{ no quadratic form!}
```

We try to reformulate our discrete problem into a quadratic form.

Let $Lu = v = v^+ - v^-$ with $v^+_j \geq 0,\, v^-_j \geq 0$ for
$j = 0, \ldots, n-1$. Then we can write:
```{math}
\min_{u,v^+,v^- \in \mathbb{R}^n} u^T A^T A u - 2f^T A u + \alpha \mathbf{1}^T v^+ + \alpha \mathbf{1}^T v^-
```

```{math}
\Rightarrow \quad \min_{y \in \mathbb{R}^{3n}} \tfrac{1}{2} x^T H x + h^T y,
```
with
```{math}
y = \begin{bmatrix} u \\ v^+ \\ v^- \end{bmatrix} \in \mathbb{R}^{3n}, \quad
  H = \begin{pmatrix} 2A^T A & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad
  h = \begin{bmatrix} -2A^T f \\ \alpha \mathbf{1} \\ \alpha \mathbf{1} \end{bmatrix}
```
under the conditions $Lu = v^+ - v^-$ and $v^+_j \geq 0,\, v^-_j \geq 0$
for $j = 0, \ldots, n$.