Deep learning methods for inverse problems

Direct mapping

The direct mapping category is used as the objective criterion in a large body of research in deep learning based inverse problems (Adler & Öktem, 2017; Antholzer, Haltmeier & Schwab, 2019; Jin et al., 2017; Kelly, Matthews & Anastasio, 2017; Anirudh et al., 2018; Zhang & Ghanem, 2018; Fan et al., 2017). These methods seek to find end-to-end solutions for

(14) $$\underset{W_1}{min}\left\{\sum _{i=1}^{N}D(\underset{\_}{z},G(\underset{\_}{m},W_1))+\lambda R(G(\underset{\_}{m},W_1))\right\}$$whereby D(·, ·) is the cost function to be minimized by a DNN $G(\underset{\_}{m},W_1)$, on the basis of optimizing DNN weights W₁. $R(G(\underset{\_}{m},W_1))$ specifies a generic analytical regularizer, to restrict the estimator to feasible solutions.

The Direct Mapping category approximates an estimator G as an inverse to the forward model F, requiring a dataset of pairs $\{({\underset{\_}{m}}_i,{\underset{\_}{z}}_i){\}}_i$ of observed measurements and corresponding target system parameters, as illustrated in Fig. 1.

This category of DNN is typically used in those cases where we have a model-based imaging system having a linear forward model $\underset{\_}{m}=F(\underset{\_}{z})$, where $\underset{\_}{z}$ is an image, so that convolution networks (CNNs) are nearly always used. As discussed earlier, for Image Restoration problems the measurements themselves are already images, however in more general contexts we may choose to project the measurements as $F^H\underset{\_}{m}$, back into the domain of $\underset{\_}{z}$, such that the CNN is trained to learn the estimator

(15) $$\underset{\_}{\hat{z}}=G(F^H\underset{\_}{m},W_1)$$

The translation invariance of F^H F, relatively common in imaging inverse problems, makes the convolutional-kernel nature of CNNs particularly suitable for serving as the estimator for these problems.

In general, the performance of direct inversion is remarkable (Lucas et al., 2018). However the receptive field (i.e., the size of the field of view the unit has over its input layer) of the CNN should be matched to the support of the point spread function (Aggarwal, Mani & Jacob, 2018). Therefore, large CNNs with many parameters and accordingly extensive amount of training time and data are often needed for the methods in this category. These DNNs are highly problem dependent and for different forward models (e.g., with different matrix sizes, resolutions, etc.) a new DNN will need to be learned.

Data consistency optimizer

The Data Consistency Optimizer category of deep learning aims to optimize data consistency as an unsupervised criterion within a variational framework (Aggarwal, Mani & Jacob, 2018; Cha, Lee & Oh, 2019):

(16) $$\underset{W_2}{min}\left\{\sum _{i=1}^{N}D(\underset{\_}{m},F(G(\underset{\_}{m},W_2)))+\lambda R(G(\underset{\_}{m},W_2))\right\}$$where, as in (14), D(·, ·) is the cost function to be minimized by DNN $G(\underset{\_}{m},W_2)$, parameterized by weights W₂, subject to regularizer $R(G(\underset{\_}{m},W_1))$. The overall picture is summarized in Fig. 2.

In contrast to (14), where the network cost function D is expressed in the space of unknowns $\underset{\_}{z}$, here (16) expresses the cost in the space of measurements $\underset{\_}{m}$, based on forward model F(·). That is, the data consistency term is no longer learning from supervised examples, rather from the forward model we obtain an unsupervised data consistency term, not needing data labels, whereby the forward model provides some form of implicit supervision.

Compared to the direct mapping category, the use of the forward model in (16) leads to a network with relatively few parameters, in part because the receptive field of the DNN need not be matched to the support of the point spread function. However, the ill-posedness of the inverse problem causes a semi-convergent behaviour (Arridge et al., 2019) using this criterion, therefore an early stopping regularization needs to be adopted in the learning process.

