Task adapted reconstruction for inverse problems

Only seeking an estimate of $x^{*}\in X$ is limiting since it does not account for the uncertainty. A more comprehensive analysis is based on introducing a $X$-valued random variable $\mathsf{x}\sim\pi^{*}$ whose true (unknown) probability distribution $\pi^{*}\in\mathscr{P}_{X}$ generates $x^{*}$. One can then rephrase the inverse problem stated earlier as the task of recovering the probability measure $\pi^{*}\in\mathscr{P}_{X}$ given data $y\in Y$ generated by $\mathsf{y}$, which is related to $x^{*}$ through the data model as in eq. 3. An important special case is when $\pi^{*}$ is parametrized by $x^{*}\in X$ in a known way, so the inverse problem reduces to the task of recovering $x^{*}\in X$.

In a Bayesian setting, one considers the posterior distribution of $\mathsf{x}$ given $\mathsf{y}=y$ up to a constant of proportionality. More precisely, consider a setting where the joint law $(\mathsf{x},\mathsf{y})\sim\mu$ can be written in terms of conditional probabilities:

Here, $\pi_{0}$ serves as a (possibly improper) prior and the last equality in eq. 5 follows from the definition of the data model as the conditional distribution of $\mathsf{y}$ given $\mathsf{x}=x^{*}$. In particular, the joint law $\mu$ in eq. 5 is proportional to the posterior, so the decomposition above exists as soon as Bayes’ theorem holds. This is the case in a rather general setting [18, Theorem 14], but a decomposition is also possible is some cases where the prior is not proper¹¹ 1 Under certain circumstances it is possible to work with improper priors on the model parameter space, e.g., by computing posterior distributions that approximate the posteriors one would have obtained using proper conjugate priors whose extreme values coincide with the improper prior..

A key point in the Bayesian setting is to explore the posterior distribution of $\mathsf{x}$ given $\mathsf{y}=y$ assuming that both $x\mapsto\pi_{0}\in\mathscr{P}_{X}$ (prior) and $x\mapsto\mathcal{M}(x)$ (data model) are known, but $x^{*}\in X$ is unknown. The data model often has an associate density $\mathcal{L}$ (data likelihood) that is known to sufficient degree of accuracy, in which case $\,\mathrm{d}\mathcal{M}(x)(y)=\mathcal{L}(y\mid x)\,\mathrm{d}y$.

A reconstruction method is formally a measurable $X$-valued mapping on $Y$, which in the statistical setting corresponds to an estimator. More precisely, $\bigl((Y,\mathfrak{S}_{Y}),\{\mathcal{M}(x)\}_{x\in X}\bigr)$ defines a statistical model parametrized by the model parameter space $X$ and a reconstruction method corresponds to a point estimator. The latter is a non-randomized decision rule for a statistical estimation problem where the model parameter space $X$ parametrizes the underlying statistical model and at the same time constitutes the decision space. The reader may here consult [60, section 3.1] for formal definitions of decision theoretic notions used here.

There are many possible reconstruction methods (estimators) so one needs a framework where these can be compared against each other. Statistical decision theory offers such a framework by associating a notion of risk to a decision rule. This quantifies the downside that comes with using a particular reconstruction method. The first step is to define the loss function on the decision space, which in our specific setting becomes a measurable mapping (see [60, Definition 3.2] for the definition in a general setting):

A common choice in imaging inverse problems is the $L^{2}$-loss, which is the squared $L^{2}$-distance. There are however alternatives that are not based on point-wise differences but on differences between high-level image features, e.g., the Wasserstein distance [3] and perceptual losses [46].

Having selected a loss function as in eq. 6 and a prior $\pi_{0}$ in eq. 5 on the model parameter space, the $\pi_{0}$-average risk (Bayes risk or expected loss) for reconstruction is given as

A natural criteria to select a reconstruction method (estimator) is to minimize Bayes risk, i.e., to select an estimator (non-randomized decision rule) that minimizes $\ForwardOp^{\dagger}\mapsto\mathcal{R}_{\pi_{0}}(\ForwardOp^{\dagger})$ in eq. 7.

