---
abstract: 'Diverse inverse problems in imaging can be cast as variational problems composed of a task-specific data fidelity term and a regularization term. In this paper, we propose a novel learnable general-purpose regularizer exploiting recent architectural design patterns from deep learning. We cast the learning problem as a discrete sampled optimal control problem, for which we derive the adjoint state equations and an optimality condition. By exploiting the variational structure of our approach, we perform a sensitivity analysis with respect to the learned parameters obtained from different training datasets. Moreover, we carry out a nonlinear eigenmode analysis, which reveals interesting properties of the learned regularizer. We show state-of-the-art performance for classical image restoration and medical image reconstruction problems.'
author:
- |
    Erich Kobler\
    Graz University of Technology\
    [erich.kobler@icg.tugraz.at]{}
- |
    Alexander Effland\
    Graz University of Technology\
    [effland@tugraz.at]{}
- |
    Karl Kunisch\
    University of Graz\
    [karl.kunisch@uni-graz.at]{}
- |
    Thomas Pock\
    Graz University of Technology\
    [pock@icg.tugraz.at]{}
bibliography:
- 'references.bib'
title: Total Deep Variation for Linear Inverse Problems
---

Introduction
============

The statistical interpretation of linear inverse problems allows the treatment of measurement uncertainties in the input data $z$ and missing information in a rigorous framework. Bayes’ theorem states that the posterior distribution $p(x\vert z)$ is proportional to the product of the data likelihood $p(z\vert x)$ and the prior $p(x)$, which represents the belief in a certain solution $x$ given the input data $z$. A classical estimator for $x$ is given by the maximum a posterior (MAP) estimator, which in a negative log-domain amounts to minimizing the variational problem $$\mathcal{E}(x)\coloneqq\mathcal{D}(x,z)+\mathcal{R}(x).$$ Here, the data fidelity term $\mathcal{D}$ can be identified with the negative log-likelihood $-\log p(z\vert x)$ and the regularization term corresponds to the negative log-probability of the prior distribution$-\log p(x)$. Assuming Gaussian noise in the data $z$, the data fidelity term naturally arises from the negative log-Gaussian, which essentially leads to a quadratic $\ell^2$-term of the form $\mathcal{D}(x,z)\coloneqq\frac{1}{2}\Vert Ax-z\Vert_2^2$. In this paper, we assume that $A$ is a task-specific linear operator (see [@ChPo16]) such as a downsampling operator in the case of single image super-resolution. While the data fidelity term is straightforward to model, the grand challenge in inverse problems for imaging is the design of a regularizer that captures the complexity of the statistics of natural images.

![ Visualization of the proposed deep regularizer ${r}$ (TDV$^3$), where colors indicate different hierarchical levels. On the highest level, ${r}(x,\theta)=w^\top \mathcal{N}(Kx)$ assigns to each pixel an energy value incorporating the local neighborhood. The function $\mathcal{N}$ (yellow) is composed of three macro-blocks (red), each representing a CNN with a U-Net type architecture. The macro-blocks consist of five micro-blocks (blue) with a residual structure on three scales. []{data-label="fig:network"}](DVN.pdf){width=".4\linewidth"}

A classical and widely used regularizer is the total variation (TV) originally proposed in [@RuOsFa92], which is based on the first principle assumption that images are piecewise constant with sparse gradients. A well-known caveat of the sparsity assumption of TV is the formation of clearly visible artifacts known as staircasing effect. To overcome this problem, the first principle assumption has later been extended to piecewise smooth images incorporating higher order image derivatives such as infimal convolution based models [@ChLi97] or the total generalized variation [@BrKu10]. Inspired by the fact that edge continuity plays a fundamental role in the human visual perception, regularizers penalizing the curvature of level lines have been proposed in [@NiMu93; @ChKa02; @ChPo19]. While these regularizers are mathematically well-understood, the complexity of natural images is only partially reflected in their formulation. For this reason, handcrafted variational methods have nowadays been largely outperformed by purely data-driven methods.

It has been recognized quite early that a proper statistical modeling of regularizers should be based on learning [@ZhWu98]. One of the most successful early approaches is the Fields of Experts (FoE) regularizer [@RoBl09], which can be interpreted as a generalization of the total variation, but builds upon learned filters and learned potential functions. While the FoE prior was originally learned generatively, it was shown in [@SaTa09] that a discriminative learning via implicit differentiation yields improved performance. A computationally more feasible method for discriminative learning is based on unrolling a finite number of iterations of a gradient descent algorithm [@Do12]. Additionally using iteration dependent parameters in the regularizer was shown to significantly increase the performance (TNRD [@ChPo17], [@Le16]). In [@KoKl17], variational networks (VNs) are proposed, which give an incremental proximal gradient interpretation of TNRD.

Interestingly, such truncated schemes are not only computationally much more efficient, but are also superior in performance with respect to the full minimization. A continuous time formulation of this phenomenon was proposed in [@EfKo19] by means of an optimal control problem, within which an optimal stopping time is learned.

An alternative approach to incorporate a regularizer into a proximal algorithm, known as plug-and-play prior [@VeSi13] or regularization by denoising [@RoEl17], is the replacement of the proximal operator by an existing denoising algorithm such as BM3D [@DaFo07]. Combining this idea with deep learning was proposed in [@MeMo17]. However, all the aforementioned schemes lack a variational structure and thus are not interpretable in the framework of MAP inference.

In this paper, we introduce a novel regularizer, which is inspired by the design patterns of state-of-the-art deep convolutional neural networks and simultaneously ensures a variational structure. We achieve this by representing the total energy of the regularizer by means of a residual multi-scale network (Figure \[fig:network\]) leveraging smooth activation functions. In analogy to [@EfKo19], we start from a gradient flow of the variational energy and utilize a semi-implicit time discretization, for which we derive the discrete adjoint state equation using the discrete Pontryagin maximum principle. Furthermore, we present a first order necessary condition of optimality to automatize the computation of the optimal stopping time. Our proposed *Total Deep Variation* (TDV) regularizer can be used as a generic regularizer in variational formulations of linear inverse problems. The major contributions of this paper are as follows:

-   The design of a novel generic multi-scale variational regularizer learned from data.

-   A rigorous mathematical analysis including a sampled optimal control formulation of the learning problem and a sensitivity analysis of the learned estimator with respect to the training dataset.

-   A nonlinear eigenmode analysis for the visualization and understanding of the learned regularizer.

-   State-of-the-art results on a number of classical image restoration and medical image reconstruction problems with an impressively low number of learned parameters.

Due to space limitations, all proofs are presented in the supplementary material.

Sampled optimal control problem
===============================

Let $x\in{\mathbb{R}}^{nC}$ be a corrupted input image with a resolution of $n=n_1\cdot n_2$ and $C$ channels. We emphasize that the subsequent analysis can easily be transferred to image data in any dimension.

The variational approach for inverse problems frequently amounts to computing a minimizer of a specific energy functional $\mathcal{E}$ of the form $\mathcal{E}(x,\theta,z)\coloneqq\mathcal{D}(x,z)+\mathcal{R}(x,\theta)$. Here, $\mathcal{D}$ is a data fidelity term and a $\mathcal{R}$ is a parametric regularizer depending on the learned training parameters $\theta\in\Theta$, where $\Theta\subset{\mathbb{R}}^p$ is a compact and convex set of training parameters. Typically, the minimizer of this variational problem is considered as an approximation of the uncorrupted ground truth image. Accordingly, in this paper we analyze different data fidelity terms of the form $\mathcal{D}(x,z)=\frac{1}{2}\Vert Ax-z\Vert_2^2$ for fixed task-dependent $A\in{\mathbb{R}}^{lC\times nC}$ and fixed $z\in{\mathbb{R}}^{lC}$.

The proposed total deep variation (TDV) regularizer is given by $\mathcal{R}(x,\theta)=\sum_{i=1}^n r(x,\theta)_i\in{\mathbb{R}}$, which is the total sum of the pixelwise deep variation defined as ${r}(x,\theta)=w^\top\mathcal{N}(Kx)\in{\mathbb{R}}^n$. Here, $K\in{\mathbb{R}}^{nm\times nC}$ is a learned convolution kernel with zero-mean constraint (i.e. $\sum_{i=1}^{nC}K_{j,i}=0$ for $j=1,\ldots,nm$), $\mathcal{N}:{\mathbb{R}}^{nm}\to{\mathbb{R}}^{nq}$ is a multiscale convolutional neural network and $w\in{\mathbb{R}}^q$ is a learned weight vector. Hence, $\theta$ encodes the kernel $K$, the weights of the convolutional layers in $\mathcal{N}$ and the weight vector $w$. The exact form of the network is depicted in Figure \[fig:network\]. In the TDV$^l$ network for $l\in{\mathbb{N}}$, the function $f$ is composed of $l$ consecutive macro-blocks $\mathbf{Ma}^1,\ldots,\mathbf{Ma}^l$ (red blocks). Each macro-block $\mathbf{Ma}^i$ ($i\in\{1,\ldots,l\}$) has a U-Net type architecture [@RoFi15] with five micro-blocks $\mathbf{Mi}_1^i,\ldots,\mathbf{Mi}_5^i$ (blue blocks) distributed over three scales with skip-connections on each scale. In addition, residual connections are added between different scales of consecutive macro-blocks whenever possible. Finally, each micro-block $\mathbf{Mi}_j^i$ for $i\in\{1,\ldots,l\}$ and $j\in\{1,\ldots,5\}$ has a residual structure defined as $\mathbf{Mi}_j^i(x)=x+K_{j,2}^i\phi(K_{j,1}^i x)$ for convolution operators $K_{j,1}^i$ and $K_{j,2}^i$. We use a smooth log-student-t-distribution of the form $\phi(x)=\frac{1}{2\nu}\log(1+\nu x^2)$, which is an established model for the statistics of natural images [@HuMu99] and has the properties $\phi'(0)=0$ and $\phi''(0)=1$. All convolution operators $K_{j,1}^i$ and $K_{j,2}^i$ in each micro-block correspond to $3\times3$ convolutions with $m$ feature channels and no additional bias. To avoid aliasing, downsampling and upsampling are implemented incorporating $3\times 3$ convolutions and transposed convolutions with stride $2$ using a blurring of the kernels proposed by [@Zh19]. Finally, the concatenation maps $2m$ feature channels to $m$ feature channels using a single $1\times 1$ convolution.