Deep regularizer

Finally, the Deep Regularizer (DR) category needs a different problem modeling scheme, since there is not a learning phase as in DM and DC. Instead, only a DNN (usually a classifier) is trained to be used as the regularizer in a variational optimization framework. That is, DR continues to optimize the data consistency term, however the overall optimization process is undertaken in the form of an analytical variational framework and uses a DNN as the regularizer (Chang et al., 2017; Li et al., 2018):

(17) $$\underset{\underset{\_}{\hat{z}}}{min}\left\{\sum _{i=1}^{N}D(\underset{\_}{m},F(\underset{\_}{\hat{z}}))+\lambda R(\underset{\_}{\hat{z}},W_3)\right\}$$

Here $R(\underset{\_}{\hat{z}},W_3)$ is a pre-trained deep regularizer, based on weight matrix W₃, usually chosen as a deep classifier (Chang et al., 2017; Li et al., 2018), discriminating the feasible solutions from non-feasible ones.

This category usually includes an analytical variational framework consisting of a data consistency term and a learned DNN to capture redundancy in the parameter space (see Fig. 3).

Since our interest is in the DNN solution of the inverse problem, and not the details of the optimization, we have chosen two fairly standard optimization approaches, a simplex/Nelder-Mead approach Singer & Nelder (2009) (DR-NM) and a Genetic Algorithm strategy (DR-GA), both based on their respective Matlab implementations. Because GA solutions may be different from one run to the next, in general we report the results averaged over multiple independent runs.

The Deep Regularizer category needs the fewest parameter settings, compared to the earlier categories; however because of the optimization based inference step it is computationally demanding.

Linear regression

We begin with an exceptionally simple inverse problem. Consider a set of one dimensional samples $\{(x^{(i)},m_y^{(i)}){\}}_{i=1}^N$, subject to noise, with some number of the training data subject to more extreme outliers, as illustrated in Fig. 4.

As an inverse problem, we need to define the forward model, which for linear regression is simply

(18) $${\underset{\_}{m}}_y=\alpha \underset{\_}{x}+\beta +\underset{\_}{\nu }.$$

Since our interest is in assessing the robustness of the resulting inverse solver, the number and behaviour of outliers should be quite irregular, to make it challenging for a network to generalize from the training data. As a result, the noise $\underset{\_}{\nu }$ is random variance, plus heavy-tailed (power law) outliers, where the number of outliers is exponentially distributed.

For this inverse problem, the unknown state is comprised of the system parameters ${\underset{\_}{z}}^T=[\alpha ,\beta ]$. Thus linear regression leads to a reconstruction problem, for which the goal is to recover the line parameters from a sample set including noisy and outlier data points.

With the problem defined, we next need to formulate an approach for each of the three solution categories. For direct mapping (DM) and data consistency (DC), the training data and DNN structures are the same, shown in Fig. 5, where the DC approach includes an additional layer which applies the given forward model of (18). We used the KERAS library, in which a Lambda layer is designed for this forward operation.

Since the problem is one-dimensional with limited spatial structure, the network contains only dense feed-forward layers. Residual blocks are used in order to allow gradient flow through the DNN and to improve training. Network training was based on 1,000 records, each of N = 500 noisy sample points.

The Deep Regularizer (DR) category needs a different problem modeling scheme, since 2there is not a learning phase as in DM and DC. Instead, only a DNN (usually a classifier) is trained to be used as the regularizer in a variational optimization framework. The DNN regularizer is given the system parameters (α, β) and determines whether they account for a feasible line. Here, we define the feasible line as a line having a tangent in some specified range. We generate a synthetic set of system parameters with associated labels for training a fully connected DNN as the regularizer for this category. Since our interest is in the DNN solution of the inverse problem, and not the details of the optimization, we have chosen two fairly standard optimization approaches, a simplex/Nelder-Mead approach Singer & Nelder (2009) and a Genetic Algorithm (GA) strategy, both based on their respective Matlab implementations. Because GA solutions may be different over multiple runs, we report the results averaged over ten independent runs.