Note here that in the finite dimensional setting, minimizing Bayes risk is the same as computing the conditional mean (posterior mean) if and only if the loss function in eq. 6 is the Bregman distance of a strictly convex non-negative differentiable functional [6]. This holds in particular when the loss function is given by the squared $L^{2}$-norm. Next, another common choice is the maximum a posteriori estimator that maximizes the posterior, so it corresponds to the most likely reconstruction given the data. On the other hand, a maximum likelihood estimator maximizes the negative log-likelihood of data, i.e., it corresponds to the model parameter that generates the most likely data. This is an unsuitable estimator in ill-posed inverse problems since it frequently leads to overfitting.

To summarize, we will henceforth consider a reconstruction method that minimizes Bayes risk and, as already mentioned, this equals the conditional mean when using a $L^{2}$-loss.

In the Bayesian setting (section 5.1), both the true model parameter and data are assumed to be generated by random variables, and the goal is to recover the conditional probability of the model parameter given data (posterior) [48, 25, 88, 18, 13]. In contrast, classical (deterministic) approaches view an inverse problem as an operator equation of the type eq. 1 [23, 49, 85, 52] where data may be generated by a random variable, but there are no statistical assumptions on model parameters.

The Bayesian viewpoint offers a more complete analysis than the classical approach that is based on solving eq. 1 in the sense that the posterior describes all possible solutions. In particular, different reconstructions can be obtained by using different estimators and there is a natural framework for uncertainty quantification, e.g., by computing Bayesian credible sets. Furthermore, small changes in the data lead to small changes in the posterior distribution in a fairly general setting [18, Theorem 16] (continuity of the posterior distribution in the Hellinger metric), so working with probability measures on the model parameter space (posterior) and adopting a suitable prior stabilizes an ill-posed inverse problem.

The posterior is, on the other hand, often quite complicated with no closed form expression. Much of the contemporary research therefore focuses on realizing the above advantages with Bayesian inversion without having access to the full posterior. Key topics are designing a “good” prior $\pi_{0}\in\mathscr{P}_{X}$ and to have computationally feasible means for exploring the posterior.

As outlined in section 5.3, there are two challenges associated with using Bayesian inversion: selecting a “good” prior (section 5.3.1) and providing a computationally feasible approach for computing suitable estimators, like the one that minimizes Bayes risk (section 5.3.2).

As we outline here, learned iterative methods address both these challenges. It makes use of techniques from machine learning, and deep neural networks in particular, which have demonstrated a remarkable capacity in capturing intricate relations from example data [57]. A key element is usage of highly parametrized generic models that can be adapted to specific decision rules, such as reconstruction by eq. 7, by training against example data. Learned iterative methods use a deep neural network to define an estimator (reconstruction method) that minimizes Bayes risk while accounting for the knowledge about how data is generated.

To give a more precise description, consider the joint law $\mu=\pi_{0}\otimes\mathcal{M}(x)$ in eq. 7 used for defining Bayes risk. In most practical applications, this joint law is unknown. Often one may however have access to the corresponding empirical measure given by supervised training data $(x_{1},y_{1}),\ldots,(x_{m},y_{m})\in X\times Y$ generated by $(\mathsf{x},\mathsf{y})\sim\mu$. This avoids introducing a handcrafted prior $\pi_{0}\in\mathscr{P}_{X}$. Furthermore, searching over all non-randomized decision rules is computationally unfeasible. Instead, we restrict our attention to those given by a (deep) neural network architecture, which are known to have large capacity (can approximate any Borel measurable mapping arbitrarily well [75]) and there are computationally feasible implementations. To summarize, we have a family of reconstruction methods $\ForwardOp_{\theta}^{\dagger}\colon Y\to X$ parametrized by a finite dimensional parameter set $\Theta$ and the optimal one is given by solving the training problem

The above approach for defining a reconstruction operator $\ForwardOp_{\theta^{*}}^{\dagger}\colon Y\to X$ is fully data driven in the sense that neither a prior on model parameter space nor a data model are handcrafted beforehand. Instead, all information is derived from the training data, which in particular does not utilize knowledge about how data is generated. This becomes a serious issue when the number of independent samples in training data are low compared to number of unknowns, which is commonly the case in imaging. Next, in many inverse problem the data model $x\mapsto\mathcal{M}(x)$ that describes how data is generated is known. Thus, it is unnecessarily pessimistic to disregard this information as in a fully data driven approach to reconstruction.

