arXiv:1710.10898v1 [cs.CV] 30 Oct 2017

Learning to solve inverse problems
using Wasserstein loss

Jonas Adler1,2 , Axel Ringh1 , Ozan Öktem1 , Johan Karlsson1
1
Department of Mathematics, KTH Royal Institute of Technology
2
Elekta Instrument AB
{ jonasadl, aringh, ozan } @kth.se, johan.karlsson@math.kth.se

Abstract
We propose using the Wasserstein loss for training in inverse problems. In particular,
we consider a learned primal-dual reconstruction scheme for ill-posed inverse
problems using the Wasserstein distance as loss function in the learning. This is
motivated by miss-alignments in training data, which when using standard mean
squared error loss could severely degrade reconstruction quality. We prove that
training with the Wasserstein loss gives a reconstruction operator that correctly
compensates for miss-alignments in certain cases, whereas training with the mean
squared error gives a smeared reconstruction. Moreover, we demonstrate these
effects by training a reconstruction algorithm using both mean squared error and
optimal transport loss for a problem in computerized tomography.

1

Introduction

In inverse problems the goal is to determine model parameters from indirect noisy observations.
Example of such problems arise in many different fields in science and engineering, e.g., in X-ray
computed tomography (CT) [27], electron tomography [28], and magnetic resonance imaging [8].
Machine learning has recently also been applied in this area, especially in imaging applications.
Using supervised machine-learning to solve inverse problems in imaging requires training data where
ground truth images are paired with corresponding noisy indirect observations. The learning provides
a mapping that associates observations to corresponding images. However, in several applications
there are difficulties in obtaining the ground truth, e.g., in many cases it may have undergone a
distortion. For example, a recent study showed that MRI images may be distorted by up to 4 mm
due to, e.g., inhomogeneities in the main magnetic field [36]. If these images are used for training,
the learned MRI reconstruction will suffer in quality. Similar geometric inaccuracies arise in several
other imaging modalities, such as Cone Beam CT and full waveform inversion in seismic imaging.
This work seeks to provide a scheme for learning a reconstruction scheme for an ill-posed inverse
problem with a Wasserstein loss by leveraging upon recent advances in efficient solutions of optimal
transport [10, 22] and learned iterative schemes for inverse problems [3]. The proposed method
is demonstrated on a computed tomography example, where we show a significant improvement
compared to training the same network using mean squared error loss. In particular, using the
Wasserstein loss instead of standard mean squared error gives a result that is more robust against
potential miss-alignment in training data.

2
2.1

Background
Inverse problems

In inverse problems the goal is to reconstruct an estimate of the signal ftrue ∈ X from noisy indirect
measurements (data) g ∈ Y assuming
g = T (ftrue ) + δg.

(1)

In the above X and Y are referred to as the reconstruction and data space, respectively. Both are
typically Hilbert or Banach spaces. Moreover T : X → Y denotes the forward operator, which
models how a given signal gives rise to data in absence of noise. Finally, δg ∈ Y is the noise
component of data. Many inverse problems of interest are ill-posed, meaning that there is no uniques
solution to (1) and hence there is no inverse to T . Typically reconstructions of ftrue are sensitive to
the data and small errors gets amplified. One way to mitigate these effects is to use regularization
[12].
Variational regularization In variational regularization one formulates the reconstruction problem as an optimization problem. To this end, one introduces a data discrepancy functional
f 7→ L(T (f ), g), where L : Y × Y → R, that quantifies the miss-fit in data space, and a regularization functional S : X → R that encodes a priori information about ftrue by penalizing undesirable
solutions. For a given g ∈ Y , this gives an optimization problem of the form
min L(T (f ), g) + λS(f ).

f ∈X

(2)