Table 2 shows the average solution found by each category over 10 independent trainings for DM and DC, and 10 independent inferences for DR. The table also reports Least-Squares (LS) results as a point of reference method, particularly to show the improvement that deep learning methods have to offer for robustness in solving inverse problems. Observe the significant difference when the DNN methods are trained with noise-free as opposed to noisy data, such that the noisy training data force the network to acquire a robustness to outliers.

For DR we trained a 5 layer MLP with dense layers of sizes 5, 4, 3, 2, 1, as the regularizer, using the generated synthetic data including feasible line parameters (in the specific range) as the positive training samples and invalid line parameters as the negative training samples. The average test accuracy of the trained regularizer is 95.70%.

We performed the Wilcoxon signed rank test, for both cases of training with noisy data (Table 3) and noise-free training (Table 4). The tables show the pairwise p-values over the 10 independent runs. A p-value in excess of 0.05 implies that the two methods are likely to stem from the same distribution; in particular, the Wilcoxon test computes the probability that the difference between the results of two methods are from a distribution with median equal to zero. Clearly all of the DNN methods are statistically significantly different from the least-squares (LS) results. For noisy training data, the statistical results in Table 3 show similar performance for DM and DC, and for DR-NM and DR-GA, the latter similarity suggesting that the specific choice of optimization methodology does not significantly affect the DR performance.

The results in Table 2 show that DM and DC significantly improve in robustness when trained with noisy data, relative to training with noise-free data. The principal difference between DM/DC vs DR is the learning phase for DM/DC, allowing us to conclude that, at least for reconstruction problems, a learning phase using noisy samples in training significantly improves the robustness of the solution. A further observation is that whereas DM and DC achieve similar performance, DC is unsupervised and DM is supervised. Thus it would appear that the forward model knowledge and the data consistency term as objective criterion for DC provide an equal degree of robustness compared to the supervised learning in DM.

For this reconstruction problem, we conclude that both DC and DM perform well, with the unsupervised DC showing strong performance both with noisy and noise-free training data.

Finding coefficients of Burgers’ PDE (inverse problems in PDEs)

To test deep learning for Inverse Problems in PDEs, we chose Burgers’ Partial Differential Equations (PDEs) as a dynamic, continuous time PDE in our experiments.

Burgers’ PDE or Bateman–Burgers equation (Bateman, 1915; Burgers, 1948) is a basic partial differential equation occurring in various areas of applied mathematics, such as fluid mechanics, nonlinear acoustics, gas dynamics, and traffic flow Xin (2009). This setup encapsulates a wide range of problems in the mathematical physics including conservation laws, diffusion processes, advection-diffusion-reaction systems, and kinetic equations (Raissi, Perdikaris & Karniadakis, 2019).

For a given field $\underset{\_}{m}(x,t)$, we consider Burgers’ equation as the forward function, defined as

(19) $$F(\underset{\_}{m}(x,t);\underset{\_}{z})={\underset{\_}{m}}_t+z_1\underset{\_}{m}{\underset{\_}{m}}_x-z_2{\underset{\_}{m}}_{xx}=0x\in \mathrm{\Omega },t\in [0,T]$$where z₁, z₂ denote the parameters of the equation and $\underset{\_}{m}(x,t)$ the state of the system, the subscripts denoting partial differentiation in either time or space. The goal is then to estimate parameters z₁, z₂, given a collection of points (Raissi, Perdikaris & Karniadakis, 2019; Bu & Karpatne, 2021).

Solutions for this type of inverse problems in the literature, including Physics Informed Neural Networks (PINNs) (Raissi, Perdikaris & Karniadakis, 2019), and Quadratic Residual Neural Networks (QRes) (Bu & Karpatne, 2021) actually use a regularized version of data loss, with $F(\underset{\_}{m};\underset{\_}{z}{)}^2$ as the regularizer. For instance, PINNs are defined as