Learned iterative schemes [1, 2] define a non-linear reconstruction operator parametrized by a deep convolutional neural network architecture that accounts for the data model, or more precisely the data likelihood. The idea is to unroll a fixed point iterative scheme relevant for solving the inverse problem and replace the explicit iterative updating rule with a learned one given by a deep convolutional residual network. The approach can be formulated as a general scheme for solving (possibly non-linear) inverse problems [1, 2, 35], see also [65, 66, 20] for a formulation that learns proximal updates in linear inverse problems. This results in a computationally feasible approach with surprisingly low requirements on training data and good generalization properties that outperforms state-of-the-art image reconstruction in computed tomography [1, 2, 35], magnetic resonance imaging [64, 65, 37, 66], photoacoustic tomography [38], and superresolution [64, 65, 20].

Let $\bigl((X,\mathfrak{S}_{X}),\{\mathbb{P}_{z}\}_{z\in\triangle}\bigr)$ be a statistical model where the model parameter space $(X,\mathfrak{S}_{X})$ is a measurable space and $\{\mathbb{P}_{z}\}_{z\in\triangle}\subset\mathscr{P}_{X}$ is some family of probability measures on $X$ parametrized by elements in $\triangle$. Next, there is a measurable space $(D,\mathfrak{S}_{D})$ (decision space) and a fixed (task adapted) loss function ([60, Definition 3.2])

with given $\featuremap\colon\triangle\to D$ and $\ell_{D}\colon D\times D\to\mathbb{R}$ (decision distance). The statistical model along with the decision space and loss function defines a statistical estimation problem. Many tasks can now be seen as an appropriate non-randomized decision rule $\mathcal{T}\colon X\to D$ (task operator).

Before proceeding, it is worth reflecting over the roles of the above sets. In our set-up, the decision making associated with the task is based on elements in the decision space $D$ whereas actual observables are elements in $X$, so the task is represented by a measurable mapping $\mathcal{T}\colon X\to D$ (task operator). Often it is more natural to formalize the task as a mapping $\featuremap\colon\triangle\to D$ where elements in the set $\triangle$ are related to those in $X$. A difficulty is that elements in $\triangle$ are not observable and the mapping relating its elements to those in $X$ is unknown. Hence, the challenge is to infer an appropriate mapping $\mathcal{T}$ given $\featuremap$ by resorting to some suitable “optimality” principle. The examples in section 6.2 will further clarify the various roles of these sets in decision making.

Just as in section 5.2, we consider a decision rule that minimizes Bayes risk. More precisely, assume $\triangle$ is itself a measurable space and consider a fixed probability measure $\eta_{0}\in\mathscr{P}_{\triangle}$ (task prior). The task operator is the non-randomized decision rule $\mathcal{T}\colon X\to D$ that minimizes the associated Bayes risk:

A difficulty is to provide a ‘reasonable’ task prior $\eta_{0}\in\mathscr{P}_{\triangle}$. Another is that $\mathbb{P}_{z}\in\mathscr{P}_{X}$ is not known. Hence, one needs to consider the joint law $\eta:=\eta_{0}\otimes\mathbb{P}_{z}$ in eq. 10 as an unknown. Note that this differs from reconstruction, where the joint law is either known (as in section 5.1), or the prior is unknown but the data likelihood is known (as in section 5.4). Since the joint law is unknown, we replace it by the empirical measure given by (supervised) training data $(z_{1},x_{1}),\ldots,(z_{m},x_{m})\in\triangle\times X$, i.e., one has i.i.d. samples generated by a $(\triangle\times X)$-valued random variable $(\mathsf{z},\mathsf{x})\sim\eta$. Furthermore, due to issues associated with computational feasibility (section 5.3.2), we consider a parametrized family of decision rules $\mathcal{T}_{\vartheta}\colon X\to D$ given by a (deep) neural network architecture. Then, the task operator is the decision rule $\mathcal{T}_{\vartheta^{*}}\colon X\to D$ parametrized by a finite dimensional parameter in $\Xi$ and the optimal one $\vartheta^{*}\in\Xi$ is given by empirical risk minimization:

We conclude with examples showing how a wide range of image processing tasks can be phrased as decision rules in a statistical decision problem.

The abstract framework in section 6.1 for formalizing a task on model parameter space is very generic and covers a wide range of possible tasks. In the following, we list concrete examples from imaging in order to show how this framework can be adapted to specific cases. To ensure a computational feasible implementation, our focus is on tasks that have been successfully addressed using techniques from deep learning. Deep learning has proven to be an efficient computational framework for many tasks, much thanks to its ability to progressively learn discriminative hierarchal features of the input data by means of training a suitable deep neural network. Hence, this limitation is not as restrictive as it may seem at a first glance, which will also become evident by the examples listed here and in section 6.3.