Here, λ acts as a trade-off parameter between the data discrepancy and regularization functional. In
many cases L is taken to be the negative data log-likelihood, e.g., L(T (f ), g) = k T (f ) − gk22 in
the case of additive white Gaussian noise. Moreover, a typical choice for regularization functional
is total variation (TV) regularization, S(f ) = k∇f k1 [33]. These regularizers typically give rise
to large scale and non-differentiable optimization problems, which requires advanced optimization
algorithms.
Learning for inverse problems In many applications, and so also in inverse problems, data driven
approaches have shown dramatic improvements over the state-of-the-art [24]. Using supervised
learning to solve an inverse problem amounts to finding a parametrized operator T †Θ : Y → X where
the parameters Θ are selected so that
g = T (ftrue ) + δg =⇒ T †Θ (g) ≈ ftrue .
For inverse problems in image processing, such as denoising and deblurring, we have Y = X and
it is possible to apply a wide range of widely studied machine learning techniques, such as fully
convolutional deep neural networks with various architectures, including fully convolutional networks
[18] and denoising auto-encoders [38].
However, in more complicated inverse problems as in tomography, the data and reconstruction spaces
are very different, e.g., their dimension after discretization may differ. For this reason, learning a
mapping from Y to X becomes nontrivial, and classical architectures that map, e.g., images to images
using convolutional networks cannot be applied as-is. One solution is to use fully-connected layers as
in [30] for very small scale tomographic reconstruction problems. A major disadvantage with such a
fully learned approach is that the parameters space has to be very high dimensional in order to be
able to learn both the prior and the data model, which often renders it infeasible due to training time
and lack of training data.
A more successful approach is to first apply some crude reconstruction operator T † : Y → X and then
use machine learning to post process the result. This separates the learning from the complications
of mapping between spaces since the operator T † can be applied off-line, prior to training. Such
an approach has been demonstrated for tomographic reconstruction in [31, 37]. Its drawback for ill
posed inverse problems is that information is typically lost by using T † , and this information cannot
be recovered by post processing.
Finally, another approach is to incorporate the forward operator T and its adjoint T ∗ into the neural
network. In these learned iterative schemes, classical neural networks are interlaced with applications
2