(20) $$G:=F(\underset{\_}{m}(x,t);\underset{\_}{z})$$with F from (19), and using a DNN to approximate $\underset{\_}{m}(x,t)$. The DNN, along with (20), form the Physics Informed Neural Network G(t, x) in which the chain rule could be used for differentiating compositions of functions using automatic differentiation (Baydin et al., 2018), which we call AutoGrad, and has the same parameters as the network representing G(t, x; W₂), albeit with different activation functions due to the action of F. The shared parameters between the neural networks m(t, x) and G(t, x; W₂) can be learned by minimizing the mean squared error loss

(21) $$J_{DM}=PGLoss+DataLoss$$where

(22) $$DataLoss={\displaystyle \frac{1}{N_m}\sum _{i=1}^{N_m}|\underset{\_}{m}(t_m^i,x_m^i;W_1)-{\underset{\_}{m}}^i{|}^2}$$and

(23) $$PGLoss={\displaystyle \frac{1}{N_G}\sum _{i=1}^{N_G}|G(t_G^i,x_G^i;W_2){|}^2}$$where $\{t_m^i,x_m^i,{\underset{\_}{m}}^i\}$ denote the initial and boundary training data on $\underset{\_}{m}(t,x)$ and $\{t_G^i,x_G^i{\}}_{i=1}^{N_G}$ specify the collocations points for G(t, x). DataLoss corresponds to the initial and boundary data while PGLoss enforces the structure imposed by (19) at a finite set of collocation points.

This definition of the loss functions makes them consistent with the objective of the Data Consistency optimizer (DC) solution category. Therefore, we include these methods within DC in our experiments.

In the case of the DM approach, we can define the loss function as

(24) $$J_{DM}(\underset{\_}{m}(x,t);\underset{\_}{z})={\displaystyle \frac{1}{N_G}\sum _{i=1}^{N_G}|z_{pred}^i-z^{GT}{|}^2+|G(t_G^i,x_G^i;W_2){|}^2}$$where zⁱ_pred, z^GT stand for the predicted parameter by the solution category and its ground truth, respectively.

Figure 6 shows the DM and DC DNNs for Burgers’ Inverse problem, where AutoGrad, the automatic differentiation component, is used for computing the needed gradients. The DC solution category only uses the PGLoss in its training procedure.

For DR, we define the loss function as

(25) $$J_{DR}={\displaystyle \frac{1}{N_G}\sum _{i=1}^{N_G}|G(t_G^i,x_G^i){|}^2+\lambda R(\underset{\_}{m};W_3)}$$where ${\displaystyle \frac{1}{N_G}{\sum }_{i=1}^{N_G}|G(t_G^i,x_G^i){|}^2}$ is the data consistency term and $R(\underset{\_}{m};W_3)$ is a deep classifier, for which we trained an MLP classifier with dense layers of size 5, 4, 3, 2, 1, trained by he available measurement states, to control the values of m to be in the specified range, provided by Bu & Karpatne (2021).

For the experiments, we used synthetic data provided by Bu & Karpatne (2021) as the training and test data, where the standard deviation of the noise is set to 1% of the data standard deviation, and for xⁱ_m, tⁱ_m we used equi-spaced values in the specified ranges. In the case of outliers, we used additive Gaussian noise with magnitude equal to 10 times the data standard deviation for 0.05% of the data. Table 5 compares the MSE between the obtained parameter values by existing methods, averaged over 5 independent training/ inferences.

The statistical analyses of the results in Table 5 are reported in Table 6. From Tables 5 and 6, it is observable that in this case it is DC which achieves the best robustness performance. The statistical analysis shows that the choice of DR optimization method does not impact the results. The results also show that the learning phase in DC significantly improves the obtained results compared with DR under the same objective.

Image denoising (restoration)