Unless otherwise stated, tasks are formulated for grey-scale images defined on a fixed domain $\Omega\subset\mathbb{R}^{n}$, i.e., $X:=L^{2}(\Omega,\mathbb{R})$. We will also assume that $X$ is a measurable space for some $\sigma$-algebra $\mathfrak{S}_{X}$. Finally, $\mathscr{M}$ denotes the space of measurable mappings, e.g., $\mathscr{M}(X,D)$ is $D$-valued measurable mappings defined on $X$.

A common trait with the examples worked out in section 6.2 is that each of them can be recast as finding an optimal decision rule in a statistical estimation problem . Furthermore, deep neural networks offer a computationally feasible implementation of estimating this decision rule from supervised training data.

There is a wide range of other tasks beyond those mentioned in section 6.2 that can be represented as a non-randomised decision rule, which in turn is efficiently parametrized by a suitable deep neural network architecture. The (incomplete) list below is based on [4, 77] and aims to show the diversity of tasks from computer vision that can be approached successfully using a suitable deep neural network architecture.

Inpainting:

This is essentially interpolation/extrapolation to recover lost or deteriorated parts of images and videos and approaches based on trainable neural networks [97].

Depixelization/super-resolution:

The task is here to upsample, i.e., to synthesize realistic details into images while enhancing their resolution [17] or to fill out information “between” pixels by increasing the resolution of the final picture, also known as the single image super-resolution problem [79].

Demosaicing:

The task here is to reconstructing a full color image from the incomplete color samples output from an image sensor [90]. Almost all digital cameras, ranging from smartphone cameras to the top-of-the-line digital SLR cameras, use a demosaicing algorithm to convert the captured sensor information into a color image.

Colourising:

The task is to apply color to grey scale photos and videos [45].

Image translation:

The task is to translate between two classes of images of the same object, e.g., computed tomography and magnetic resonance imaging images [95].

Object recognition:

This visual classification task involves localization, detection and classification and this can be seen as an example of constellation models, which are a general class of model that describe objects in terms of a set of parts and their spatial relations [77, section 20.5]. An integrated framework based on deep convolutional networks for detection and localizatio was already introduced in [86], see also [40] for recognition and [26] for scene understanding. The most well-known use case is recognition of multiple faces in an image where statistical shape models play a central role [100, 22, 94] An analogous task relevant for clinical image guided diagnostics is detecting melanoma from images of skin lesions [24, 36].

Non-rigid image registration:

The task here is to deform a template image in a “natural” way so that it matches a target image. This is a key step in spatiotemporal imaging and deep neural networks have been utilized for this purpose [29, 98].

Parametric regression:

The task here is to statistically determine relationships among variables (parameters) in a statistical model. This is important in clinical diagnostics where one seeks to determine risk factors and biomarkers of diagnostic and prognostic value associated with clinical progression and severity of specific diseases. An example of image guided regression is estimating cardiovascular risk factors, such as age, from retinal fundus photographs [76]. Another is predicting patient overall survival as in the 2018 Multimodal Brain Tumor Segmentation Challenge based on BRATS imaging data [68], and predicting scores for Alzheimer’s disease from imaging data [67, 59].

The abstract framework in section 7.1 allows one to perform any of the above tasks (and those in section 6.2) jointly with a reconstruction step for solving an inverse problem. Examples of the latter are deconvolution, Fourier inversion (magnetic resonance imaging imaging), or more elaborate schemes for various types of tomographic imaging. A key success factor is access to suitable training data, another is usage of a differentiable loss function along with a trainable differentiable parametrization (deep neural network architecture) of the task operator.

We will present a generic framework for task adapted reconstruction that is computationally feasible and adaptable to specific inverse problems and tasks. A key element is to formalize both the reconstruction and task as non-randomized decision rules within a statistical estimation problem.

More precisely, our starting point is the inverse problem in section 5 where the data model in eq. 3 is known. Following section 5.2, the reconstruction step can be seen as a decision rule in a statistical estimation problem defined by the statistical model $\bigl((Y,\mathfrak{S}_{Y}),\{\mathcal{M}(x)\}_{x\in X}\bigr)$, decision space $(X,\mathfrak{S}_{X})$, and loss $\ell_{X}\colon X\times X\to\mathbb{R}$. Given a prior $\pi_{0}\in\mathscr{P}_{X}$ allows us to define a reconstruction method as a non-randomized decision rule that minimizes the $\pi_{0}$-average risk (Bayes risk), i.e., as a mapping that solves