A common strategy for the minimization of the energy $\mathcal{E}$ is the incorporation of a *gradient flow* [@AmGi08] posed on a finite time interval $(0,T)$, which reads as $$\begin{aligned}
\dot{\tilde{x}}(t)=&-\nabla_1\mathcal{E}(\tilde{x}(t),\theta,z)=f(\tilde{x}(t),\theta,z)\\
\coloneqq&-A^\top(A\tilde{x}(t)-z)-\nabla_1\mathcal{R}(\tilde{x}(t),\theta)\end{aligned}$$ for $t\in(0,T)$, where $\tilde{x}(0)=x_{\mathrm{init}}$ for a fixed *initial value* $x_{\mathrm{init}}\in{\mathbb{R}}^{nC}$. Here, $\tilde{x}:[0,T]\to{\mathbb{R}}^{nC}$ denotes the differentiable *flow of $\mathcal{E}$* and $T>0$ refers to the *time horizon*. We assume that $T\in[0,{T_\mathrm{max}}]$ for a fixed ${T_\mathrm{max}}>0$. To formulate optimal stopping, we exploit the reparametrization $x(t)=\tilde{x}(tT)$, which results for $t\in(0,1)$ in the equivalent gradient flow formulation $$\dot{x}(t)=Tf(x(t),\theta,z),\qquad x(0)=x_{\mathrm{init}}.$$

Following [@EHa19; @LiCh17; @LiHa18], we cast the training process as a *sampled optimal control problem* with control parameters $\theta$ and $T$. For fixed $N\in{\mathbb{N}}$, let $(x_{\mathrm{init}}^i,y^i,z^i)_{i=1}^N\in({\mathbb{R}}^{nC}\times{\mathbb{R}}^{nC}\times{\mathbb{R}}^{lC})^N$ be a collection of $N$ triplets of initial-target image pairs, and observed data independently drawn from a task-dependent fixed probability distribution. In additive Gaussian image denoising, for instance, a ground truth image $y^i$ is deteriorated by noise $n^i\sim\mathcal{N}(0,\sigma^2)$ and a common choice is $x_{\mathrm{init}}^i=z^i=y^i+n^i$. Let ${l}:{\mathbb{R}}^{nC}\to{\mathbb{R}}_0^+$ be a convex, twice continuously differentiable and coercive (i.e. $\lim_{\Vert x\Vert_2\to\infty}{l}(x)=+\infty$) function. We will later use ${l}(x)=\sqrt{\Vert x\Vert_1^2+\varepsilon^2}$ for $\varepsilon>0$ and ${l}(x)=\frac{1}{2}\Vert x\Vert_2^2$. Then, the sampled optimal control problem reads as $$\inf_{T\in[0,{T_\mathrm{max}}],\,\theta\in\Theta}
\left\{J(T,\theta)\coloneqq\frac{1}{N}\sum_{i=1}^N{l}(x^i(1)-y^i)\right\}
\label{eq:objectiveFunction}$$ subject to the state equation for each sample $$\dot{x}^i(t)=Tf(x^i(t),\theta,z^i),\qquad x^i(0)=x_{\mathrm{init}}^i\label{eq:gradientFlowSample}$$ for $i=1,\ldots,N$ and $t\in(0,1)$. The next theorem ensures the existence of solutions to this optimal control problem.

The minimum in  subject to the side conditions  is attained.

Discretized optimal control problem
===================================

In this section, we present a novel fully discrete formulation of the previously introduced sampled optimal control problem. The state equation  is discretized using a semi-implicit scheme resulting in $$x_{s+1}^i=x_s^i-\tfrac{T}{S}A^\top(Ax_{s+1}^i-z^i)-\tfrac{T}{S}\nabla_1\mathcal{R}(x_s^i,\theta)\in{\mathbb{R}}^{nC}
\label{eq:stateEquationDiscrete}$$ for $s=0,\ldots,S-1$ and $i=1,\ldots,N$, where the *depth* $S\in{\mathbb{N}}$ is a priori fixed. This equation is equivalent to $x_{s+1}^i=\widetilde{f}(x_s^i,T,\theta,z^i)$ with $$\widetilde{f}(x,T,\theta,z)\coloneqq({\mathrm{Id}}+\tfrac{T}{S}A^\top A)^{-1}(x+\tfrac{T}{S}(A^\top z-\nabla_1\mathcal{R}(x,\theta))).
\label{eq:discreteRightHandSide}$$ The initial state satisfies $x_0^i=x_{\mathrm{init}}^i\in{\mathbb{R}}^{nC}$. Then, the discretized sampled optimal control problem reads as $$\inf_{T\in[0,{T_\mathrm{max}}],\,\theta\in\Theta}
\left\{J_S(T,\theta)\coloneqq\frac{1}{N}\sum_{i=1}^N{l}(x_S^i-y^i)\right\}
\label{eq:discreteOptimalControl}$$ subject to $x_{s+1}^i=\widetilde{f}(x_s^i,T,\theta,z^i)$. Following the discrete Pontryagin maximum principle [@Ha66; @LiHa18], the associated discrete adjoint state $p_s^i$ is given by $$p_s^i=({\mathrm{Id}}-\tfrac{T}{S}\nabla_1^2\mathcal{R}(x_s^i,\theta))({\mathrm{Id}}+\tfrac{T}{S}A^\top A)^{-1}p_{s+1}^i
\label{eq:adjointEquationDiscrete}$$ for $s=S-1,\ldots,0$ and $i=1,\ldots,N$, and the terminal condition $p_S^i=-\frac{1}{N}\nabla{l}(x_S^i-y^i)$. For further details, we refer the reader to the supplementary material.

The next theorem states an exactly computable condition for the optimal stopping time, which is of vital importance for the numerical optimization:

Let $(\overline{T},\overline{\theta})$ be a stationary point of $J_S$ with associated states $\overline{x}_s^i$ and adjoint states $\overline{p}_s^i$ satisfying  and subject to the initial conditions $\overline{x}_0^i=x_{\mathrm{init}}^i$ and the terminal conditions $\overline{p}_S^i=-\frac{1}{N}\nabla{l}(\overline{x}_S^i-y^i)$. We further assume that $\nabla \widetilde{f}(\overline{x}_s^i,\overline{T},\overline{\theta},z^i)$ has full rank for all $i=1,\ldots,N$ and $s=0,\ldots,S$. Then, $$-\frac{1}{N}\sum_{s=0}^{S-1}\sum_{i=1}^N\langle\overline{p}_{s+1}^i,({\mathrm{Id}}+\tfrac{\overline{T}}{S}A^\top A)^{-1}(\overline{x}_{s+1}^i-\overline{x}_s^i)\rangle=0.
\label{eq:optimalTDiscrete}$$

Note that  is the derivative of the discrete Lagrange functional minimizing $J_S$ subject to the discrete state equation. The proof and a time continuous version of this theorem are presented in the supplementary material.

An important property of any learning based method is the dependency of the learned parameters with respect to different training datasets [@GeBi92].