We now consider an image denoising problem, following the steps described in “Linear Regression” for regression. We consider real and synthetic images, including 5 classes and 1,200 training images, 400 test images per class, from the Linnaeus dataset (Chaladze & Kalatozishvili, 2017) as real data, and synthesized 5,000 texture images generated by sampling from stationary periodic kernels, as synthetic data.

The synthetic images are generated using an FFT method (Fieguth, 2010), based on a thin-plate second-order Gauss–Markov random field kernel

(26) $$𝒫=\left[\begin{array}{ccccc}0& 0& 1& 0& 0\\ 0& 2& -8& 2& 0\\ 1& -8& 20+{\alpha }^2& -8& 1\\ 0& 2& -8& 2& 0\\ 0& 0& 1& 0& 0\\ & & & & \end{array}\right]$$such that a texture T is found by inverting the kernel in the frequency domain,

(27) $$T=FF{T}_2^{-1}\left(\sqrt{1\oslash FF{T}_2(𝒫)}\odot FF{T}_2(W)\right),$$with $\odot$, $\oslash$ as element-by-element multiplication and division, W as unit-variance white noise, and with the kernel $P$ zero-padded to the intended size of T. Further details about this approach can be found in Fieguth (2010).

Parameter α², affecting the central element of the kernel $P$, effectively determines the texture spatial correlation-length in T, as

(28) $${\alpha }^2={10}^{4-lo{g}_{10}u}$$for process correlation length, u, measured in pixels. We set u to be a random integer in the range [10, 200] in our experiments.

All images are set to be 64 × 64 in size, with pixel values normalized to [0, 1]. Pixels are corrupted by additive Gaussian noise, with an exponentially distributed number of outliers.

The inverse problem is a restoration problem, having the objective of restoring the original image from its noisy/outlier observation. The linear forward model is

(29) $$\underset{\_}{m}=\underset{\_}{z}+\underset{\_}{\nu }$$for measured, original, and added noise, respectively. The Gaussian noise ν has zero mean and random variance, and an exponential number of pixels become outliers, their values replaced with a uniformly distributed random intensity value.

We used 5,000 training samples and 500 test samples for the learning and evaluation phases of the DM and DC approaches. The DNN structure for both DM and DC is the same and is shown in Fig. 7. In the case of DC, we design a DNN layer to compute the forward function. Since we are dealing with input images, both as measurements and system state, we design a fully convolutional DNN in an encoder-decoder structure, finding the main structures in the image through encoding and recovering the image via decoding. Since there may be information loss during encoding, we introduce skip connections to help preserve desirable information.

The DR category needs a pre-trained regularizer which determines whether the prediction is a feasible texture image. We trained a classifier for texture discrimination, generated using (27), from ordinary images gathered from the web, as the regularizer. Both GA and Nelder-Mead optimizers are used.

We use peak signal to noise ratio (PSNR) as the evaluation criterion, computed as

(30) $$PSNR(I^{pred},I^{GT})=20\cdot lo{g}_{10}max(I^{pred})-10\cdot lo{g}_{10}MSE,$$

(31) $$MSE={\displaystyle \frac{1}{n}\sum _{i,j}(I_{i,j}^{GT}-I_{i,j}^{pred}{)}^2}$$where I^GT_i,j, I^pred_i,j are the (i,j)^th pixel in the ground-truth and predicted images, respectively. Note that in the DR case, since the input and output of the model are 64 * 64 = 4,096 images, the GA optimization routine was unable to find the solution in a reasonable time, therefore we do not avoid report any DR-GA results for this problem.

As a reference point, we also report results obtained by the non-local means (NLM) filter (Buades, Coll & Morel, 2011), to give insight into the amount of improvement of deep learning inverse methods over a well-established standard in image denoising.

Figure 8 shows results based on synthetic textures. Each row in the figure shows a sample image associated with a particular correlation length noise standard deviation. The DM approach offers by far the best reconstruction among the DNN methods, and outperforms NLM in terms of PSNR. The time complexity of GA in DR-GA makes it inapplicable to problems of significant size (even though the images were still quite modest in size).