Likewise, following section 6.1 a task becomes a decision rule in a statistical estimation problem defined by the statistical model $\bigl((X,\mathfrak{S}_{X}),\{\mathbb{P}_{z}\}_{z\in\triangle}\bigr)$, decision space $(D,\mathfrak{S}_{D})$, and loss given by eq. 9 with known $\featuremap\colon\triangle\to D$ (feature extraction map) and $\ell_{D}\colon D\times D\to\mathbb{R}$ (decision distance). Given a task prior $\eta_{0}\in\mathscr{P}_{\triangle}$, the non-randomized decision rule representing the task (task operator) can be seen as a minimizer to the $\eta_{0}$-average risk (Bayes risk):

We now have the following three approaches to task adapted reconstruction.

Sequential approach:

A sequential approach starts with determining the reconstruction operator, and then uses it to define the the task operator. It is based on assuming that statistical assumptions for $(\mathsf{z},\mathsf{x})\sim\eta_{0}\otimes\mathbb{P}_{z}$ and $(\mathsf{x},\mathsf{y})\sim\pi_{0}\otimes\mathcal{M}(x)$ are consistent, e.g., by assuming that $\mathbb{P}_{z}$ is the push forward of $\mathcal{M}(x)$ through the reconstruction operator²² 2 One may here consider alternate assumptions, like assuming that $\pi_{0}\in\mathscr{P}_{X}$ can be obtained by marginalizing the measure $\eta_{0}\otimes\mathbb{P}_{z}$ over $\triangle$ using $\eta_{0}\in\mathscr{P}_{\triangle}$.. The task adapted reconstruction is then given as

$$\widehat{\mathcal{T}}\circ\hskip 2.0pt\widehat{\hskip-2.0pt\ForwardOp}^{\dagger}\colon Y\to D$$

(18)

where $\hskip 2.0pt\widehat{\hskip-2.0pt\ForwardOp}^{\dagger}\in\mathscr{M}(Y,X)$ solves eq. 16 and $\widehat{\mathcal{T}}\in\mathscr{M}(X,D)$ solves eq. 17.

End-to-end approach:

The fully end-to-end approach ignores the distinction between reconstruction and the task. Assuming $(\mathsf{z},\mathsf{y})\sim\nu$ for some measure $\nu\in\mathscr{P}_{\triangle\times Y}$, the task adapted reconstruction is here given as $\widehat{\EndToEnd}\colon Y\to D$ that solves

$$\widehat{\EndToEnd}\in\argmin_{\EndToEnd\in\mathscr{M}(Y,D)}\biggl\{\Expect_{\eta_{0}\otimes\mathbb{P}_{z}}\Bigl[\ell_{D}\bigl(\featuremap(\mathsf{z}),\EndToEnd(\mathsf{y})\bigr)\Bigr]\biggr\}\quad\text{where $(\mathsf{z},\mathsf{x})\sim\nu$.}$$

(19)

Joint approach:

The joint approach is a middle-way between the sequential and end-to-end approaches. It is based on assuming that there is a joint law $(\mathsf{z},\mathsf{x},\mathsf{y})\sim\sigma$, which by the chain rule in probability can be written in terms of conditional probabilities:

$$\,\mathrm{d}\sigma(z,x,y)=\,\mathrm{d}\pi(y\mid z,x)\,\mathrm{d}\pi(x\mid z)\,\mathrm{d}\pi(z).$$

Assume that $x$ is a sufficient statistic for $y$, i.e., $\,\mathrm{d}\pi(y\mid z,x)=\,\mathrm{d}\pi(y\mid x)$. Combined with $(\mathsf{z},\mathsf{x})\sim\eta_{0}\otimes\mathbb{P}_{z}$ and $(\mathsf{x},\mathsf{y})\sim\pi_{0}\otimes\mathcal{M}(x)$, we obtain

$$\,\mathrm{d}\sigma(z,x,y)=\,\mathrm{d}\mathcal{M}(x)(y)\,\mathrm{d}\mathbb{P}_{z}(x)\,\mathrm{d}\eta_{0}(z).$$