\[thm:sensitivity\] Let $(T,\theta),(\widetilde{T},\widetilde{\theta})$ be two pairs of control parameters obtained from two different training datasets. We denote by $x,\widetilde x\in({\mathbb{R}}^{nC})^{(S+1)}$ two solutions of the state equation with the same observed data $z$ and initial condition $x_{\mathrm{init}}$, i.e. $$x_{s+1}=\widetilde{f}(x_s,T,\theta,z),\quad
\widetilde{x}_{s+1}=\widetilde{f}(\widetilde{x}_s,\widetilde{T},\widetilde{\theta},z)$$ for $s=1,\ldots,S-1$ and $x_0=\widetilde{x}_0=x_{\mathrm{init}}$. Let $B(T)\coloneqq{\mathrm{Id}}+\frac{T}{S}A^\top A$ and  $L_\mathcal{R}$ be the Lipschitz constant of $\mathcal{R}$, i.e. $$\Vert\nabla_1\mathcal{R}(x,\theta)-\nabla_1\mathcal{R}(\widetilde{x},\widetilde{\theta})\Vert_2\leq L_\mathcal{R}
\left\Vert
\begin{pmatrix}
x\\
\theta
\end{pmatrix}
-
\begin{pmatrix}
\widetilde{x}\\
\widetilde{\theta}
\end{pmatrix}
\right\Vert_2
\label{eq:Lipschitz}$$ for all $x,\widetilde{x}\in{\mathbb{R}}^{nC}$ and all $\theta,\widetilde{\theta}\in\Theta$. Then, $$\begin{aligned}
\Vert x_{s+1}-\widetilde{x}_{s+1}\Vert_2
\leq&\Vert B(T)^{-1}-B(\widetilde{T})^{-1}\Vert_2
\Big(\Vert x_s\Vert_2+\tfrac{T}{S}\Vert A^\top z\Vert_2
+\tfrac{T}{S}\Vert\nabla_1\mathcal{R}(x_s,\theta)\Vert_2\Big)\notag\\
&+\Vert B(\widetilde{T})^{-1}\Vert_2
\Big(\Vert x_s-\widetilde{x}_s\Vert_2
+\tfrac{\vert T-\widetilde{T}\vert}{S}\Vert A^\top z\Vert_2+\tfrac{\vert T-\widetilde{T}\vert}{S}\Vert\nabla_1\mathcal{R}(\widetilde{x}_s,\widetilde{\theta})\Vert_2\notag\\
&+\tfrac{T}{S}L_\mathcal{R}\Vert(x,\theta)^\top-(\widetilde{x},\widetilde{\theta})^\top\Vert_2\Big).
\label{eq:sensitivity}\end{aligned}$$

Hence, this theorem provides a computable upper bound for the norm difference of two states with observed value $z$ and initial value $x_{\mathrm{init}}$ evaluated at the same step $s$, which amounts to a sensitivity analysis w.r.t. the training data.

Numerical results
=================

In this section, we elaborate on the training and optimization, and we analyze the numerical results for four exemplary applications of TDVs: image denoising, CT and MRI reconstruction, and single image super-reconstruction.

Training and optimization
-------------------------

For all models considered, we solely use 400 images taken from the BSDS400 dataset [@MaFo01] for training, where we apply a data augmentation by flipping and rotating the images by multiples of $90^\circ$. We compute the approximate minimizers of the discretized sampled optimal control problem  using the stochastic ADAM optimizer [@KiBa15] on batches of size $32$ with a patch size $96\times 96$. The zero-mean constraint of the kernel $K$ is enforced by a projection after each iteration step. For image denoising, we incorporate the squared $\ell^2$-loss function ${l}(x)=\frac{1}{2}\Vert x\Vert_2^2$ and the learning rate $4\cdot 10^{-4}$, whereas for single image super-resolution the regularized $\ell^1$-loss ${l}(x)=\sqrt{\Vert x\Vert_1^2+\varepsilon^2}$ with $\varepsilon=10^{-3}$ and the learning rate $10^{-3}$ is used. In all cases, the first and second order momentum variables of the ADAM optimizer are set to $0.9$ and $0.999$ and $10^6$ training steps are performed. Throughout all experiments we set $\nu=9$, the number of feature channels $m=32$, and $S=10$ if not otherwise stated.

Image denoising
---------------

In the first task, we analyze the performance of the TDV regularizer for additive white Gaussian denoising. To this end, we create the training dataset by uniformly drawing a ground truth patch $y^i$ from the BSDS400 dataset and add Gaussian noise $n^i\sim\mathcal{N}(0,\sigma^2)$. Consequently, we set $x_{\mathrm{init}}^i=z^i=y^i+n^i$ and the linear operator $A$ coincides with the identity matrix ${\mathrm{Id}}$.

Table \[tab:awgn\] lists the average PSNR values of classical image test datasets for varying noise levels $\sigma\in\{15,25,50\}$. In the penultimate column, the PSNR values of our proposed TDV regularizer with three macro-blocks solely trained for $\sigma=25$ (denoted by TDV$^3_{25}$) are presented. To apply the TDV$^3_{25}$ model to different noise levels, we first rescale the noisy images $\widehat{x}_{\mathrm{init}}^i=\widehat{z}^i=\tfrac{25}{\sigma}z^i$, then apply the learned scheme , and obtain the results via $x_S^i=\tfrac{\sigma}{25}\widehat{x}_S^i$. In the last column, the PSNR values of the proposed TDV regularizer with three macro-blocks (denoted by TDV$^3$) *individually* trained for each specific noise levels are shown. For all considered noise levels and datasets, state-of-the-art FOCNet [@JiLi19] achieves the best results in terms of PSNR score. The proposed TDV$^3_{25}$ and TDV$^3$ models have comparable PSNR values. The adaption of the TDV$^3$ model to specific noise levels further increases the PSNR value. To conclude, we emphasize that we achieve results on par with FOCNet with less than 1% of the trainable parameters, which highlights the potential of our approach. Compared to FOCNet, the *inherent variational structure* of our model allows for a *deeper mathematical analysis*, that we elaborate on below.

Figure \[fig:energyLandscape\] depicts surface plots of the deep variation $[-1,1]\ni(\xi_1,\xi_2)\mapsto{r}(\xi_1 x+\xi_2 n)_i$ of TDV$_{25}^3$ evaluated at four prototypic patches $x$ of size $49\times 49$, where $i$ is the index of the center pixel marked by the red points. Here, $n$ refers to a Gaussian noise with standard deviation $\sigma=25$. Hence, the surface plots visualize the local regularization energy in the image contrast direction and a random noise direction. The deep variation is smooth with commonly only a single local minimizer and the shape significantly varies depending on the pixel neighborhood.

![Surface plots of the deep variation $[-1,1]\ni(\xi_1,\xi_2)\mapsto{r}(\xi_1 x+\xi_2 n)_i$ of four patches–each evaluated at the red center pixel using TDV$^3$ with BSDS68 and $\sigma=25$.[]{data-label="fig:energyLandscape"}](figure0.pdf){width=".5\linewidth"}

Figure \[fig:PSNRDenoising\] (top) visualizes the average PSNR values $S\mapsto\frac{1}{N}\sum_{i=1}^N\mathrm{PSNR}(x_S^i,y^i)$ on the BSDS68 test dataset for $\sigma=25$. In addition, the second plot $S\mapsto\overline{T}$ (bottom) shows the learned optimal stopping time as a function of the depth. In all experiments, TDV regularizers with more macro-blocks perform significantly better. Moreover, the average PSNR value saturates beyond the depth $S=10$. The optimal stopping time converges for large $S$ in all models. Consequently, larger depth values $S$ lead to a finer time discretization of the trajectories and depth values $S$ exceeding 10 yield no further improvement.

![Plots of the functions $S\mapsto\frac{1}{N}\sum_{i=1}^N\mathrm{PSNR}(x_S^i,y^i)$ (top) and $S\mapsto\overline{T}$ (bottom) for TDV$^1$ (blue), TDV$^2$ (orange) and TDV$^3$ (green).[]{data-label="fig:PSNRDenoising"}](S_PSNR_T.pdf){width=".5\linewidth"}

To address the significance of the optimal stopping time $\overline{T}$, we evaluate the PSNR values (top) and the first order condition  (bottom) as a function of the stopping time in Figure \[fig:OptimalStoppingTime\] for TDV$_{25}^3$ trained for $S=10$. Each black dashed curve represents a single image of the BSDS68 test dataset and the red curve is the average among all test samples. Initially, all PSNR curves monotonically increase up to a unique maximum value located near the optimal stopping time $\overline{T}=0.0297$ and decrease for larger values of $T$, which is exactly determined by the first order condition . We stress that all curves peak around a very small neighborhood, which results from an overlap of all curves near the zero crossing in the second plot.

![Plots of the functions $T\mapsto\mathrm{PSNR}(x_S^i,y^i)$ (top) and the first order condition $T\mapsto-\sum_{s=0}^S \langle p_{s+1}^i,x_{s+1}^i-x_s^i\rangle$ (bottom) for $i=1,\ldots,68$ using TDV$^3$ for BSDS68 and $\sigma=25$. The averages across the samples are depicted by the red curves.[]{data-label="fig:OptimalStoppingTime"}](FOC_plot.pdf){width=".5\linewidth"}

A qualitative and quantitative analysis of the impact of the stopping time for TDV$_{25}^3$ trained with $S=10$ on a standard test image is depicted in Figure \[fig:OptimalStoppingSequence\]. Starting from the noisy input image, the restored image sequence for increasing $S\in\{5,10,15,20\}$ gradually removes noise. Beyond the optimal stopping time, the model smoothes out fine details. We emphasize that even high-frequency patterns like the vertical stripes on the chimney are consistently preserved. Thus, the proposed model generates a stable and interpretable transition from noisy images to cartoon-like images and the optimal stopping time determines the best intermediate result in terms of PSNR.

![From left to right: Ground truth, noisy input with noise level $\sigma=25$ and resulting output of TDV$^3$ for $(S,T)\in\{(5,\frac{\overline{T}}{2}),(10,\overline{T}),(15,\frac{3\overline{T}}{2}),(20,2\overline{T})\}$, where the optimal stopping time is $\overline{T}=0.0297$. Note that the best image is framed in red.[]{data-label="fig:OptimalStoppingSequence"}](figure1.pdf){width="\linewidth"}