The Wilcoxon signed rank test was performed on the DM, DC and DR-(Nelder-Mead) results. The statistical analysis of the obtained results gave a p value of 0.002 for each pairwise comparison, implying a statistically significant difference, thus the very strong performance of DM in Fig. 8 is validated.

In the case of real images, Fig. 9 shows the visual results obtained by DM, DC and DR-NM for seven test samples.

The statistical analysis is consistent with the results from the synthetic texture case, which is that all pairwise Wilcoxon tests led to a conclusion of statistically significant differences, with p values well below 0.05.

From the results in Figs. 8 and 9 and their respective statistical analyses, we conclude that:

• For image denoising as a prototype for restoration problems, which have the same measurement and system parameter spaces, the concentration of the loss function on the true parameters (as in DM) provides better information and leads to a more effective estimator having greater robustness than the measurements themselves (as in DC).

• DR-(Nelder-Mead) performed poorly, even though it optimizes data consistency, like DC, however we believe that the learning phase in DC, compared to DR, provides knowledge for its inference and allows DC to be more robust than DR for restoration inverse problems.

3D shape inverse rendering (reconstruction)

We now wish to test a 3D shape inverse rendering (IR) (Aldrian & Smith, 2012) problem, for which a 3D morphable model (3DMM) (Blanz & Vetter, 1999) describes the 3D shape of a human face $\underset{\_}{s}$. This model is based on extracting eigenfaces ${\underset{\_}{s}}_i$, usually using PCA, from a set of 3D face shapes as the training data, then to obtain new faces as a weighted combination z_i of the eigenfaces. The 3D shape model reconstructs a 3D face in homogeneous coordinates as

(32) $$\underset{\_}{s}=\underset{\_}{\overline{s}}+\sum _{i=1}^n{z}_i{\underset{\_}{s}}_i,$$where $\underset{\_}{\overline{s}}$ is the mean shape of the 3DMM, and z_i the weight of eigenface ${\underset{\_}{s}}_i$. We use the Besel Face Model (Aldrian & Smith, 2012) as the 3DMM in this experiment for which there are N = 54,390 3D points in each face shape and 199 eigenfaces. We can therefore rewrite (32) as

(33) $${\underset{\_}{s}}_N={\underset{\_}{\overline{s}}}_N+{\underset{\_}{z}}^T\ast S_N$$where S is the tensor of 199 eigenfaces. In our experiments each face is characterized by 72 standard landmarks, shown in Fig. 10, which are normalized and then presented to the system as the measurements. Therefore we actually only care about L = 72 out of N = 54,390 3D points in the 3DMM. This experiment tackles the reconstruction of a 3D human face by finding the weights $\underset{\_}{z}$ of the 3DMM from its input 2D landmarks. We generated training data from the 3DMM by assigning random values to the 3DMM weights, resulting in a 3D human face, and rendered the obtained 3D shape into a 2D image using orthographic projection.

The measurement noise consists of small perturbations of the 2D landmarks, with outliers as much larger landmark perturbations. We add zero-mean Gaussian noise having a standard deviation of 3 × 10³ in the training data and 5 × 10³ in the test data. Outliers are much larger, with a standard deviation of 5 × 10⁴ added to 10 of the 72 landmarks in 10% of the training data and 20% of the test data. Landmark point coordinates are in the range [−8 × 10⁴, 8 × 10⁴], so the outlier magnitudes are very large.

Let subscript _L represent the the set of landmark point indices, in which case the forward model is the orthographic projection

(34) $$\underset{\_}{m}=C{\underset{\_}{s}}_L+\underset{\_}{\nu }C=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\end{array}\right]$$such that C converts from homogeneous 3D to homogeneous 2D coordinates, and the measurement noise is

(35) $$\underset{\_}{\nu }\sim 0.9N(0,3\times {10}^{3}I)+0.1N(0,5\times {10}^{4}I)$$as noise and outliers associated with the projection operator. Since the goal of this inverse problem is to estimate $\underset{\_}{z}$ in the 3DMM for a given 3D shape, we write (34) as