Next, we introduce a joint loss that interpolates between the sequential case and the end-to-end case. Specifically, we let $\ell_{\mathrm{joint}}\colon(X\times D)\times(X\times D)\to\mathbb{R}$ be given as

$$\ell_{\mathrm{joint}}\bigl((x,d),(x^{\prime},d^{\prime})\bigr):=(1-C)\ell_{X}(x,x^{\prime})+C\ell_{D}(d,d^{\prime})\quad\text{for fixed $C\in[0,1]$.}$$

(20)

Task adapted reconstruction is now given by eq. 18 where the operators jointly solve

$$(\hskip 2.0pt\widehat{\hskip-2.0pt\ForwardOp}^{\dagger},\widehat{\mathcal{T}})\in\!\!\argmin_{\begin{subarray}{c}\mathcal{T}\in\mathscr{M}(X,D)\\ \ForwardOp^{\dagger}\in\mathscr{M}(Y,X)\end{subarray}}\Expect_{\sigma}\Bigr[\ell_{\mathrm{joint}}\Bigl(\bigl(\mathsf{x},\featuremap(\mathsf{z})\bigr),\bigl(\ForwardOp^{\dagger}(\mathsf{y}),\mathcal{T}\circ\ForwardOp^{\dagger}(\mathsf{y})\bigr)\Bigr)\Bigr].$$

(21)

Note first that in the limit $C\to 0$, the joint approach becomes the sequential one. Next, it may seem sufficient to only consider the loss $\ell_{D}$ in eq. 21, i.e., to set $C=1$ in eq. 20, which recovers the end-to-end approach. There is however a problem with non-uniqueness in this case since

This non-uniqueness does not arise when $C<1$, so incorporating a loss term associated with the reconstruction may act as a regularizer. This also indicates that the limit $C\to 1$ does not necessarily coincide with the case $C=1$.

In section 5.3.1 we mentioned the difficulty to select an appropriate prior $\pi_{0}\in\mathscr{P}_{X}$ for Bayesian inversion whereas the measure $\mathcal{M}(x)\in\mathscr{P}_{Y}$ is often known by the data model. Furthermore, both measures $\eta_{0}\in\mathscr{P}_{\triangle}$ and $\mathbb{P}_{z}\in\mathscr{P}_{X}$ must be considered as unknown for most tasks. Hence any realistic scenario would contain $\sigma$ as an unknown. An option is to replace these measures by their empirical counterparts given by suitable supervised training data.

Another concern is computational feasibility. The optimizations in eqs. 16, 17, 19 and 21 are taken over all measurable mappings between relevant spaces, which is computationally unfeasible. This can be addressed by considering parametrized sets of measurable mappings as done in sections 5.4 and 6.1. More precisely, we employ a learned iterative scheme to parametrize a family of reconstruction methods $\ForwardOp_{\theta}^{\dagger}\colon Y\to X$ since this parametrization includes knowledge about the data model (section 5.4). Likewise, decision rules associated with the task are given by an appropriate parametrized family of mappings $\mathcal{T}_{\vartheta}\colon X\to D$. Finally, the approach for the end-to-end setting is to directly parametrize $\EndToEnd_{\vartheta}\colon Y\to D$.

Using such parametrizations allows one to reformulate eqs. 16 and 17 as eq. 25, eq. 19 as eq. 25, and eq. 21 as eq. 27. A key aspect for the implementation is to use stochastic gradient descent based methods for finding appropriate parameters by approximately solving the empirical versions of eqs. 25, 19 and 27. This requires that the above parametrizations are differentiable, which in particular requires using differentiable loss-functions.

Depending on the type of supervised training data, we can now pursue either of the three approaches listed in section 7.1.

Sequential approach:

Here we have access to two sets of supervised training data that are coupled:

$$\begin{split}(x_{i},y_{i})\in X\times Y&\quad\text{generated by $(\mathsf{x},\mathsf{y})\sim\pi_{0}\otimes\mathcal{M}(x)$ for $i=1,\ldots,m$}\\ (z_{i},x_{i})\in\triangle\times X&\quad\text{generated by $(\mathsf{z},\mathsf{x})\sim\eta_{0}\otimes\mathbb{P}_{z}$ for $i=1,\ldots,m$.}\end{split}$$

(23)

The coupling is that $x_{i}$’s in the second data set (bottom one) are reconstructions obtained from $y_{i}$’s in the first data set (top one). This ensures consistency with the statistical assumptions mentioned before for the sequential approach. The reconstruction is then given by the mapping

