6.4. Variational deconvolution

6.4.1. Classical Tikhonov regularisation

We have the following optimisation problem in the continuous setting:

\[\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:

(6.3)\[ \min_{u \in \mathbb{R}^n} \|Au - f\|^2_2 + \lambda \|u\|^2_2.\]

6.4.1.1. 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.

Lemma 6.5 (Optimality conditions)

Any minimizer \(u \in \R^n\) of the discrete variational deconvolution problem (6.3) has to satisfy the following optimality condition (also known as normal equations):

\[u \ = \ (A^T A + \alpha I)^{-1} A^T f.\]

Proof. We begin by rewriting the variational problem (6.3) using the inner product of \(\R^n\).

\[\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}\):

\[\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:

\[\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 (6.3) via the normal equations in Lemma 6.5 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.

Lemma 6.6 (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:

\[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)\).

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:

\[A^T \: = \: (V^T)^T D^T U^T \: = \: VD^T U^T.\]

Multiplying these two identities yields:

\[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:

\[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:

\[(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:

\[\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:

   \[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:

   \[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.

6.4.2. Generalised Tikhonov regularisation

We have the following optimisation problem in the continuous setting:

\[\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:

\[\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.,

\[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!):

   \[T_\alpha^G(f) = (A^T A + \alpha L^T L)^{-1} A^T f.\]

2. Stacked form:

   \[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\).

6.4.3. Total Variation regularisation

We have the following optimisation problem in the continuous setting:

\[\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\).

6.4.3.1. Primal-dual optimization

In the discrete setting we aim to solve the variational problem:

\[\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

\[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

\[\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:

\[\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

\[\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:

\[\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:

\[\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:

\[\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\):

\[\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:

\[\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:

\[\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}\]

6.4.3.2. Alternative approach:

In the discrete setting we want to solve:

\[\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:

\[\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^-\]

\[\Rightarrow \quad \min_{y \in \mathbb{R}^{3n}} \tfrac{1}{2} x^T H x + h^T y,\]

with

\[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\).