(36) $$\underset{\_}{m}=C({\underset{\_}{\overline{s}}}_L+S_L\ast \underset{\_}{z})+\underset{\_}{\nu }$$

For the DM and DC solutions we generated 4,000 sample faces as training data, using the Besel face model (Aldrian & Smith, 2012) as the 3DMM. The DR regularizer is a pre-trained classifier which discriminates a feasible 3D shape from random distorted versions of it.

In DC we implemented the forward function layer as described in Aldrian & Smith (2012), with the resulting DM and DC DNN shown in Fig. 11, where we used feed-forward layers because the system input is the vectorized 72 2D homogeneous coordinates and its output a weight vector. We design an encoder-decoder structure for DNNs, so as to map the 2D coordinates to a low dimensional space and to recover the parameters from that low dimensional representation.

For the DR regularizer we trained a five layer MLP classifier to discriminate between a 3D face shape, generated by BFM, and randomly generated 3D point clouds as negative examples.

Figure 12 shows visual results obtained by each solution category, where heat maps visualize the point-wise error magnitude relative to the ground truth. The visual results show that the DM and DC methods can capture the main features in the face (including eye, nose, mouth) better than the DR variants, however the differences between DM and DC seem to be negligible.

To validate our observations, the numerical results and respective statistical analyses are shown in Tables 7 and 8. Table 7 lists the RMSE values for each solution category. We used 10 out of sample faces in the BFM model as test cases for reporting the results. In the case of DR (Nelder-Mead) we set the start point, i.e., z₀, as a random value and report the averaged result over 10 independent runs. Note that the RMSE values are expected to be relatively large, since each 3D face shape provided by BFM is a point cloud of 53,490 3D coordinates in the range [−8 × 10⁴, 8 × 10⁴]. As a point of comparison, we computed the average RMSE between a set of 500 generated 3D faces and 1,000 random generated faces, to have a sense of RMSE normalization to random prediction. The average RMSE for random prediction is 1.28 × 10⁴, a factor of two to four times larger than the RMSE values reported in Table 7.

Table 8 shows the results of the Wilcoxon p values for statistical significance in the difference between reported values in Table 7, where we consider a p value threshold of 0.07.

Based on the preceding numerical results and statistical analysis, we claim the following about each solution category facing with Reconstruction inverse problems:

• Broadly, for training and test data not involving outliers, the overall performance of the methods is similar, with DM outperforming. This observation shows that the learning phase is not crucial in the presence of noise, and methods which concentrate on the test data can achieve equal performance compared to trainable frameworks.

• In cases involving outliers the performance of the methods is more distinct, but with the DM and DC methods, having a learning phase for optimizing their main objective term, outperforming the DR variants. We conclude that a learning phase is important to make methods robust to outliers.

• In the case of DR, the results show similar performance of the GA and NM optimization schemes, with GA outperforming NM. This observation encourages the reader to use optimization methods with more exploration power Eftimov & Korošec (2019), the ability of an optimization method to search broadly across the whole solution space, for DR solutions to reconstruction problems.

• In all cases, we can observe that although DC is unsupervised, its performance when solving reconstruction inverse problems is near to that of DM, even outperforming DM in the case of outliers. Therefore, it is possible to solve reconstruction problems even without label information in the training phase.

• One interesting observation is that while 3D shape inverse rendering is a complex reconstruction problem, the results for each solution category are qualitatively similar to the very different and far simpler inverse problem of linear regression, where DC similarly outperformed training data containing noisy and outlier samples.

Single object tracking (dynamic estimation)

Up to this point we have investigated deep learning approaches applied to static problems. We would now like to examine a dynamic inverse problem, that of single-object tracking.