$$\mathcal{T}_{\vartheta^{*}}\circ\ForwardOp_{\theta^{*}}^{\dagger}\colon Y\to D$$

(24)

where $\theta^{*}\in\Theta$ solves eq. 8 and $\vartheta^{*}\in\Xi$ solves eq. 11, i.e.,

$$\begin{split}\theta^{*}&\in\argmin_{\theta\in\Theta}\Bigl\{\frac{1}{m}\sum_{i=1}^{m}\ell_{X}\bigl(x_{i},\ForwardOp_{\theta}^{\dagger}(y_{i})\bigr)\Bigr\}\\ \vartheta^{*}&\in\argmin_{\vartheta\in\Xi}\Bigl\{\frac{1}{m}\sum_{i=1}^{m}\ell_{D}\bigl(\featuremap(z_{i}),\mathcal{T}_{\vartheta}(x_{i})\bigr)\Bigr\}.\end{split}$$

(25)

Note here that the only requirement for the sequential approach is that the assumptions $(\mathsf{z},\mathsf{x})\sim\eta_{0}\otimes\mathbb{P}_{z}$ and $(\mathsf{x},\mathsf{y})\sim\pi_{0}\otimes\mathcal{M}(x)$ are jointly consistent. The learned task operator given by solving for $\vartheta^{*}$ in eq. 25 is only well defined for input taken from the support of it’s training data, so it may fail when applied to data it has never seen. This is especially the case if new data has a different statistical assumption. Hence, it is important to ensure the range of the reconstruction operator is contained in the support of the elements $x\in X$ used to train the task. In most practical implementations, this is ensured by simply letting $x_{i}$’s in $(z_{i},x_{i})$ in eq. 23 be the output of the learned reconstruction operator $\ForwardOp_{\theta^{*}}^{\dagger}\colon X\to Y$.

End-to-end approach:

Supervised training data is here of the form

$$(z_{i},y_{i})\in\triangle\times X\quad\text{generated by $(\mathsf{z},\mathsf{x})\sim\eta_{0}\otimes\mathbb{P}_{z}$ for $i=1,\ldots,m$.}$$

(26)

The reconstruction is given by $\EndToEnd_{\vartheta}\colon Y\to D$ where $\vartheta^{*}\in\Xi$ solves

$$\vartheta^{*}\in\argmin_{\vartheta\in\Xi}\Bigl\{\frac{1}{m}\sum_{i=1}^{m}\ell_{D}\bigl(\featuremap(z_{i}),\EndToEnd_{\vartheta}(y_{i})\bigr)\Bigr\}.$$

(27)

Joint approach:

In this approach we assume access to supervised training data that jointly involves the reconstruction and task:

$$(x_{i},y_{i},z_{i})\in X\times Y\times\triangle\quad\text{generated by $(\mathsf{x},\mathsf{y},\mathsf{z})\sim\sigma$ for $i=1,\ldots,m$.}$$

(28)

Given such data, the corresponding reconstruction method can be defined as in eq. 24 where $(\theta^{*},\vartheta^{*})\in\Theta\times\Xi$ solves the following joint empirical loss minimization:

$$\begin{split}(\theta^{*},\vartheta^{*})&\in\!\!\!\argmin_{(\theta,\vartheta)\in\Theta\times\Xi}\biggl\{\frac{1}{m}\sum_{i=1}^{m}\ell_{\mathrm{joint}}\Bigl(\bigl(x_{i},\featuremap(z_{i})\bigr),\bigl(\ForwardOp_{\theta}^{\dagger}(\mathsf{y}_{i}),\mathcal{T}_{\vartheta}\circ\ForwardOp_{\theta}^{\dagger}(\mathsf{y}_{i})\bigr)\Bigr)\biggr\}\end{split}$$

(29)

where $\ell_{\mathrm{joint}}\colon(X\times D)\times(X\times D)\to\mathbb{R}$ is the joint loss in eq. 20. Note that one may in addition have access to separate sets of training data of the form eq. 23 and eq. 28. In such case, it is possible to first pre-train by solving for (8) and (11) separately, and use the resulting outcomes to initialize an algorithm for solving eq. 29.

In the following we demonstrate performance of the task adapted reconstruction scheme for eq. 24 that is based on solving eq. 29. All cases involve tomographic reconstruction from 2D parallel beam data and as tasks, we consider MNIST classification and segmentation.