To analyze the local behavior of TDV$_{25}^3$, we examine the variational problem $$(\overline{x},\overline{\lambda})\in\operatorname*{argmax}_{x\in{\mathbb{R}}^{nC},\lambda\in{\mathbb{R}}}-\mathcal{R}(x,\theta)+\frac{\lambda}{2}(\Vert x\Vert_2^2-\Vert x_{\mathrm{init}}\Vert_2^2),
\label{eq:constrainedOptimization}$$ which is a constrained optimization problem on the hypersphere $\Vert x\Vert_2=\Vert x_{\mathrm{init}}\Vert_2$ for a given input image $x_{\mathrm{init}}\in{\mathbb{R}}^{nC}$. The optimality condition of  with respect to $x$ is equivalent to $\nabla_1\mathcal{R}(\overline{x},\theta)=\overline{\lambda}\overline{x}$, which shows that $(\overline{x},\overline{\lambda})$ is actually a nonlinear eigenpair of $\nabla_1\mathcal{R}$. Figure \[fig:eigenpairs\] depicts ten triplets of input images, eigenmodes and eigenvalues $(x_{\mathrm{init}},\overline{x},\overline{\lambda})$, where the output images are computed using accelerated projected gradient descent. As a result, the algorithm generates cartoon-like eigenmodes and energy minimizing patterns are hallucinated due to the constraint. For instance, novel stripe patterns are prolongated in the body and facial region of the woman in the third column or in the body parts of the zebra in the fifth column. Furthermore, contours are axis aligned like the roof in the eighth column or the ship and the logo in the last image.

![Triplets of input images $x_{\mathrm{init}}$ (top), eigenmodes $\overline{x}$ (bottom) and eigenvalues $\overline{\lambda}$ associated with the TDV$_{25}^3$ regularizer.[]{data-label="fig:eigenpairs"}](figure2.pdf){width="\linewidth"}

![ Comparison of RMSE along two sequences computed with two different parameter sets with the estimate . []{data-label="fig:sensitivity"}](sensitivity.pdf){width=".5\linewidth"}

In Figure \[fig:sensitivity\], we compare the root mean squared error (RMSE) between the image sequences $x_s$ and $\widetilde{x}_s$ generated using TDV$^1$ with parameter sets $(T,\theta)$ and $(\widetilde{T},\widetilde{\theta})$ with the upper bound estimate in Theorem \[thm:sensitivity\]. Here, $(T,\theta)$ are computed by training on the entire BSDS400 training set, whereas we solely use 10 randomly selected images of this training set to obtain $(\widetilde{T},\widetilde{\theta})$. The local Lipschitz constant $L_\mathcal{R}$ in  is estimated directly from the sequences. As a result, the RMSE of the norm differences is roughly constant along $S$ and the upper bound estimate is efficient. Further details of the sensitivity analysis are presented in the supplementary material.

Computed tomography reconstruction
----------------------------------

To demonstrate the broad applicability of the proposed TDV, we perform two-dimensional computed tomography (CT) reconstruction using the TDV$^3_{25}$ regularizer trained for image denoising and $S=10$. We stress that the regularizer is applied *without* any additional training of the parameters.

The task of computed tomography is the reconstruction of an image given a set of projection measurements called sinogram, in which the detectors of the CT scanner measure the intensity of attenuated X-ray beams. Here, we use the linear attenuation model introduced in [@HaMu18], where the attenuation is proportional to the intersection area of a triangle, which is spanned by the X-ray source and a detector element, and the area of an image element. In detail, the sinogram $z$ of an image $x$ is computed by $z=A_R x$, where $A_R$ is the lookup-table based area integral operator of [@HaMu18] for $R$ angles and $768$ projections. Typically, a fully sampled acquisition consists of $2304$ angles. For this task, we consider the problem of angular undersampled CT [@ChTa08], where only a fraction of the angles are measured. We use a 4-fold ($R=576$) and 8-fold ($R=288$) angular undersampling to reconstruct a representative image of the MAYO dataset [@McBa17]. To account for an imbalance of regularization and data fidelity, we manually scale the data fidelity term by $\lambda>0$, i.e. $\mathcal{D}(x,z)\coloneqq\frac{\lambda}{2}\Vert A_Rx-z\Vert_2^2$. The resulting smooth variational problem is optimized using accelerated gradient descent with Lipschitz backtracking. Further details of the CT reconstruction task and the optimization scheme are included in the supplementary material.

We present qualitative and quantitative results for CT reconstruction in Figure \[fig:CT\] for a single abdominal CT image. As an initialization, we perform $50$ steps of a conjugate gradient method on the data fidelity term (first and last column). Using the proposed regularizer TDV$^3_{25}$, we are able to suppress the undersampling artifacts while preserving the fine vessels in the liver. This highlights that the learned regularizer can be effectively applied as a generic regularizer for linear inverse problems without any transfer learning.

![Conjugate gradient reconstruction for 4/8-fold angular undersampled CT task (first/fifth image), results obtained by using the TDV$^3_{25}$ regularizer with $\lambda=500\cdot10^3$ for 4-fold (second image) and $\lambda=10^6$ for 8-fold undersampling (fourth image), and fully sampled reference reconstruction using the TDV$^3_{25}$ regularizer with $\lambda=125\cdot10^3$ (third image).[]{data-label="fig:CT"}](figure3.pdf){width="\linewidth"}

Magnetic resonance imaging reconstruction
-----------------------------------------

Next, the flexibility of our regularizer TDV$_{25}^3$ learned for denoising and $S=10$ is shown for accelerated magnetic resonance imaging (MRI) *without* any further adaption of $\theta$.

In accelerated MRI, k–space data is acquired using $N_C$ parallel coils, each measuring a fraction of the full k–space to reduce acquisition time [@HaKl18]. Here, we use the data fidelity term $\mathcal{D}(x,\{z_i\}_{i=1}^{N_C})=\frac{\lambda}{2}\sum_{i=1}^{N_C}\Vert M_RFC_ix-z_i\Vert_2^2$, where $\lambda>0$ is a manually adjusted weighting parameter, $M_R$ is a binary mask for $R$-fold undersampling, $F$ is the discrete Fourier transform, and $C_i$ are sensitivity maps computed following [@UeLa14]. We use 4-fold and 6-fold Cartesian undersampled MRI data to reconstruct a sample knee image. Again, we minimize the resulting variational energy by accelerated gradient descent with Lipschitz backtracking. Further details of this task and the optimization are in the supplementary material.

Figure \[fig:MRI\] depicts qualitative results and PSNR values for the reconstruction of 4-fold and 6-fold undersampled k–space data. The first and last columns show the initial images obtained by applying the adjoint operator to the undersampled data. Due to the undersampling in k–space both images are severely corrupted by backfolding artifacts, which are removed by applying the proposed regularizer without losing fine details.

![Zero filling initialization for acceleration factors $R\in\{4,6\}$ (first/fifth image), output using the TDV$^3_{25}$ regularizer with $\lambda=1000$ for $R=4$ (second image) and $\lambda=1500$ for $R=6$ (fourth image), and fully sampled reference (third image).[]{data-label="fig:MRI"}](figure4.pdf){width="\linewidth"}

Single image super-resolution
-----------------------------

For single image super-resolution with scale factor $\gamma\in\{2,3,4\}$, we start with a full resolution ground truth image patch $y^i\in{\mathbb{R}}^{nC}$ uniformly drawn from the BSDS400 dataset. To obtain the linear downsampling operator $A$, we implement the adjoint operator matching MATLAB^^’s bicubic upsampling operator `imresize`, which is an implementation of a scale factor-dependent interpolation convolution kernel in conjunction with a stride. The observed low resolution image is $z^i=Ay^i\in{\mathbb{R}}^{nC/\gamma^2}$. As an initialization, $3$ iteration steps of the conjugate gradient method for the variational problem $x_{\mathrm{init}}^i=\operatorname*{argmin}_{x\in{\mathbb{R}}^{nC}}\Vert Ax-z^i\Vert_2^2$ are applied. The matrix inverse in  is approximated in each step by $7$ steps of a conjugate gradient scheme. Here, all results are obtained by training a TDV$^3$ regularizer for each scale factor individually.

In Table \[tab:SR\], we compare our approach with several state-of-the-art networks of similar complexity and list average PSNR values of the Y-channel in the YCbCr color space over test datasets. For the BSDS100 dataset, our proposed method achieves similar results as OISR-LF-s [@HeMo19] with only one third of the trainable parameters. Figure \[fig:OptimalStoppingSequenceSR\] depicts a restored sequence of images for the single image super-resolution task with scale factor $4$ using TDV$^3$ for a representative sample image of the Set14 dataset. Starting from the low resolution initial image, interfaces are gradually sharpened and the best quality is achieved for $T=\overline{T}$. Beyond this point, interfaces are artificially intensified.

![ From left to right: High resolution, low resolution with scale factor $4$, and resulting output of TDV$^3$ for $(S,T)\in\{(5,\frac{\overline{T}}{2}),(10,\overline{T}),(15,\frac{3\overline{T}}{2}),(20,2\overline{T})\}$, where the optimal stopping time is $\overline{T}=0.098$. Note that the best image is framed in red. []{data-label="fig:OptimalStoppingSequenceSR"}](figure5.pdf){width="\linewidth"}