The classical approach for tracking is the Kalman Filter (KF) (Fieguth, 2010) and its many variations, all based on a predictor-corrector framework, meaning that the filter alternates between prediction (asserting the time-dynamics) and correcting (asserting information based on the measurements). For the inverse problem under study, we consider the current location estimation (filtering) in a two dimensional environment. Synthetic object tracking problems, as considered here, are studied in a variety of object tracking papers (Kim et al., 2019; Choi & Christensen, 2013; Black, Ellis & Rosin, 2003; Lyons & Benjamin, 2009), where the specific tracking problem in this section is inspired from the approach of Fraccaro et al. (2017), Vermaak, Lawrence & Pérez (2003).

The inverse problem task is to estimate the current ball location, given the noisy measurement in the corresponding time step and the previous state of the ball. Formally, we denote the measured ball location by ${\underset{\_}{m}}^t$, and the system state, the current location of the ball, as ${\underset{\_}{z}}^t$. The graphical model in Fig. 13 illustrates the problem definition of the tracking problem, where the objective of the inverse problem is to address the dashed line, the inference of system state from corresponding measurement.

To perform the experiments, we generate the training and test sets similar to Fraccaro et al. (2017) except that we assume that our measurements are received from a detection algorithm, which detects the ball location from input images having a size of 32 × 32 pixels, and that the movement of the ball is non-linear.

In each training and test sequence the ball starts from a random location in the 2D environment, with a random speed and direction, and then moving for 30 time steps. The dynamic of the generated data includes changing the ball location ${\underset{\_}{z}}^t$ and its velocity ${\underset{\_}{v}}^t$ as

(37) $${\underset{\_}{z}}^t={\underset{\_}{v}}^{(t-1)}\mathrm{\Delta }t+{\underset{\_}{z}}^{(t-1)}$$

(38) $${\underset{\_}{v}}^t={\underset{\_}{v}}^{(t-1)}-(c({\underset{\_}{v}}^{(t-1)}{)}^2\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}({\underset{\_}{v}}^t))$$where c is a constant and is set to 0.001. In our data, collisions with walls are fully elastic and the velocity decreases exponentially over time. In this simulation, the training and testing data-sets contain 10,000 and 3,000 sequences of 30-time steps, respectively.

The training measurement noise is

(39) $$\underset{\_}{\nu }\sim 0.95N(0,0.2I)+0.05N(0,10I),$$a mixture model of Gaussian noise with 5% outliers. The testing noise is similar,

(40) $$\underset{\_}{\nu }\sim 0.85N(0,0.4I)+0.15N(0,10I)$$with a higher likelihood of outliers.

The inverse problem is single-target tracking for which the dynamic of the model is unknown. The inverse problem of interest is to find ${\underset{\_}{z}}^t$ in

(41) $${\underset{\_}{z}}^t=G({\underset{\_}{z}}^{(t-1)},m_t)$$

As shown in Fig. 13, we can model our problem as a first order Markov model where the current measurement is independent of others given the current system state. The forward model is then defined as

(42) $${\underset{\_}{m}}^t=F({\underset{\_}{z}}^t)=C{\underset{\_}{z}}^t+\underset{\_}{\nu },C=I,\underset{\_}{\nu }\sim N(0,\sigma )$$

We can model Markov models using Recurrent Neural Networks (RNN) (Krishnan, Shalit & Sontag, 2017; Hafner et al., 2019; Rangapuram et al., 2018; Coskun et al., 2017). The DNN structure for DM and DC solution categories is shown in Fig. 14, in which the LSTM layers lead the learning process to capture the time state and dynamic information in the data sequences.

We design the regularizer of the DR category as a classifier to classify location feasibility—those locations lying within the border of the 2D environment. Figure 15 shows the positive and negative samples which we used to train the DR regularizer.

As before, we used GA (DR-GA) and Nelder-Mead (DR-NM) algorithms as optimizers for DR. In the case of using Nelder-Mead, the results vary as a function of starting point ${\underset{\_}{z}}_0$, and found that using the last sequence measurement as the starting point empirically gave the best result for DR-NM.