of the forward and backward operator, thus allowing for the learned reconstruction operator to work
directly from data without having to learn the data model. For example, in [39] an alternating direction
method of multipliers (ADMM)-like scheme for Fourier inversion is learned and [32] consider solving
inverse problems typically arising in image restoration by a learned gradient-descent scheme. In [4]
this later approach is shown to be applicable to large scale tomographic inversion. Finally, in [3] they
apply learning in both spaces X and Y , yielding a Learned Primal-Dual scheme, and show that it
outperforms learned post-processing for reconstruction of medical CT images.
Loss functions for learning Once the Θ parametrization of T †Θ is set, the parameters are typically
chosen by minimization of some loss functional L. Without doubt, the most common loss function is
the mean squared error, also called L2 loss, given by
h
i
L(Θ) = Ef,g k T †Θ (g) − fk22 .
(3)
It has however been noted that it is sub-optimal for imaging, and a range of other loss functions have
been investigated. These include the classical `p norms and the structural similarity index (SSIM)
[40], as well as more complex losses such as perceptual losses [20] and adversarial networks [26].
Recently, optimal mass transport has also been considered as loss function for classification [14] and
generative models [5]. In this work we consider using optimal transport for training a reconstruction
scheme for ill-posed inverse problems.
2.2

Optimal mass transport and Sinkhorn iterations

In optimal mass transport the aim is to transform one distribution into another by moving the mass in
a way that minimizes the cost of the movement. For an introduction and overview of the topic, see,
e.g., the monograph [35]. Lately, the area has attracted a lot of research [10, 11, 9] with applications
to, e.g., signal processing [17, 15, 19, 13] and inverse problems [7, 22].
The optimal mass transport problem can be formulated as follows: let Ω ⊂ Rd be a compact set, and
let µ0 and µ1 be two measures, defined on Ω, with the same total mass. Given a cost c : Ω × Ω → R+
that describes the cost for transporting a unit mass from one point to another, find a (mass preserving)
transference plan M that is as cheap as possible. Here, the transference plan characterizes how to
move the mass of µ0 in order to deform it into µ1 . Letting the transference plan be a nonnegative
measure dM on the space Ω × Ω yields a linear programming problem in the space of measures:
Z
T (µ0 , µ1 ) =
min
c(x0 , x1 )dM (x0 , x1 )
(4)
dM ≥0
(x0 ,x1 )∈Ω×Ω
Z
subject to µ0 (x0 )dx0 =
dM (x0 , x1 ),
x ∈Ω
Z 1
µ1 (x1 )dx1 =
dM (x0 , x1 ).
x0 ∈Ω

Although this formulation is only defined for measures µ0 and µ1 with the same total mass, it can also
be extended to handle measures with unbalanced masses [15, 9]. Moreover, under suitable conditions
one can define the Wasserstein metrics Wp using T , by taking c(x0 , x1 ) = d(x0 , x1 )p for p ≥ 1
and where d is a metric on Ω, and Wp (µ0 , µ1 ) := T (µ0 , µ1 )1/p [35, Definition 6.1]. As the name
indicates, Wp is a metric on the set of nonnegative measures on Ω with fixed mass [35, Theorem 6.9],
and T is weak∗ continuous on this set. One important property is that T (and thus also Wp ) does not
only compare objects point by point, as standard Lp metrics, but instead quantifies how the mass is
moved. This makes optimal transport natural for quantifying uncertainty and modelling deformations
[19, 21].
One way to solve the optimal transport problem in applications is to discretize Ω and solve the
corresponding finite-dimensional linear programming problem. In this setting the two measures are
represented by point masses on the discretization grid, i.e., by two vectors µ0 , µ1 ∈ Rn+ where the
element [µk ]i corresponds to the mass in the point x(i) ∈ Ω for i = 1, . . . , n and k = 0, 1. Moreover,
a transference plan is represented by a matrix M ∈ Rn×n
where the value mij := [M ]ij denotes
+
the amount of mass transported from point x(i) to x(j) . The associated cost of a transference plan
Pn
is i,j=1 cij mij = trace(C T M ), where [C]ij = cij = c(x(i) , x(j) ) is the transportation cost from
3

x(i) to x(j) , and by discretizing the constraints we get that M is a feasible transference plan from µ0
to µ1 if the row sums of M is µ0 and the column sums of M is µ1 . The discrete version of (4) thus
takes the form
T (µ0 , µ1 ) =

min

M ≥0

trace(C T M )

subject to µ0 = M 1n

(5)
µ1 = M T 1n ,

where M ≥ 0 denotes element-wise non-negativity of the matrix. However, even though (5) is a
linear programming problem it is in many cases computationally infeasible due to the vast number of
variables. Since M ∈ Rn×n
the number of variables is n2 , and thus if one seek to solve the optimal
+
transport problem between two 512 × 512 images this results in more than 6 · 1010 variables.
One approach for addressingP
this problem was proposed by Cuturi [10] that introduces an entropic
n
regularizing term D(M ) = i,j=1 (mij log(mij ) − mij + 1) for approximating the transference
plan, so the resulting perturbed optimal transport problem reads as
min

M ≥0

trace(C T M ) + εD(M )

subject to µ0 = M 1n

(6)

µ1 = M T 1n .

One can show that an optimal solution to (6) is of the form
M = diag(u)Kdiag(v),

(7)
Rn+

where K = exp(−C/ε) (point-wise exponential) is known, and u, v ∈
are unknown. This shows
that the solution is parameterized by only 2n variables. Moreover, the two vectors can be computed
iteratively by so called Sinkhorn iterations, i.e., alternatingly compute u and v that matches µ0 and
µ1 respectively. This is summarizied in Algorithm 1 where denotes elementwise multiplication
and ./ elementwise division. The procedure has been shown to have a linear convergence rate, see
[10] and references therein.
Moreover, when the underlying cost c(x0 , x1 ) is translation invariant the discretized cost matrix C,
and thus also the transformation K, gets a Toeplitz-block-Toeplitz structure. This structure can be
used in order to compute Kv and K T u efficiently using the fast Fourier transform in O(n log n),
instead of naive matrix-vector multiplication in O(n2 ) [22]. This is crucial for applications in imaging
since for images of size 512×512 pixels one would have to explicitly store and multiply with matrices
of size 262144 × 262144 .
Algorithm 1 Sinkhorn iterations for computing entropy-regularized optimal transport [10]
1:
2:
3:
4:
5:
6:

3

Input C, ε, µ0 , µ1
initialize v0 > 0 and K = exp(−C/ε)
for i = 1, . . . , N do
ui ← µ0 ./(Kvi−1 )
vi ← µ1 ./(K T ui )
return uTN (K C)vN

Learning a reconstruction operator using Wasserstein loss

In this work we propose to use entropy regularized optimal transport (6) to train a reconstruction
operator, i.e., to select the parameters as
h
i
Θ∗ ∈ arg min Ef,g T (T †Θ (g), f) .
(8)
Θ

This should give better results when data g is not aligned with the ground truth f . To see this,
consider the case when f is a point mass. In that case training the network with the L2 loss (3)
will (in the ideal case) result in a perfect reconstruction composed with a convolution that “smears”
the reconstruction over the area of possible miss-alignment. On the other hand since optimal mass
transport does not only compare objects point-wise, the network will (in the ideal case) learn a perfect
reconstruction combined with a movement of the object to the corresponding barycenter (centroid) of
the miss-alignment. These statements are made more precise in the following propositions. Formal
definitions and proofs are deferred to the appendix.
4

Proposition 1. Let g ∈ L2 (Rn ), let τ be a Rn -valued random variable with probability measure
dP (t), and let gτ (x) := g(x − τ ). Then there exists a function f ∈ L2 (Rn ) that minimizes
Eτ [kf − gτ k22 ], and this f has the form
Z
f (x) = (dP ∗ g)(x) :=
g(x − t)dP (t).
Rn

Proposition 2. Let δ(x) be the Dirac delta function on Rn , let τ be a Rn -valued random variable
with probability measure dP (t), and let δτ (x) := δ(x − τ ). Then
Z Z
Z

Eτ [T (δτ , µ)] =
c(t, x)dP (t) dµ(x) for µ ≥ 0 and
dµ(x) = 1
n
Rn
Rn
|R
{z
}
:=F (x)

and Eτ [T (δτ , µ)] = ∞ otherwise. Furthermore, finding a µ that minimizes Eτ [T (δτ , µ)] is equivalent
to finding the global minimizers to F (x). In particular, if (i) the probability measure dP is symmetric
around its mean, (ii) the underlying cost c is of the form c(t, x) = d(x − t), where d is convex and
symmetric, and (iii) dP and d are such that
Z
Z
d(x − t)dP (t)
and
∂d(x − t)dP (t)
Rn

Rn
n

both exist and are finite for all x ∈ R , then µ(x) = δ(x − E[τ ]) is an optimal solution. Furthermore,
if d is also strictly convex, then this is the unique minimizer.
To illustrate Propositions 1 and 2 we consider the following example.
Example 3. Let τ be uniformly distributed on [−1, 1], and let c(x0 , x1 ) = (x0 − x1 )2 . This gives
Z
1 1
1
F (x) =
(x − t)2 dt = + x2 ,
2 −1
3
which has minimum x = 0, and hence the (unique) minimizer to Eτ [T (δτ , µ)] is µ(x) = δ(x). For
the L2 case with the uniform distribution, the minimizer of Eτ [kf − gτ k22 ] is the smoothed function
g ∗ 12 χ[−1,1] .
The most common choice of distance c is to use the squared norm c(x0 , x1 ) = kx0 − x1 k2 , as in
the previous example. In this case the result of Proposition 2 can be strengthened, as shown in the
following example.
Example 4. Let τ be a Rn -valued random variable with probability measure dP (t) with finite first
and second moments, and let c(x0 , x1 ) = kx0 − x1 k2 . This gives
Z
F (x) =
(x − t)2 dP (t) = x2 − 2x E[τ ] + E[τ 2 ],
Rn

which has a unique global minimum in x = E[τ ] and hence µ(x) = δ(x − E[τ ]).

4

Implementation and evaluation

We use the recently proposed learned primal-dual structure in [3] for learning a reconstruction
operator T †Θ for solving the inverse problem in (1). In this algorithm, a sequence of small blocks work
alternatingly in the data (dual) space Y and the reconstruction (primal) space X and are connected
using the forward operator T and its adjoint T ∗ . The algorithm works with any differentiable operator
T , but we state the version for linear operators for simplicity in algorithm 2.
Algorithm 2 Learned Primal-Dual reconstruction algorithm
1: Initialize f0 ∈ X Nprimal , h0 ∈ U Ndual
2: for i = 1, . . . , I do

(2)
3:
hi ← ΓΘd hi−1 , T (fi−1 ), g
i

(1)

fi ← ΛΘpi fi−1 , T ∗ (hi )
(1)
†
5: T Θ (g) := fI

4:



5

7

32

32

5
6

32

32

5

apply
apply
copy
3x3 conv + PReLU
3x3 conv

Figure 1: Network architecture used to solve the inverse problem. Dual and primal iterates are in
blue and red boxes, respectively. Several arrows pointing to the same box indicates concatenation.
The initial values f0 , h0 enter from the left, while the data g is supplied to the dual iterates.

→
−
(a) Phantom

→
−
(b) Translated phantom

(c) Data

Figure 2: Example of data generation process used for training and validation, where 2a shows
an example phantom, 2b is the phantom with a random translation and 2c is the data (sinogram)
corresponding to 2b with additive white noise on top. The pair (gi , fi ) = (2c, 2a) is what is used in
the training.

The method was implemented using ODL [2], ASTRA [34], and TensorFlow [1]. We used the
reference implementation1 with default parameters, i.e., the number of blocks in the primal and dual
space was I = 10, and the number of primal and dual variables was set to Nprimal = Ndual = 5.
Moreover, the blocks used a residual structure and had three layers of 3×3 convolutions with 32 filters.
PReLU nonlinearities were used. Thus, this corresponds to a residual CNN with convolutional depth
of 10 · 2 · 3 = 60, as shown in graphical format in fig. 1. We used zero initial values, f0 = h0 = 0.
We compare a learned reconstruction operator of this form when trained using L2 loss (3) and using
optimal transport loss (8). Moreover, the evaluation is done on a problem similar to the evaluation
problem in [3, 4], i.e., on a problem in computed tomography. More specifically, training is done on
data that consists of randomly generated circles on a domain of 512 × 512 pixel, and the forward
operator T is the ray transform [27]. What makes this an ill-posed problem is that the data acquisition
is done from only 30 views with 727 parallel lines. Moreover, the source of noise is two-fold in
this set-up: (i) the pairs (gi , fi ) of data sets and phantoms are not aligned, meaning that the data is
computed from a phantom with a random change in position. This random change is independent
for the different circles, and for each circle it is a shift which is uniformly distributed over [−40, 40]
pixels, both in up-down and left-right direction. (ii) on the data computed from the shifted phantom,
5% additive Gaussian noise was added. For an example, see fig. 2.
The optimal mass transport distance computed with Sinkhorn iterations was used as loss function,
where we used the transport cost


4
4
c(x1 , x2 ) = 1 − e−kx1 −x2 k /80 .
1

https://github.com/adler-j/learned_primal_dual

6

(a) Phantom

(b) Translated phantom (c) Mean squared error loss (d) Optimal transport loss

Figure 3: In 3a we show the validation phantom, which was generated from the same training set
but not used in training, in 3b the translated phantom from which the validation data was computed,
in 3c a reconstruction with neural network trained using mean squared error loss (3), and in 3d a
reconstruction with neural network trained using optimal mass transport loss (8).

This was chosen since it heavily penalizes large movements, while not diverging to infinity which
causes numerical instabilities. Moreover, c(x1 , x2 )1/4 is in fact a metric on R2 (see lemma 6 in the
appendix) and thus W4 (µ0 , µ1 ) := T (µ0 , µ1 )1/4 gives rise to a Wasserstein metric on the space of
images, where T (µ0 , µ1 ) is the optimal mass transport distance with the transport cost c(x1 , x2 ).
Since this cost is translation invariant, the matrix-vector multiplications Ku and K T v can be done
with fast Fourier transform, as mentioned in section 2.2, and this was implemented in Tensorflow. We
used 10 Sinkhorn iterations with entropy regularization ε = 10−3 , to approximate the optimal mass
transport. Automatic differentiation was used to back-propagate the result during training.
Since the optimal mass transport function (6) is only finite for marginals µ0 and µ1 with the
same total mass, in the training we normalize the output of the reconstruction T †Θ (g) with
mass(f )/mass(T †Θ (g)). This makes T †Θ (g) invariant with respect to the total mass, which is
undesirable. To compensate for this, a small penalization on the error in total mass was added to the
loss function.
The training also followed [3] closely. In particular, we used 2 · 104 batches of size 1, using the ADAM
optimizer [23] with default values except for β2 = 0.99. The learning rate (step length) used was
cosine annealing [25] with initial step length 10−3 . Moreover, in order to improve training stability
we performed gradient norm clipping [29] with norms limited to 1. The convolution parameters
were initialized using Xavier initialization [16], and all biases were initialized to zero. The training
took approximately 3 hours using a single Titan X GPU. The source code used to replicate these
experiments are available online 2
Results are presented in fig. 3. As can be seen, the reconstruction using L2 loss “smears” the reconstruction to an extent where the shape is impossible to recover. On the other hand, the reconstruction
using the Wasserstein loss retains the over-all global shape of the object, although relative and exact
positions of the circles are not recovered.

5

Conclusions and future work

In this work we have considered using Wasserstein loss to train a neural network for solving ill-posed
inverse problems in imaging where data is not aligned with the ground truth. We give a theoretical
motivation for why this should give better results compared to standard mean squared error loss, and
demonstrate it on a problem in computed tomography. In the future, we hope that this method can be
applied to other inverse problems and to other problems in imaging such as segmentation.

2

https://github.com/adler-j/wasserstein_inverse_problems

7

Appendix: Deferred definition and proofs
Proof of Proposition 1. To show that f (x) = (dP ∗ g)(x) ∈ L2 (Rn ) minimizes Eτ [kf − gτ k22 ] we
expand the expression and use Fubini’s theorem to get


Z Z
Z Z
2
2
2
Eτ [kf − gτ k2 ] =
|f (x) − gt (x)| dx dP (t) =
(f (x) − gt (x)) dP (t) dx.
Rn

Rn

Rn

Rn

R

dP (t)dt = 1, this can be written as
2
Z 
Z
Eτ [kf − gτ k22 ] =
f (x) −
gt (x)dP (t) dx + c,

Rearranging terms and using that

Rn

Rn

Rn

where c is a constant. Using this it follows that the minimizing f is of the form
Z
f (x) =
gt (x)dP (t).
Rn

To see that f ∈ L2 (Rn ) we note that, by using Fubini’s theorem, we have
2

Z Z
Z Z Z
2
kf k2 =
gt (x)dP (t) dx =
gs (x)gt (x)dx dP (s)dP (t)
Rn
Rn
Rn Rn
Rn

Z Z  Z
1
≤
gs (x)2 + gt (x)2 dx dP (s)dP (t) = kgk22 < ∞
2
Rn Rn
Rn
where the first inequality is the arithmetic-geometric mean inequality. This completes the proof.
Definition 5. Let h : Rn → R. A subgradient to h in a point y is a vector v so that
∀ x ∈ Rn .

h(x) ≥ h(y) + hv, x − yi,

The set of all subgradients in a point y is called the subdifferential of h at y, and is denoted by
n
R∂h(y). This is a set-valued operator, and for any measure dν on R we define (dν ∗ ∂h)(y) :=
∂h(t)dν(y
−
t)
to
be
the
set-valued
operator
Rn
Z

y 7→
v(t)dν(y − t) ∈ Rn | v(t) ∈ ∂h(t) .
Rn

Proof of Proposition 2. We consider finding the marginal µ that minimize Eτ [T (δτ , µ)]. Without
loss of generality we assume that τ is zero-mean, since otherwise we simply consider τ − E[τ ] which
is a zero-mean random variable. First we note that T (δt , µ) is only finite for nonegative measures µ
with total mass 1, and hence Eτ [T (δτ , µ)] is only finite for such measures. Second, for such a µ we
have
Z
T (δt , µ) =
c(t, x)dµ(x),
Rn

since one needs to transport all mass in µ into the point t where δt has its mass. Using this and
expanding the expression for the expectation gives that

Z
Z Z
Eτ [T (δτ , µ)] =
T (δt , µ)dP (t) =
c(t, x)dP (t) dµ(x),
Rn

Rn

Rn

where we have used Fubini’s theorem in the last step. This completes the first half of the statement.
To prove the second half of the statement, note that the optimal µ have support only in the global
minimas of the function
Z
Z
F (x) :=
d(x − t)dP (t) =
d(t)dP (x − t),
Rn

Rn

which by assumption exists and is finite. Now, since d is convex we have that
d(x) ≥ d(y) + h∂d(y), x − yi,
8

for all x, y ∈ Rn ,

(9)

and convolving this inequality with dP gives the inequality
Z
F (x) ≥ F (y) + h
∂d(y − t)dP (t), x − yi,

for all x, y ∈ Rn ,

(10)

Rn

where all terms exist and are bounded by assumption. This shows that
Z
∂d(y − t)dP (t) ⊂ ∂F (y).
Rn

Now, since d is symmetric we have that ∂d is anti-symmetric, i.e., that ∂d(x) = −∂d(−x), since
d(x) = d(−x) ≥ d(−y) + h∂d(−y), −x + yi = d(y) + h−∂d(−y), x − yi.
Therefore
Z
∂F (0) ⊃
∂d(−t)dP (t) 3 0,
Rn

where the last inclusion follows since dP is symmetric and ∂d is anti-symmetric. Now, since
0 ∈ ∂F (0) we have that x = 0 is a global minimizer to F (x) [6, Theorem 16.2], and thus one
optimal solution to the problem is µ(x) = δ(x). Now, if d is strictly convex, the inequality (9) is
strict for x 6= y, and thus (10) is also strict, which shows that the optimal solution is unique.
Lemma 6. Let k · k be a norm on Rm . Then
d(x1 , x2 ) = 1 − e−kx1 −x2 k

n

 n1

.

is a metric on Rm for n ≥ 1.
Proof. It is easily seen that d(x1 , x2 ) is symmetric, nonnegative, and equal to zero if only if x1 = x2 .
Thus we only need to verify that the triangle inequality holds. To this end we note that if
n 1
n 1
n 1
1 − e−(a+b) n ≤ 1 − e−a n + 1 − e−b n , for all a, b ≥ 0,
(11)
for all n ≥ 1, then by taking a = kx1 − x2 k, b = kx2 − x3 k, and using the triangle inequality for
the norm k · k we have that
n 1
n 1
d(x1 , x3 ) = 1 − e−kx1 −x3 k n ≤ 1 − e−(kx1 −x2 k+kx2 −x3 k) n
n 1
n 1
≤ 1 − e−kx1 −x2 k n + 1 − e−kx2 −x3 k n = d(x1 , x2 ) + d(x2 , x3 ).
Therefore we will show that (11) holds for all n ≥ 1, and to do so we will
(i) show that if a function g : R+ → R+ fulfills g(0) = 0, g(x)0 ≥ 0, g 00 (x) ≤ 0 for all
x ∈ R+ , then g(x1 + x2 ) ≤ g(x1 ) + g(x2 ),
n

1

(ii) show that for x ≥ 0 the map x 7→ (1 − e−x ) n fulfills the assumptions in (i) for any n ≥ 1.
To show (i) we note that
x1 +x2

Z x1
Z x1 +x2
g 0 (t)dt =
g 0 (t)dt +
g 0 (t)dt
0
0
x1
Z x1
Z x2
0
0
≤
g (t)dt +
g (t)dt = g(x1 ) + g(x2 ),
Z

g(x1 + x2 ) =

0

0

where the inequality uses that g 0 (t) ≥ g 0 (x + t) for any x, t ≥ 0 since g 00 (x) ≤ 0 for all x ≥ 0.
n

1

To show (ii), let g(x) := (1 − e−x ) n and observe that g(0) = 0. Differentiating g twice gives
n

g 0 (x) =

n

1

e−x (1 − e−x ) n xn−1
1 − e−xn
=:hn (xn )

z
}|
{
1
xn n
n
xn
xn
−xn n
n−2
(1
−
e
)
x
ne
x
−
x
+
e
−
ne
+
n
−
1
g 00 (x) = −
.
(exn − 1)2
For x ≥ 0 we see that g 0 (x) ≥ 0 for all n ≥ 1. Moreover, for x ≥ 0 we see that g 00 (x) ≤ 0 for all
x ≥ 0 and for all n ≥ 1 if and only if hn (xn ) ≥ 0. With the change of variable xn = y, we thus
want to show that hn (y) ≥ 0 for all y ≥ 0 and all n ≥ 1. To see this we note that hn (0) = 0 and that
h0n (y) = ney y + ey − 1 ≥ 0 for all y ≥ 0 and n ≥ 1.
This shows (ii), and hence completes the proof.
9

References
[1] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. Corrado, A. Davis, J. Dean, M. Devin,
et al. Tensorflow: Large-scale machine learning on heterogeneous distributed systems. arXiv preprint
arXiv:1603.04467, 2016.
[2] J. Adler, H. Kohr, and O. Öktem. ODL - a Python framework for rapid prototyping in inverse problems.
2017. In preparation, KTH Royal Institute of Technology. Code and documentation available online:
https://github.com/odlgroup/odl.
[3] J. Adler and O. Öktem. Learned primal-dual reconstruction. arXiv preprint arXiv:1707.06474, 2017.
[4] J. Adler and O. Öktem. Solving ill-posed inverse problems using iterative deep neural networks. Inverse
Problems, 2017.
[5] M. Arjovsky, S. Chintala, and L. Bottou. Wasserstein GAN. arXiv preprint arXiv:1701.07875, 2017.
[6] H.H. Bauschke and P.L. Combettes. Convex analysis and monotone operator theory in Hilbert spaces.
Springer, New York, 2011.
[7] J.-D. Benamou, G. Carlier, M. Cuturi, L. Nenna, and G. Peyré. Iterative Bregman projections for regularized
transportation problems. SIAM Journal on Scientific Computing, 37(2):A1111–A1138, 2015.
[8] R.W. Brown, Y.-C. N. Cheng, E.M. Haacke, M.R. Thompson, and R. Venkatesan. Magnetic Resonance
Imaging: Physical Principles and Sequence Design. John Wiley & Sons Ltd, 2014.
[9] L. Chizat, G. Peyré, B. Schmitzer, and F.-X. Vialard. Unbalanced optimal transport: geometry and
Kantorovich formulation. arXiv preprint arXiv:1508.05216, 2015.
[10] M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural
Information Processing Systems, pages 2292–2300, 2013.
[11] M. Cuturi and G. Peyré. A smoothed dual approach for variational Wasserstein problems. SIAM Journal
on Imaging Sciences, pages 320–343, 2016.
[12] H.W. Engl, M. Hanke, and A. Neubauer. Regularization of inverse problems. Kluwer Academic Publisher,
2000.
[13] B. Engquist and B.D. Froese. Application of the Wasserstein metric to seismic signals. Communications
in Mathematical Sciences, 12(5), 2014.
[14] C. Frogner, C. Zhang, H. Mobahi, M. Araya, and T.A. Poggio. Learning with a Wasserstein loss. In
Advances in Neural Information Processing Systems, pages 2053–2061, 2015.
[15] T.T. Georgiou, J. Karlsson, and M.S. Takyar. Metrics for power spectra: an axiomatic approach. IEEE
Transactions on Signal Processing, 57(3):859–867, 2009.
[16] X. Glorot and Y. Bengio. Understanding the difficulty of training deep feedforward neural networks. In
Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages
249–256, 2010.
[17] S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent. Optimal mass transport for registration and warping.
International Journal of computer vision, 60(3):225–240, 2004.
[18] V. Jain and S. Seung. Natural image denoising with convolutional networks. In D. Koller, D. Schuurmans,
Y. Bengio, and L. Bottou, editors, Advances in Neural Information Processing Systems 21, pages 769–776.
Curran Associates, Inc., 2009.
[19] X. Jiang, Z.-Q. Luo, and T.T. Georgiou. Geometric methods for spectral analysis. IEEE Transactions on
Signal Processing, 60(3):1064–1074, 2012.
[20] J. Johnson, A. Alahi, and L. Fei-Fei. Perceptual Losses for Real-Time Style Transfer and Super-Resolution,
pages 694–711. Springer International Publishing, Cham, 2016.
[21] J. Karlsson and T.T. Georgiou. Uncertainty bounds for spectral estimation. IEEE Transactions on Automatic
Control, 58(7):1659–1673, 2013.
[22] J. Karlsson and A. Ringh. Generalized Sinkhorn iterations for regularizing inverse problems using optimal
mass transport. arXiv preprint arXiv:1612.02273, 2016.
[23] D. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
[24] Y. LeCun, Y. Bengio, and G. Hinton. Deep learning. Nature, 521:436—444, 2015.
[25] I. Loshchilov and F. Hutter.
arXiv:1608.03983, 2016.

SGDR: stochastic gradient descent with restarts.

arXiv preprint

[26] M. Mardani, E. Gong, J.Y. Cheng, S. Vasanawala, G. Zaharchuk, M.T. Alley, N. Thakur, S. Han, W.J.
Dally, J.M. Pauly, and L. Xing. Deep generative adversarial networks for compressed sensing automates
MRI. CoRR, abs/1706.00051, 2017.

10

[27] F. Natterer and F. Wübbeling. Mathematical methods in image reconstruction. SIAM, 2001.
[28] O. Öktem. Mathematics of electron tomography. In O. Scherzer, editor, Handbook of Mathematical
Methods in Imaging, pages 937–1031. Springer, New York, NY, 2015.
[29] R. Pascanu, T. Mikolov, and Y. Bengio. Understanding the exploding gradient problem. CoRR,
abs/1211.5063, 2012.
[30] P. Paschalis, N. D. Giokaris, A. Karabarbounis, G. K. Loudos, D. Maintas, C. N. Papanicolas,
V. Spanoudaki, Ch. Tsoumpas, and E. Stiliaris. Tomographic image reconstruction using artificial neural
networks. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers,
Detectors and Associated Equipment, 527(1–2):211–215, 2004.
[31] D. M. Pelt and K. J. Batenburg. Fast tomographic reconstruction from limited data using artificial neural
networks. IEEE Transactions on Image Processing, 22(12):5238–5251, 2013.
[32] P. Putzky and M. Welling. Recurrent inference machines for solving inverse problems. Submitted to ICLR
2017, Toulon, France, April 24–26, 2017. Report available from https://openreview.net/pdf?id=HkSOlP9lg,
2017.
[33] L.I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. Physica D:
Nonlinear Phenomena, 60(1):259–268, 1992.
[34] W. van Aarle, W.J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. De Beenhouwer, K.J.
Batenburg, and J. Sijbers. Fast and flexible X-ray tomography using the ASTRA toolbox. Optics express,
24(22):25129–25147, 2016.
[35] C. Villani. Optimal transport: old and new, volume 338. Springer Science & Business Media, 2008.
[36] A. Walker, G. Liney, P. Metcalfe, and L. Holloway. MRI distortion: considerations for MRI based
radiotherapy treatment planning. Australasian Physical & Engineering Sciences in Medicine, 37(1):103–
113, Mar 2014.
[37] T. Würfl, F. C. Ghesu, V. Christlein, and A. Maier. Deep learning computed tomography. In S. Ourselin,
L. Joskowicz, M. Sabuncu, G. Unal, and W. Wells, editors, MICCAI 2016: Medical Image Computing
and Computer-Assisted Intervention – MICCAI 2016, volume 9902 of Lecture Notes in Computer Science,
pages 432–440. Springer-Verlag, 2016.
[38] J. Xie, L. Xu, and E. Chen. Image denoising and inpainting with deep neural networks. In F. Pereira,
C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, Advances in Neural Information Processing
Systems 25, pages 341–349. Curran Associates, Inc., 2012.
[39] Y. Yang, J. Sun, H. Li, and Z. Xu. Deep ADMM-Net for compressive sensing MRI. In D. D.
Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 29, pages 10–18. Curran Associates, 2016. Report available from
http://papers.nips.cc/paper/6406-deep-admm-net-for-compressive-sensing-mri.pdf.
[40] H. Zhao, O. Gallo, I. Frosio, and J. Kautz. Loss functions for image restoration with neural networks.
IEEE Transactions on Computational Imaging, 3(1):47–57, March 2017.

11