Conclusion
==========

In this paper, we have introduced total deep variation regularizers, which are motivated by established deep network architectures. The inherent variational structure of TDVs enables a rigorous mathematical understanding encompassing an optimality condition for optimal stopping, a nonlinear eigenmode analysis, and a sensitivity analysis that yields computable upper error bounds. For image denoising and single image super-resolution, our model generates state-of-the-art results with an impressively low number of trainable parameters. Moreover, to underline the versatility of TDVs for generic linear inverse problems, we successfully demonstrated their applicability for the challenging CT and MRI reconstruction tasks without requiring any additional training.

Notation and mathematical preliminaries
=======================================

In this section, we briefly introduce the notation and give a short presentation of the mathematical preliminaries required for the proofs. All results can be found in [@AdFo03; @Te12].

Let $\Omega\subset{\mathbb{R}}^l$ be a bounded domain. We denote by $C^0(\overline\Omega,{\mathbb{R}}^n)$ the *set of continuous functions* and by $C^k(\overline\Omega,{\mathbb{R}}^n)$, $k\in{\mathbb{N}}$, the *set of $k$ times continuously differentiable functions* mapping form $\overline\Omega$ to ${\mathbb{R}}^n$. The corresponding norms are $$\begin{aligned}
\Vert f\Vert_{C^0(\overline\Omega,{\mathbb{R}}^n)}&\coloneqq\max_{x\in\overline\Omega}\Vert f(x)\Vert_2,\notag\\
\Vert f\Vert_{C^k(\overline\Omega,{\mathbb{R}}^n)}&\coloneqq\sum_{\vert\alpha\vert\leq k}\Vert D^\alpha f\Vert_{C^0(\overline\Omega,{\mathbb{R}}^n)},\end{aligned}$$ respectively. Here, we set $D^\alpha f=(\frac{\partial^{\alpha_1}f}{\partial x_1^{\alpha_1}},\ldots,\frac{\partial^{\alpha_n}f}{\partial x_{n}^{\alpha_n}})^\top$ for a multi-index $\alpha=(\alpha_1,\ldots,\alpha_n)\in{\mathbb{N}}_0^n$ with $\vert\alpha\vert=\sum_{i=1}^n\alpha_i$. We define the *Hölder space* as $$C^{k,s}(\overline{\Omega},{\mathbb{R}}^n)\coloneqq\{f\in C^k(\overline{\Omega},{\mathbb{R}}^n):\Vert f\Vert_{C^{k,s}(\overline{\Omega},{\mathbb{R}}^n)}<\infty\}$$ for $k\in{\mathbb{N}}$ and $s\in(0,1]$, where the Hölder norm is $$\Vert f\Vert_{C^{k,s}(\overline{\Omega},{\mathbb{R}}^n)}\coloneqq\Vert f\Vert_{C^k(\overline{\Omega},{\mathbb{R}}^n)}+\sup_{\substack{x,y\in\overline{\Omega}\\x\neq y}}\frac{\Vert f(x)-f(y)\Vert_2}{\Vert x-y\Vert_2^s}.$$ For $1\leq p<\infty$, the *Lebesgue space* $L^p(\Omega,{\mathbb{R}}^n)$ is the set of all measurable mappings from $\Omega$ to ${\mathbb{R}}^n$ such that $\Vert f\Vert_{L^p(\Omega,{\mathbb{R}}^n)}<\infty$, where $$\Vert f\Vert_{L^p(\Omega,{\mathbb{R}}^n)}\coloneqq\left(\int_\Omega\Vert f(x)\Vert_2^p{\,\mathrm{d}}x\right)^\frac{1}{p}.$$ The *Sobolev space* $H^m(\Omega,{\mathbb{R}}^n)$ with $m\in{\mathbb{N}}$ is the set of all functions $f\in L^2(\Omega,{\mathbb{R}}^n)$ such that for all $\alpha\in{\mathbb{N}}_0^n$ with $\vert\alpha\vert\leq m$ a function $f^{(\alpha)}\in L^2(\Omega,{\mathbb{R}}^n)$ exists satisfying $$\int_\Omega\langle D^\alpha\eta,f\rangle{\,\mathrm{d}}x=(-1)^{\vert\alpha\vert}\int_\Omega\langle\eta,f^{(\alpha)}\rangle{\,\mathrm{d}}x$$ for all smooth functions $\eta:\Omega\to{\mathbb{R}}^n$ with compact support. The space $H^m(\Omega,{\mathbb{R}}^n)$ is endowed with the norm $$\Vert f\Vert_{H^m(\Omega,{\mathbb{R}}^n)}\coloneqq\sum_{\vert\alpha\vert\leq m}\Vert f^{(\alpha)}\Vert_{L^2(\Omega,{\mathbb{R}}^n)}.$$ In what follows, we list several theorems for later reference:

\[thm:maximumDomainExistence\] We assume that for $x_{\mathrm{init}}\in{\mathbb{R}}^n$, $g\in C^0({\mathbb{R}}\times{\mathbb{R}}^n,{\mathbb{R}}^n)$ and each $T>0$ there exist constants $C_1(T),C_2(T)\in(0,\infty)$ such that $$\Vert g(t,x)\Vert_2\leq C_1(T)+C_2(T)\Vert x\Vert_2$$ for all $(t,x)\in[0,T]\times{\mathbb{R}}^n$. Then all solutions of the initial value problem $\dot{x}(t)=g(t,x)$ with $x(0)=x_{\mathrm{init}}$ are defined for all $t\geq 0$.

See [@Te12 Theorem 2.17].

\[thm:Gronwall\] Suppose that the functions $y,\alpha,\beta\in C^0([0,T],{\mathbb{R}})$ satisfy $$y(t)\leq\alpha(t)+\int_0^t\beta(s)y(s){\,\mathrm{d}}s$$ for all $t\in[0,T]$. If $\alpha$ is monotonically increasing (i.e. $\alpha(s)\leq\alpha(t)$ for $s\leq t\in[0,T]$) and $\beta\geq0$, then $$y(t)\leq\alpha(t)\exp\left(\int_0^t\beta(s){\,\mathrm{d}}s\right)$$ for all $t\in[0,T]$.

See [@Te12 Lemma 2.7].

Next, we present special cases of important results in functional analysis tailored to our particular application.

A sequence $\{g_k\}_{k\in{\mathbb{N}}}$ in $H^m(\Omega,{\mathbb{R}}^n)$ converges weakly to $g\in H^m(\Omega,{\mathbb{R}}^n)$ (denoted by $g_k\rightharpoonup g$) if for every element $g^\ast$ of the dual space of $H^m(\Omega,{\mathbb{R}}^n)$ $g^\ast(g_k)\to g^\ast(g)$ as $k\to\infty$.

\[thm:compactness\] Let $I\subset{\mathbb{R}}$ be a bounded and open interval and $\{g_k\}_{k\in{\mathbb{N}}}$ in $H^m(I,{\mathbb{R}}^n)$ be a uniformly bounded sequence in $H^m(I,{\mathbb{R}}^n)$, i.e. there exists a constant $C>0$ such that $\Vert g_k\Vert_{H^m(I,{\mathbb{R}}^n)}<C$ for all $k\in{\mathbb{N}}$. Then there exists a subsequence $\{g_{k_l}\}_{l\in{\mathbb{N}}}$ of $\{g_k\}_{k\in{\mathbb{N}}}$ and $g\in H^m(I,{\mathbb{R}}^n)$ such that $g_{k_l}\rightharpoonup g$ in $H^m(I,{\mathbb{R}}^n)$ as $l\to\infty$.

See [@Al16 Chapter 8].

\[thm:SobolevEmbedding\] Let $I\subset{\mathbb{R}}$ be a bounded and open interval. If $k,m\in{\mathbb{N}}$ and $\alpha\in(0,1)$ satisfy $m-\frac{1}{2}>k+\alpha$, then the operator ${\mathrm{Id}}:H^m(I,{\mathbb{R}}^n)\to C^{k,\alpha}(\overline{I},{\mathbb{R}}^n)$ is continuous and compact. In particular, for every $g\in H^m(I,{\mathbb{R}}^n)$ there exists a function $\hat{g}\in C^{k,\alpha}(\overline{I},{\mathbb{R}}^n)$ that coincides with $g$ almost everywhere such that $$\Vert \hat{g}\Vert_{C^{k,\alpha}(\overline{I},{\mathbb{R}}^n)}\leq C\Vert g\Vert_{H^m(I,{\mathbb{R}}^n)}.$$ The constant $C$ solely depends on $k$, $m$, $I$ and $\alpha$.

See [@Al16 Theorem 10.13].

Existence of solutions
======================

In this section, we briefly recall the functional analytic setting of the optimal control problem. Afterwards, we present the proof for the existence of minimizers of the optimal control problem.

Throughout this section, we consider the data fidelity term $\mathcal{D}(x,z)=\frac{1}{2}\Vert Ax-z\Vert_2^2$ for a fixed task-dependent $A\in{\mathbb{R}}^{lC\times nC}$ and fixed $z\in{\mathbb{R}}^{lC}$. The total deep variation $\mathcal{R}$ is parametrized by $\theta\in\Theta$, where $\Theta\subset{\mathbb{R}}^p$ is a compact and convex set of admissible training parameters. We highlight that the building blocks of $\mathcal{R}$ are convolutional layers at different scales and the smooth log-student-t-distribution of the form $\phi(x)=\frac{1}{2\nu}\log(1+\nu x^2)$ for $\nu>0$. In particular, this structure implies $$\begin{aligned}
\Vert\nabla_1\mathcal{R}(x,\theta)\Vert_2&\leq C_\mathcal{R}(\theta,\nu)\Vert x\Vert_2
\label{eq:growthEstimateR}\end{aligned}$$ for a constant $C_\mathcal{R}(\theta,\nu)>0$ depending on $\theta$ and $\nu$. We define the energy functional as $$\mathcal{E}(x,\theta,z)\coloneqq\mathcal{D}(x,z)+\mathcal{R}(x,\theta).$$ Then, the gradient flow for the finite time interval $(0,T)$ is $$\begin{aligned}
\dot{\tilde{x}}(t)=&-\nabla_1\mathcal{E}(\tilde{x}(t),\theta,z)=f(\tilde{x}(t),\theta,z)\\
\coloneqq&-A^\top(A\tilde{x}(t)-z)-\nabla_1\mathcal{R}(\tilde{x}(t),\theta)\end{aligned}$$ for $t\in(0,T)$, where $\tilde{x}(0)=x_{\mathrm{init}}$ for a fixed initial value $x_{\mathrm{init}}\in{\mathbb{R}}^{nC}$. Here, $\tilde{x}:[0,T]\to{\mathbb{R}}^{nC}$ denotes the differentiable flow of $\mathcal{E}$, and $T>0$ refers to the time horizon. We assume that $T\in[0,{T_\mathrm{max}}]$ for a fixed ${T_\mathrm{max}}>0$. To formulate optimal stopping, we exploit the reparametrization $x(t)=\tilde{x}(tT)$, which results for $t\in(0,1)$ in the equivalent gradient flow formulation $$\dot{x}(t)=Tf(x(t),\theta,z),\qquad x(0)=x_{\mathrm{init}}.$$ For fixed $N\in{\mathbb{N}}$, let $(x_{\mathrm{init}}^i,y^i,z^i)_{i=1}^N\in({\mathbb{R}}^{nC}\times{\mathbb{R}}^{nC}\times{\mathbb{R}}^{lC})^N$ be a collection of $N$ triplets of initial-target image pairs, and observed data independently drawn from a task-dependent fixed probability distribution. Let ${l}:{\mathbb{R}}^{nC}\to{\mathbb{R}}_0^+$ be a convex, twice continuously differentiable and coercive (i.e. $\lim_{\Vert x\Vert_2\to\infty}{l}(x)=+\infty$) function. Then, the sampled optimal control problem reads as $$\inf_{T\in[0,{T_\mathrm{max}}],\,\theta\in\Theta}
\left\{J(T,\theta)\coloneqq\frac{1}{N}\sum_{i=1}^N{l}(x^i(1)-y^i)\right\}
\label{eq:objectiveFunctionApp}$$ subject to the state equation for each sample $$\dot{x}^i(t)=Tf(x^i(t),\theta,z^i),\qquad x^i(0)=x_{\mathrm{init}}^i\label{eq:gradientFlowSampleApp}$$ for $i=1,\ldots,N$ and $t\in(0,1)$. The next theorem ensures the existence of solutions to this optimal control problem.

The minimum in  subject to the side conditions  is attained.

Without restriction, we consider the case $N=1$ and omit the superscript. Taking into account  we can estimate as follows: $$\Vert Tf(x,\theta,z)\Vert_2
\leq T(\Vert A\Vert_2\Vert z\Vert_2)+T(\Vert A\Vert_2^2+C_\mathcal{R}(\theta,\nu))\Vert x\Vert_2.
\label{eq:upperEstimateRHS}$$ Thus, Theorem \[thm:maximumDomainExistence\] implies that $[0,1]$ is contained in the maximum domain of existence of the state equation  for any initial value $x_{\mathrm{init}}\in{\mathbb{R}}^n$. Let $(T_j,\theta_j)\in[0,{T_\mathrm{max}}]\times\Theta$ be a minimizing sequence for $J$ with an associated state $x_j\in C^1([0,1],{\mathbb{R}}^n)$ such that  with $(T,\theta)$ replaced by $(T_j,\theta_j)$ holds true. Due to the compactness of $[0,{T_\mathrm{max}}]$ and $\Theta$ we can deduce the existence of a subsequence (not relabeled) such that $(T_j,\theta_j)\to(T,\theta)\in[0,{T_\mathrm{max}}]\times\Theta$. Then, Theorem \[thm:Gronwall\] in conjunction with  implies the uniform boundedness of $\Vert x_j(t)\Vert_2$ and thus also $\Vert\dot{x}_j(t)\Vert_2$ for all $t\in[0,1]$ and all $j\in{\mathbb{N}}$. In particular, the sequence $\{x_j\}_{j\in{\mathbb{N}}}$ is uniformly bounded in the Hilbert space $H^1((0,1),{\mathbb{R}}^n)$. Thus, taking into account Theorem \[thm:compactness\] we can infer that a subsequence (not relabeled) converges weakly in $H^1((0,1),{\mathbb{R}}^n)$ to $x\in H^1((0,1),{\mathbb{R}}^n)$. Theorem \[thm:SobolevEmbedding\] ensures the compact embedding of $H^1((0,1),{\mathbb{R}}^n)$ into $C^0([0,1],{\mathbb{R}}^n)$, which implies $x_j\to x$ in $C^0([0,1],{\mathbb{R}}^n)$ and $x(0)=x_{\mathrm{init}}$. Since $f$ is a smooth function, we can even conclude $x\in C^1([0,1],{\mathbb{R}}^n)$ due to the general theory of ordinary differential equations presented in [@Ha80 Chapter I]. Note that $x$ is actually the unique solution to . This follows from the convergence of $x_j\to x$ in $C^0([0,1],{\mathbb{R}}^n)$ and the smoothness of $f$. Finally, the theorem follows from the continuity of ${l}$, i.e. $$\lim_{j\to\infty}J(T_j,\theta_j)=\lim_{j\to\infty}{l}(x_j(1)-y)={l}(x(1)-y)=J(T,\theta),$$ where we used the uniform convergence of $x_j\to x$ in $C^0([0,1],{\mathbb{R}}^n)$.

Discretized optimal control problem
===================================

In this section, we present details for the fully discrete formulation of the sampled optimal control problem. The state equation  is discretized using a semi-implicit scheme resulting in $$x_{s+1}^i=x_s^i-\tfrac{T}{S}A^\top(Ax_{s+1}^i-z^i)-\tfrac{T}{S}\nabla_1\mathcal{R}(x_s^i,\theta)\in{\mathbb{R}}^{nC}
\label{eq:stateEquationDiscreteApp}$$ for $s=0,\ldots,S-1$ and $i=1,\ldots,N$, where the *depth* $S\in{\mathbb{N}}$ is a priori fixed. This equation is equivalent to $x_{s+1}^i=\widetilde{f}(x_s^i,T,\theta,z^i)$ with $$\begin{aligned}
\widetilde{f}(x,T,\theta,z)&\coloneqq B(T)^{-1}(x+\tfrac{T}{S}(A^\top z-\nabla_1\mathcal{R}(x,\theta))),\notag\\
B(T)&\coloneqq{\mathrm{Id}}+\frac{T}{S}A^\top A.\end{aligned}$$ The initial state satisfies $x_0^i=x_{\mathrm{init}}^i\in{\mathbb{R}}^{nC}$. Then, the discretized sampled optimal control problem reads as $$\inf_{T\in[0,{T_\mathrm{max}}],\,\theta\in\Theta}
\left\{J_S(T,\theta)\coloneqq\frac{1}{N}\sum_{i=1}^N{l}(x_S^i-y^i)\right\}
\label{eq:discreteOptimalControlApp}$$ subject to $x_{s+1}^i=\widetilde{f}(x_s^i,T,\theta,z^i)$. We define the Lagrange functional $L_S:({\mathbb{R}}^{nC})^{N(S+1)}\times[0,{T_\mathrm{max}}]\times\Theta\times({\mathbb{R}}^{lC})^N\times({\mathbb{R}}^{nC})^{N(S+1)}\to{\mathbb{R}}$ as follows: $$(x,T,\theta,z,p)\mapsto J_S(T,\theta)+\sum_{i=1}^N\langle p_0^i,x_0^i-x_{\mathrm{init}}^i\rangle
+\sum_{s=0}^{S-1}\sum_{i=1}^N\langle p_{s+1}^i,x_{s+1}^i-\widetilde{f}(x_s^i,T,\theta,z^i)\rangle.$$ The discretization of the associated adjoint equation is now implied by the discrete Pontryagin maximum principle:

Let $(\overline{T},\overline{\theta})$ be a pair of optimal control parameters for  with the state equation $\{\overline{x}_s^i\}_{s=0,\ldots,S}^{i=1,\ldots,N}$. We define the Hamiltonian $$\begin{aligned}
H&:{\mathbb{R}}^{nC}\times{\mathbb{R}}^{nC}\times[0,{T_\mathrm{max}}]\times\Theta\times{\mathbb{R}}^{lC}\to{\mathbb{R}},\notag\\
(x,p,T,\theta,z)&\mapsto\langle p,\widetilde{f}(x,T,\theta,z)\rangle.\end{aligned}$$ If it is further assumed that $\nabla \widetilde{f}(\overline{x}_s^i,\overline{T},\overline{\theta},z^i)$ has full rank for all $i=1,\ldots,N$ and $s=0,\ldots,S$, then there exists an adjoint process $\{\overline{p}_s^i\}_{s=0,\ldots,S}^{i=1,\ldots,N}$ such that $$\begin{aligned}
\overline{x}_{s+1}^i&=\nabla_2 H(\overline{x}_s^i,\overline{p}_{s+1}^i,\overline{T},\overline{\theta},z^i),\notag\\
\overline{x}_0^i&=x_{\mathrm{init}}^i,\notag\\
\overline{p}_s^i&=\nabla_1 H(\overline{x}_s^i,\overline{p}_{s+1}^i,\overline{T},\overline{\theta},z^i),\notag\\
\overline{p}_S^i&=-\frac{1}{N}\nabla{l}(\overline{x}_S^i-y^i)\label{eq:PMP}.\end{aligned}$$ Finally, the solution is optimal in the sense that $$\sum_{i=1}^N H(\overline{x}_s^i,\overline{p}_{s+1}^i,\overline{T},\overline{\theta},z^i)
\geq \sum_{i=1}^N H(\overline{x}_s^i,\overline{p}_{s+1}^i,T,\theta,z^i)$$ for all $T\in[0,{T_\mathrm{max}}]$ and $\theta\in\Theta$.

For a complete proof, see [@LiHa18 Theorem 3] and [@Ha66]. In what follows, we give a heuristic justification for the definition of the Hamiltonian. To this end, we consider the Lagrange functional $L_S$ with multipliers $\{p_s^i\}_{s=0,\ldots,S}^{i=1,\ldots,N}$: $$\begin{aligned}
L_S(x,T,\theta,z,p)
=\frac{1}{N}\sum_{i=1}^N{l}(x_S^i-y^i)+\sum_{i=1}^N\langle p_0^i,x_0^i-x_{\mathrm{init}}^i\rangle
+\sum_{s=0}^{S-1}\sum_{i=1}^N\left(\langle p_{s+1}^i,x_{s+1}^i\rangle-H(x_s^i,p_{s+1}^i,T,\theta,z^i)\right).\notag\end{aligned}$$ Taking the derivatives with respect to $\{x_s^i\}_{s=0,\ldots,S}^{i=1,\ldots,N}$ and $\{p_s^i\}_{s=0,\ldots,S}^{i=1,\ldots,N}$ yields .

Thus, readily implies $$\begin{aligned}
\overline{p}_s^i&=\nabla_1 H(\overline{x}_s^i,\overline{p}_{s+1}^i,\overline{T},\overline{\theta},z^i)
=({\mathrm{Id}}-\tfrac{\overline{T}}{S}\nabla_1^2\mathcal{R}(\overline{x}_s^i,\overline{\theta}))({\mathrm{Id}}+\tfrac{\overline{T}}{S}A^\top A)^{-1}\overline{p}_{s+1}^i,\notag\\
\overline{p}_S^i&=-\frac{1}{N}\nabla{l}(\overline{x}_S^i-y^i).
\label{eq:adjointEquationDiscreteApp}\end{aligned}$$ The next theorem states an exact computable condition for the optimal stopping time, which is of vital importance for the numerical optimization:

Let $(\overline{T},\overline{\theta})$ be a stationary point of $J_S$ with associated states $\overline{x}_s^i$ and adjoint states $\overline{p}_s^i$ satisfying  and . We further assume that $\nabla \widetilde{f}(\overline{x}_s^i,\overline{T},\overline{\theta},z^i)$ has full rank for all $i=1,\ldots,N$ and $s=0,\ldots,S$. Then, we have $$-\frac{1}{N}\sum_{s=0}^{S-1}\sum_{i=1}^N\langle\overline{p}_{s+1}^i,B(\overline{T})^{-1}(\overline{x}_{s+1}^i-\overline{x}_s^i)\rangle=0.
\label{eq:optimalTDiscreteApp}$$

We compute the derivative of $\widetilde{f}(x_s^i,T,\theta,z^i)$ with respect to $T$: $$\begin{aligned}
\nabla_2\widetilde{f}(x_s^i,T,\theta,z^i)
=&\tfrac{1}{S}B(T)^{-1}\Big(A^\top z-\nabla_1\mathcal{R}(x_s^i,\theta)
-A^\top AB(T)^{-1}(x_s^i+\tfrac{T}{S}(A^\top z-\nabla_1\mathcal{R}(x_s^i,\theta)))\Big)\notag\\
=&\frac{1}{T}B(T)^{-1}(x_{s+1}^i-x_s^i),\end{aligned}$$ where we use $$\tfrac{{\,\mathrm{d}}}{{\,\mathrm{d}}T}(B(T)^{-1})=-B(T)^{-1}\left(\tfrac{{\,\mathrm{d}}}{{\,\mathrm{d}}T}B(T)\right)B(T)^{-1}.$$ The existence of the Lagrange multiplier $p$ is readily implied by the rank condition of $\nabla\widetilde{f}$ [@ItKu08 Chapter 1]. At the optimum $(\overline{x},\overline{T},\overline{\theta},z,\overline{p})$ the Lagrange functional $L_S$ satisfies the following optimality condition with respect to $T$: $$\nabla_2 L_S(\overline{x},\overline{T},\overline{\theta},z,\overline{p})
=-\sum_{s=0}^{S-1}\sum_{i=1}^N\langle\overline{p}_{s+1}^i,\nabla_2\widetilde{f}(\overline{x}_s^i,\overline{T},\overline{\theta},z^i)\rangle
=-\sum_{s=0}^{S-1}\sum_{i=1}^N\langle\overline{p}_{s+1}^i,\tfrac{1}{\overline{T}}B(\overline{T})^{-1}(\overline{x}_{s+1}^i-\overline{x}_s^i)\rangle,$$ which readily implies .

In the case of image denoising, the scaling $B(\overline{T})$ in the optimality condition  is neglected.

An important property of any learning based method is the dependency of the learned parameters with respect to different training datasets [@GeBi92].

Let $(T,\theta),(\widetilde{T},\widetilde{\theta})$ be two pairs of control parameters obtained from two different training datasets. We denote by $x,\widetilde x\in({\mathbb{R}}^{nC})^{(S+1)}$ two solutions of the state equation with the same observed data $z$ and initial condition $x_{\mathrm{init}}$, i.e. $$x_{s+1}=\widetilde{f}(x_s,T,\theta,z),\quad
\widetilde{x}_{s+1}=\widetilde{f}(\widetilde{x}_s,\widetilde{T},\widetilde{\theta},z)$$ for $s=1,\ldots,S-1$ and $x_0=\widetilde{x}_0=x_{\mathrm{init}}$. Let $L_\mathcal{R}$ be the Lipschitz constant of $\mathcal{R}$, i.e. $$\Vert\nabla_1\mathcal{R}(x,\theta)-\nabla_1\mathcal{R}(\widetilde{x},\widetilde{\theta})\Vert_2\leq L_\mathcal{R}
\left\Vert
\begin{pmatrix}
x\\
\theta
\end{pmatrix}
-
\begin{pmatrix}
\widetilde{x}\\
\widetilde{\theta}
\end{pmatrix}
\right\Vert_2$$ for all $x,\widetilde{x}\in{\mathbb{R}}^{nC}$ and all $\theta,\widetilde{\theta}\in\Theta$. Then, $$\begin{aligned}
\Vert x_{s+1}-\widetilde{x}_{s+1}\Vert_2
\leq&\Vert B(T)^{-1}-B(\widetilde{T})^{-1}\Vert_2
\Big(\Vert x_s\Vert_2+\tfrac{T}{S}\Vert A^\top z\Vert_2\notag\\
&+\tfrac{T}{S}\Vert\nabla_1\mathcal{R}(x_s,\theta)\Vert_2\Big)+\Vert B(\widetilde{T})^{-1}\Vert_2
\Big(\Vert x_s-\widetilde{x}_s\Vert_2\notag\\
&+\tfrac{\vert T-\widetilde{T}\vert}{S}\Vert A^\top z\Vert_2+\tfrac{\vert T-\widetilde{T}\vert}{S}\Vert\nabla_1\mathcal{R}(\widetilde{x}_s,\widetilde{\theta})\Vert_2
+\tfrac{T}{S}L_\mathcal{R}\Vert(x,\theta)^\top-(\widetilde{x},\widetilde{\theta})^\top\Vert_2\Big).\end{aligned}$$

Let $g(x,T,\theta)\coloneqq x+\frac{T}{S}A^\top z-\frac{T}{S}\nabla_1\mathcal{R}(x,\theta)$. By definition, $$\begin{aligned}
x_{s+1}&=B(T)^{-1}g(x_s,T,\theta),\notag\\
\widetilde{x}_{s+1}&=B(\widetilde{T})^{-1}g(\widetilde{x}_s,\widetilde{T},\widetilde{\theta}).\end{aligned}$$ Thus, we can estimate the upper bound as follows: $$\begin{aligned}
\Vert x_{s+1}-\widetilde{x}_{s+1}\Vert_2
=&\Vert B(T)^{-1}g(x_s,T,\theta)-B(\widetilde{T})^{-1}g(\widetilde{x}_s,\widetilde{T},\widetilde{\theta})\Vert_2\notag\\
\leq&\Vert B(T)^{-1}-B(\widetilde{T})^{-1}\Vert_2
\Big(\Vert x_s\Vert_2+\tfrac{T}{S}\Vert A^\top z\Vert_2
+\tfrac{T}{S}\Vert\nabla_1\mathcal{R}(x_s,\theta)\Vert_2\Big)+\Vert B(\widetilde{T})^{-1}\Vert_2
\Big(\Vert x_s-\widetilde{x}_s\Vert_2\notag\\
&+\tfrac{\vert T-\widetilde{T}\vert}{S}\Vert A^\top z\Vert_2+\tfrac{\vert T-\widetilde{T}\vert}{S}\Vert\nabla_1\mathcal{R}(\widetilde{x}_s,\widetilde{\theta})\Vert_2
+\tfrac{T}{S}L_\mathcal{R}\Vert(x_s,\theta)^\top-(\widetilde{x}_s,\widetilde{\theta})^\top\Vert_2\Big).\notag\end{aligned}$$ This concludes the proof.

Hence, this theorem provides a computable upper bound for the norm difference of two states with observed value $z$ and initial value $x_{\mathrm{init}}$ evaluated at the same step $s$, which amounts to a sensitivity analysis w.r.t. the training data.

Additional numerical results
============================

Figure \[fig:denoising\] depicts four additional results for the image denoising task computed with TDV$_{25}^3$ for a noise level $\sigma=25$. The best results in terms of the PSNR value are obtained for $S=10$ and are highlighted in red.

![From left to right: Ground truth, noisy input with noise level $\sigma=25$ and resulting output of TDV$^3$ for $(S,T)\in\{(5,\frac{\overline{T}}{2}),(10,\overline{T}),(15,\frac{3\overline{T}}{2}),(20,2\overline{T})\}$, where the optimal stopping time is $\overline{T}=0.0297$. Note that the best image is framed in red.[]{data-label="fig:denoising"}](figure6.pdf){width="\linewidth"}

Figure \[fig:sr\] shows four image sequences associated with image super-resolution with a scale factor $4$. The algorithm generates a transition from the initial upsampled ($S=0$) to a cartoon-like image with strongly pronounced edges ($S=20$). Here, the optimal stopping time is image dependent and determines the best image in terms of PSNR value.

![ From left to right: High resolution, low resolution with scale factor $4$, and resulting output of TDV$^3$ for $(S,T)\in\{(5,\frac{\overline{T}}{2}),(10,\overline{T}),(15,\frac{3\overline{T}}{2}),(20,2\overline{T})\}$, where the optimal stopping time is $\overline{T}=0.098$. []{data-label="fig:sr"}](figure7.pdf){width="\linewidth"}

Transferring TDV to different linear inverse problems
=====================================================

In this section, we elaborate on how the prior information incorporated in the TDV regularizer can be transferred to different linear inverse problems *without* any adaption of the learnable parameters.

As before, we consider minimizing the variational energy $$\mathcal{E}(x,\theta,z)\coloneqq\mathcal{D}(x,z)+\mathcal{R}(x,\theta),
\label{eq:energy}$$ in which the problem specific data fidelity term $$\mathcal{D}(x,z) = \frac{\lambda}{2}\Vert Ax-z\Vert_2^2$$ is used to estimate an image $x\in{\mathbb{R}}^{nC}$ from observed measurements $z\in{\mathbb{R}}^{lC}$. The linear operator $A\in{\mathbb{R}}^{lC\times nC}$ models the data acquisition process of the inverse problem at hand. The scalar weight $\lambda>0$ allows to manually balance the weighting between data fidelity and regularization. To incorporate prior knowledge in the reconstruction process , we consider the smooth TDV regularizer $\mathcal{R}(x,\theta)$. Here, we assume that the TDV parameters $\theta$ have already been determined by a previous training process for a different inverse problem (e.g. image denoising) as presented in the paper. As a result, we are able to transfer the regularization strategies obtained by training from data on a different task to another specific restoration/reconstruction problem by just using the learned regularizer.

For given measurements $z\in{\mathbb{R}}^{lC}$ and previously computed trainable parameters $\theta\in\Theta$ we retrieve the restored image $x$ via $$\operatorname*{argmin}_{x\in{\mathbb{R}}^{nC}}\mathcal{E}(x,\theta,z).
\label{eq:min}$$ Since both the data fidelity term and the regularization term are smooth, any convergent first-order gradient method is suitable for optimizing . Algorithm \[algo:agdlb\] lists the applied accelerated gradient descent method with Lipschitz backtracking. We adopt a Lipschitz backtracking strategy to approximate the local Lipschitz constant. Algorithm \[algo:agdlb\] is known to converge to stationary points of  if the number of iterations $K$ is chosen large enough [@ChPo16; @PoSa16].

To sum up, the idea of reusing prior knowledge extracted from data for a different inverse problem by means of a regularizer enables a simple and effective approach that does not require any further training. Next, we demonstrate the applicability of this idea by applying a regularizer learned for gray-scale image denoising to computed tomography (CT) and magnetic resonance imaging (MRI) reconstruction.

Computed tomography
===================

To showcase the effectiveness of the transfer idea introduced in the previous section, we consider the linear inverse problem of angular undersampled two dimensional computed tomography (CT). The task of CT is to reconstruct an image given projection measurements for different acquisition angles. All these measurements are stacked along the angular dimension to form the sinogram $z$. In practice typically $R=2304$ acquisition angles with $768$ projections are considered, thus $z\in{\mathbb{R}}^{768R}$. In the case of angular undersampled CT [@ChTa08] only a fraction of the acquisition angles are measured. Here, the goal is to reconstruct the image with a similar quality by incorporating prior knowledge in the form of regularization.

The angular undersampled CT problem can be addressed by considering the data fidelity term $$\mathcal{D}(x,z)=\frac{\lambda}{2}\Vert A_Rx - z\Vert_2^2,$$ where $A_R: {\mathbb{R}}^{R\cdot768\times768\cdot768}$ is the linear projection operator acquiring $R$ angles introduced in [@HaMu18] and $z=A_Ry$ is the measured sinogram of the true image $y$. We use the TDV$_{25}^3$ regularizer trained for image denoising and $S=10$ without changing its learnable parameters.

We evaluate the qualitative performance of the proposed regularizer transfer approach on a sample slice of the MAYO dataset [@McBa17]. The reconstructed images from the fully-sampled, 4/8-fold undersampled sinogram are depicted in Figure 8 in the paper. It is clearly visible that the TDV$_{25}^3$ regularizer yields reconstructed images of similar quality. Consequently, the regularizer is able to remove the undersampling artifacts while preserving fine details.

Magnetic resonance imaging
==========================

In this section, we demonstrate the wide-ranging usability of the proposed approach by transferring the TDV regularizer trained for image denoising to parallel accelerated MRI, which exhibits strong and structured undersampling artifacts.

In detail, in parallel accelerated MRI only a subset of the k-space data of each coil is measured during the acquisition process (for details see [@HaKl18]). The task in parallel accelerated MRI is the reconstruction of the scanned image $x$ given the undersampled k-space data $z_i$ of each receiver coil $i=1,\ldots,N_C$. To account for multiple coils, we change the data fidelity term to $$\mathcal{D}(x,\{z\}_{i=1}^{N_C})=\frac{\lambda}{2}\sum_{i=1}^{N_C}\Vert M_RFC_ix - z_i\Vert_2^2.$$ The linear operator $C_i\in{\mathbb{C}}^{n\times n}$ weights the complex image estimate by the coil sensitivity map of the $i^{th}$ coil, which was computed by [@UeLa14], $F\in{\mathbb{C}}^{n\times n}$ is the Fourier transform, and $M_R\in{\mathbb{C}}^{n\times n}$ is an binary mask for $R$-fold undersampling. Here, we consider Cartesian undersampled k-space data as in [@HaKl18]. As in the previous sections, the scalar weight $\lambda$ needs to be adapted to balance data fidelity and regularization. Despite the large difference between noise artifacts and the undersampling artifacts introduced by Cartesian undersampling, we use the TDV$_{25}^3$ regularizer trained for image denoising and $S=10$ to regularize this linear inverse problem.

We perform a qualitative evaluation of the proposed approach on a representative slice of an undersampled MRI knee image. The slice has a resolution of $n=320\cdot320$ pixels, and $N_C=15$ receiver coils were used during the acquisition. Figure 9 in the paper depicts the results of the accelerated parallel MRI problem using the TDV$_{25}^3$ regularizer. Although the TDV regularizer was not trained to account for undersampling artifacts, all these artifacts are removed in the reconstructions and only some details in the bone are lost. This highlights the versatility and effectiveness of the proposed TDV regularizer since both challenging medical inverse problems can be properly addressed *without* any fine-tuning of the learned parameters.
