Year: 2019

Abstract:
The penalized least squares (PLS) is a classic method to inverse problems, where a regularization term is added to stabilize the solution. Optimal transport (OT) is another mathematical framework for computer vision tasks that provides means to transport one measure to another at minimal cost. The cycle-consistent generative adversarial network (cycleGAN) is a recent extension of GAN to learn target distributions with less mode collapsing behavior. Although similar in that no supervised training is required, the algorithms look different, so the mathematical relationship between these approaches is not clear. In this article, we provide an important advance to unveil the missing link. Specifically, we reveal that a cycleGAN architecture can be derived as a dual formulation of the OT problem, if the consistency constraint of PLS is enforced as a regularization term of the OT problem. This suggests that cycleGAN can be considered stochastic generalization of classical PLS approaches. Our derivation is so general that various types of cycleGAN architectures can be easily derived by merely changing the transport cost. As proofs of concept, this paper provides a novel cycleGAN architecture for unsupervised learning in accelerated magnetic resonance imaging (MRI) and deconvolution microscopy problems, which confirm the efficacy and the flexibility of the theory.

heading: 1 INTRODUCTION
text: Inverse problems are ubiquitous in computer vision, biomedical imaging, scientific discovery, etc. In inverse problems, a noisy measurement y ∈ Y from an unobserved image x ∈ X is modeled by
y = Hx+ w , (1)
where w is the measurement noise, and H : X 7→ Y is the measurement operator. In inverse problems originating from physics, the measurement operator is usually represented by an integral equation:
Hx(r) := ∫ Rd h(r, r′)x(r′)dr′, r ∈ D ⊂ Rd, d = 2, 3, (2)
where h(r, r′;x) is an integral kernel. Then, the inverse problem is formulated as an estimation problem of the unknown x from the measurement y. It is well known that inverse problems are ill-posed. A classical strategy to mitigate the ill-posedness is the penalized least squares (PLS) approach:
x̂ = arg min x c(x; y) := ‖y −Hx‖q +R(x) (3)
for q ≥ 1, whereR(x) is a regularization (or penalty) function (l1, total variation (TV), etc.) (Chaudhuri et al., 2014; Sarder & Nehorai, 2006; McNally et al., 1999). In some inverse problems, the measurement operator H is not well defined, so both the unknown operator H and the image x should be estimated.
Recently, deep learning approaches with supervised training have become the mainstream approaches for inverse problems because of their excellent and ultra-fast reconstruction performance.
For example, in low-dose x-ray computed tomography (CT) denoising problems, a convolutional neural network (CNN) is trained to learn the relationship between the noisy image y and the matched noiseless (or high-dose) label images x (Kang et al., 2017). In the context of (3), the supervised neural network can be understood as directly learning the operation x̂ = arg minx c(x; y). Unfortunately, in many applications, matched label data are not available. Therefore, unsupervised training without matched reference data has become an important research topic.
Recently, the generative adversarial network (GAN) has attracted significant attention in the machine learning community by providing a way to generate target data distribution from random distribution (Goodfellow et al., 2014). In particular, Arjovsky et al. (2017) proposed the so-called Wasserstein GAN (W-GAN), which is closely related to the mathematical theory of optimal transport (OT) (Villani, 2008; Peyré et al., 2019). In OT, for two given probability measures supported on the X and Y spaces, one pays a cost for transporting one measure to another. Then, the minimization of the average transportation cost provides an unsupervised way of learning the transport map between the two measures. Unfortunately, these GAN approaches often generate artificial features due to mode collapsing, so cycle-consistent GAN (cycleGAN) (Zhu et al., 2017), which imposes one-toone correspondence, has been extensively investigated (Kang et al., 2019; Lu et al., 2017). Recently, cycleGAN was generalized to bridge two different domains by generating a continuous sequence of intermediate domains flowing from one domain to the other (Gong et al., 2019).
Although classical PLS, OT, and cycleGAN share the commonality of unsupervised learning which does not require matched training data, there is no mathematical theory to systematically link these seemingly different approaches. Therefore, one of the main contributions of this paper is to unveil the missing link between these methods. In particular, we reveal that a cycleGAN architecture can be derived as a dual formulation of the OT problem when the consistency term in PLS is enforced as a regularization term for the OT problem. Moreover, this framework is so general that various cycleGAN architectures can be easily obtained by changing the PLS cost. Accordingly, our framework provides a principled way of deriving unsupervised neural networks for various inverse problems.
As proofs of concept, we provide novel cycleGAN architectures for unsupervised learning in two physical inverse problems: accelerated MRI and deconvolution microscopy. In these applications, we show that only one CNN generator is required, which significantly reduces the training complexity. The experimental results confirm that the proposed unsupervised learning approaches can successfully provide accurate inversion results without any matched reference. All the proofs of the main theoretical results can be found in the Appendix.

heading: 2 RELATED WORKS

heading: 2.1 OPTIMAL TRANSPORT (OT)
text: OT compares two measures in a Lagrangian framework (Villani, 2008; Peyré et al., 2019). Formally, we say that T : X 7→ Y transports µ ∈ P (X ) to ν ∈ P (Y), if
ν(B) = µ ( T−1(B) ) , for all ν-measurable sets B, (4)
which is often simply represented by ν = T#µ, where T# is often called the push-forward operator. Monge’s original OT problem (Villani, 2008; Peyré et al., 2019) is then to find a transport map T that transports µ to ν at the minimum total transportation cost:
min T
M(T ) := ∫ X c(x, T (x))dµ(x), subject to ν = T#µ
However, this is usually computationally expensive due to the nature of combinatorial assignment. Kantorovich relaxed the assumption to consider probabilistic transport that allows mass splitting from a source toward several targets (Villani, 2008; Peyré et al., 2019). Specifically, Kantorovich introduced a joint measure π ∈ P (X × Y) and the associated cost c(x, y), x ∈ X , y ∈ Y such that the original problem can be relaxed as
min π
K(π) := ∫ X×Y c(x, y)dπ(x, y) (5)
subject to π(A× Y) = µ(A), π(X ×B) = ν(B) for all measurable sets A ∈ X and B ∈ Y . Here, the last two constraints come from the observation that the total amount of mass removed from any measurable set has to equal the marginals (Villani,
2008; Peyré et al., 2019). Another important advantage of Kantorovich formulation is the dual formulation as stated in the following theorem: Theorem 1 (Kantorovich duality theorem). (Villani, 2008, Theorem 5.10, p.57-p.59) Let (X , µ) and (Y, ν) be two Polish probability spaces (separable complete metric space) and let c : X × Y → R be a continuous cost function, such that |c(x, y)| ≤ cX (x) + cY(y) for some cX ∈ L1(µ) and cY ∈ L1(ν), where L1(µ) denotes the set of 1-Lipschitz functions with the measure µ. Then, there is duality:
min πΠ(µ,ν) ∫ X×Y c(x, y)dπ(x, y) = max ϕ∈L1(µ) {∫ X ϕ(x)dµ(x) + ∫ Y ϕc(y)dν(y) } and the above maximum is taken over the so-called Kantorovich potential ϕ whose c-transform ϕc(y) = supx(c(x, y)− ϕ(x)) is properly defined.

heading: 2.2 PLS WITH DEEP LEARNING PRIOR
text: Recently, PLS frameworks using a deep learning prior have been extensively studied (Zhang et al., 2017; Aggarwal et al., 2018) thanks to their similarities to the classical regularization theory. The main idea of these approaches is to utilize a pre-trained neural network to stabilize the inverse solution. For example, in model based deep learning architecture (MoDL) (Aggarwal et al., 2018), the problem can be formulated as
min x c(x; y,Θ,H) = ‖y −Hx‖2 + λ‖x−QΘ(x)‖2 (6)
for some regularization parameter λ > 0, where QΘ(x) is a pre-trained CNN with the network parameter Θ and the input x. This problem is usually solved in an alternating minimization framework:
xn+1 = arg min x ‖y −Hx‖2 + λ‖x− zn‖2 , zn+1 = QΘ(xn+1) (7)
Another type of inversion approach using a deep learning prior is the so-called deep image prior (DIP) (Ulyanov et al., 2018). Rather than using an explicit prior, the deep neural network architecture itself is used as a regularization by restricting the solution space:
min Θ
c(Θ; y,H) = ‖y −HQΘ(z)‖2 (8)
where z is a random vector. Then, the final solution becomes x = QΘ∗(z) with Θ∗ being the estimated network parameters. Similarly, generative approaches (Van Veen et al., 2018; Bora et al., 2017; Wu et al., 2019) either estimate the random variable z from (8) by fixing Θ or attempt to estimate both the random z and the network weight Θ.

heading: 3 MAIN CONTRIBUTIONS
text: One of the basic assumptions in the aforementioned PLS formulation with a deep learning prior is that the measurement y is fixed and one is interested in finding the unknown x. In this paper, we relax the assumption to consider situations where both of them are random samples from joint probability measure π(x, y). Specifically, we propose a new PLS cost function with a novel deep learning prior as follows:
c(x, y; Θ,H) = ‖y −Hx‖q + λ‖GΘ(y)− x‖p (9)
with p, q ≥ 1, where the first two terms before the semi-colon are the random vectors. Here, the regularization parameter λ takes care of the dimensional difference of x and y, with λ = 1 for the same dimension, but for simplicity, in the rest of the paper, we assume λ = 1.
The new PLS cost in (9) has very unique properties. Suppose that the global minimum with c(x, y; Θ,H) = 0 exists. The corresponding global minimizer should then satisfy
y = Hx, x = GΘ(y) (10)
If this holds for all y, then GΘ is an inverse operator ofH. In practice, due to the limited capacity of the neural network, the global minimum with c(x, y; Θ,H) = 0 may never be achieved. Even in this
case, thanks to the symmetric form of PLS, the cost in (9) has an important implication: the neural network GΘ is now estimated by enforcing the consistency y = Hx as a regularization term. This implies that the roles of the consistency and regularization terms in the PLS can be interchanged in our formulation. Another important advantage is that the unknown image x can be obtained as an output of the feedforward neural network GΘ(y) for the given measurement y. This makes the inversion procedure much simpler. One may say that the PLS cost function for the DIP in (8) can produce a different variation of unsupervised learning. However, the resulting formulation is not a feed-forward neural network. In fact, to obtain the unknown image x, an additional optimization step is required even after the neural network is trained. See Appendix C for more details.
Note that Θ,H are the network and measurement system parameters that should be estimated. These parameters can be found by minimizing the average transport cost for all combinations of x ∈ X and y ∈ Y with respect to the joint measure π(x, y):
T(Θ,H) := min π ∫ X×Y c(x, y; Θ,H)dπ(x, y) (11)
where the minimum is taken over all joint distributions whose marginal distributions with respect to X and Y are µ and ν, respectively. Below, we show that the average transportation cost in (11) has an interesting decomposition, which leads to the cycleGAN architecture. Lemma 1. If the mapping GΘ : Y 7→ X is single-valued, then the average transportation cost T(Θ,H) in (11) can be decomposed as
T(Θ,H) = `cycle(Θ,H) + `OT ′(Θ,H) (12) where
`cycle(Θ,H) = min π ∫ Γ(GΘ)∪Γ(H) c(x, y; Θ,H)dπ(x, y) (13)
`OT ′(Θ,H) = min π ∫ X×Y\Γ(GΘ)∪Γ(H) c(x, y; Θ,H)dπ(x, y) (14)
and Γ(GΘ), Γ(H) denote the graphs of GΘ,H, respectively: Γ(GΘ) := {(x, y) ∈ X × Y | x = GΘ(y), y ∈ Y}, Γ(H) := {(x, y) ∈ X × Y | y = Hx, x ∈ X}
Note the nature of the sets Γ(GΘ) and Γ(H) and `cycle. Since both Γ(GΘ) and Γ(H) are graphs of functions, they are low-dimensional manifolds embedded in the ambient space X ×Y . In this sense, their relative measures become zero and `cycle seems to be a negligible part of the whole T. This term, however, can contribute a significant portion of the total loss T depending on the existence of the singularity in the distribution. We defer the discussion to after the following proposition, which explains why we name the first term `cycle. Proposition 1.
`cycle(Θ,H) := min π ∫ Γ(GΘ)∪Γ(H) c(x, y; Θ,H)dπ(x, y)
= ∫ Γ(GΘ)∪Γ(H) c(x, y; Θ,H)π∗(dx, dy)
= ∫ X ‖x−GΘ(Hx)‖pρ(x)µ(dx) + ∫ Y ‖y −HGΘ(y)‖qσ(y)ν(dy) (15)
where π∗ denotes the optimum measure and ρ and σ are measurable functions such that 0 ≤ ρ(x) ≤ 1 µ-a.e., and 0 ≤ σ(y) ≤ 1 ν-a.e.
The remaining term in (12), which corresponds to the GAN term, can be obtained as follows: Proposition 2. For p = q = 1, the cost `OT ′(Θ,H) in (14) can be computed as
`OT ′(Θ,H) = max ϕ ∫ X ϕ(x)[1− ρ(x)]dµ(x)− ∫ Y ϕ(GΘ(y))[1− σ(y)]dν(y) (16)
+ max ψ ∫ Y ψ(y)[1− σ(y)]dν(y)− ∫ X ψ(Hx)[1− ρ(x)]dµ(x) (17)
where ϕ and ψ are 1-Lipschitz Kantorovich potential functions.
The proof of Propositions 1 and 2 can be found in Appendix B. The proof idea comes from the computation of c(x, y; Θ,H) on subsets Γ(GΘ),Γ(H), and Γ(GΘ) ∩ Γ(H) by considering the singularity of the optimum measure π∗ on these sets. In this regard, it is important to know how singularity of π∗ is distributed on these manifolds. If the targets are perfectly matched as shown in Fig. 1(a), σ(y) := π∗({x = GΘ(y)}|y) = 1; ν-a.e so that π∗(·|y) is distributed like a point mass of magnitude σ(y) on {x = GΘ(y)}. This also leads to ρ(x) := π∗({y = Hx}|x) = 1;µ-a.e. In this case, `OT ′ = 0 and GΘ(y) is a perfect inverse to compute x. On the other hand, if the joint distribution π(x, y) is absolutely continuous as shown in Fig. 1(b), then ρ(x) = 0;µ-a.e. and σ(y) = 0; ν-a.e. In this case, `cycle = 0 and the contribution of the cycle consistency term can be neglected. The most interesting case arises when the joint measure π(x, y) has significant singularities on the manifold Γ(GΘ) ∪ Γ(H) so that 0 < σ(y) = π∗({x = GΘ(y)}|y) < 1 (see Fig. 1(c)). In fact, unsupervised learning applies to such situations where once the mappings are found, a significant portion of the data in X and Y can be paired with the mappings, even though unpaired data remain. This is when the cycle consistency terms plays an important role in unsupervised learning. Then, GΘ(y) is a desirable approximation of x with high probability.
Now, the final step is implementing the Kantorovich potential using CNNs with parameters Φ and Ξ, i.e. ϕ := ϕΦ and ψ := ψΞ. By collecting all terms in the transportation cost T(Θ,H) , the dual formulation results in a cycleGAN cost function:
min Θ,H T(Θ,H) = min Θ,h max Φ,Ξ `(Θ,H; Φ,Ξ) (18)
where
`(Θ,H; Φ,Ξ) = `cycle(Θ,H) + `GAN (Θ,H; Φ,Ξ) (19) where `cycle(Θ,H) denotes the cycle-consistency loss in (15) and `GAN (Θ, h; Φ,Ξ) is the GAN loss given by:
`GAN (Θ,H; Φ,Ξ) = ∫ X ϕΦ(x)[1− ρ(x)]dµ(x)− ∫ Y ϕΦ(GΘ(y))[1− σ(y)]dν(y) (20)
+ ∫ Y ψΞ(y)[1− σ(y)]dν(y)− ∫ X ψΞ(Hx)[1− ρ(x)]dµ(x)
Here, the Kantorovich 1-Lipschitz potentials ϕ := ϕΦ and ψ := ψΞ correspond to the W-GAN discriminators. Specifically, ϕΦ tries to find the difference between the true image x and the generated image GΘ(y), whereas ψ := ψΞ attempts to find the fake measurement data that are generated by the synthetic measurement procedure Hx. Furthermore, if ρ and σ are constant almost everywhere, they become simple scaling factors for the cycle-consistency terms:
`cycle(Θ,H) = ρ ∫ X ‖x−GΘ(Hx)‖pdµ(x) + σ ∫ Y ‖y −HGΘ(y)‖qdν(y) (21)
Note that our formulation using (19) with (20) and (21) has many unique features compared to the standard cycleGAN (Zhu et al., 2017).
• The standard cycleGAN (Zhu et al., 2017) usually requires two generators in the form of deep neural networks. On the other hand, our formulation often requires a single generator using a deep neural network, if the other generator is given by the forward model Hx. This significantly reduces the number of weights, which makes the training much more stable. This can be viewed as the consistency term from the forward model working as a regularization term for OT.
• When the solution of the dual problem achieves the global minimum with diminishing cost, then the dual solution x = GΘ(y) is equivalent to the primal solution (10), validating the dual approach. This is not the case in the conventional cycleGAN with two deep generators.
• While we only consider p = q = 1 due to the simple c-transform ϕc(x) = −ϕ(x) for 1-Lipschitz ϕ, the use of the general PLS cost would be interesting, and it may lead to an interesting variation of the cycleGAN architecture. This could be done using a regularized version of OT (Peyré et al., 2019).
• The parameters ρ(x) and σ(y) in (1) and (20) originate from the singularity of the optimal joint measure. However, the estimation of ρ(x) and σ(y) is not feasible in practice, so the constant hyperparameters in (21) are usually used, and the hyper parameters should be selected by trial and error, which is the same as in the standard cycleGAN. Nevertheless, the probability interpretation may lead to a more simplified search for the hyper-parameters.

heading: 4 APPLICATIONS TO INVERSE PROBLEMS

heading: 4.1 ACCELERATED MRI
text: In accelerated magnetic resonance imaging (MRI), the goal is to recover high-quality MR images from sparsely sampled k-space data to reduce the acquisition time. This problem has been extensively studied using compressed sensing (Lustig et al., 2007), but recently, deep learning approaches have been the main research interest due to their excellent performance and significantly reduced run-time complexity (Hammernik et al., 2018; Han & Ye, 2019). A standard method for neural network training for accelerated MRI is based on supervised learning, where the MR images from fully sampled k-space data are used as references and subsampled k-space data are used as the input for the neural network. Unfortunately, in accelerated MRI, high-resolution fully sampled k-space data are very difficult to acquire due to the long scan time. Therefore, the need for unsupervised learning without matched reference data is increasing.
In accelerated MRI, the forward measurement model can be described as
x̂ = PΩFx+ w (22)
whereF is the 2-D Fourier transform andPΩ is the projection to Ω that denotes k-space sampling indices. To implement every step of the algorithm as image domain processing, (22) can be converted to the image domain forward model by applying the inverse Fourier transform:
y = F−1PΩFx+ F−1w (23)
This results in the following cost function for the PLS formulation:
c(x, y; Θ) = ‖y −F−1PΩFx‖+ ‖GΘ(y)− x‖ (24)
where we set λ = 1 in (9) since the dimensions of x and y are the same. Usually, the sampling mask Ω is known so that the forward mapping for the inverse problem is deterministic.
The schematic diagram of the associated cycleGAN architecture is illustrated in Fig. 2, whose generator and discriminator architectures are shown in Fig. 6 in Appendix D. Note that we just need a
single generator and discriminator, as the mapping from the clean to aliased images is deterministic for a given sampling pattern. As for the loss function, we use (19) with the σ = ρ = 0.5 for the cycle consistency term in (21). For GAN implementation, we use the W-GAN in (20) with the Lipschitz penalty loss (Gulrajani et al., 2017) to ensure that the Kantorovich potential becomes 1-Lipschitz. The detailed description of the training procedure is provided in Appendix D. The reconstruction results in Fig. 3 and quantitative comparison results for all test sets in Table 1 clearly show that the proposed unsupervised learning method using dedicated cycleGAN successfully recovered fine details without matched references. Moreover, the performance is even better than that of the standard cycleGAN and is comparable with that of supervised learning.

heading: 4.2 DECONVOLUTION MICROSCOPY
text: Deconvolution microscopy is extensively used to improve the resolution of widefield fluorescent microscopy. However, classical deconvolution approaches require the measurement or estimation of the point spread function (PSF) and are usually computationally expensive. Recently, CNN approaches have been extensively studied as fast and high-performance alternatives. Unfortunately, the CNN approaches usually require matched high-resolution images for supervised training. Therefore, we are interested in developing unsupervised deep neural networks for deconvolution microscopy. Mathematically, a blurred measurement can be described as
y = h ∗ x+ w , (25)
where h is the PSF. Here, we consider the blind deconvolution problem where both the unknown PSF h and the image x should be estimated. This results in the following cost function for the PLS
formulation: c(x, y; Θ, h) = ‖y − h ∗ x‖+ ‖GΘ(y)− x‖. (26)
The corresponding cycleGAN architecture is illustrated in Fig. 4. In contrast to the conventional cycleGAN approaches that require two generators, the proposed cycleGAN approach needs only a single generator, and the blur image can be generated using a linear convolution layer corresponding to the PSF, which significantly improves the robustness of network training. However, we still need two discriminators, since both linear and deep generators should be estimated. The detailed descriptions of the network architecture and training are given in Appendix D.
Fig. 5 shows lateral views of the deconvolution results of microtubule samples by various methods. Here, input images are degraded by blur and noise. The supervised learning and the standard cycleGAN with two generators showed better contrast and removed blur; however, the structural continuity was not preserved. On the other hand, in our cycleGAN approach, blur and noise were successfully removed, and the continuity of the structure was preserved. The quantitative comparison with known ground-truth data is provided in Appendix D, which clearly confirms the superiority of the proposed method,

heading: 5 CONCLUSIONS
text: In this paper, we presented a general design principle for a cycleGAN architecture for various inverse problems. Specifically, the proposed architecture was derived as a dual formulation of an OT problem that enforces the data consistency of PLS as a regularization term for the OT problem. As proofs of concept, we designed novel cycleGAN architectures for accelerated MRI and deconvolution microscopy examples, providing accurate reconstruction results without any matched reference data. Given the generality of our design principle, we believe that our method can be an important platform for unsupervised learning for inverse problems.

heading: A ADDITIONAL THEOREMS AND LEMMA
text: Theorem 2 (Optimality is inherited by restriction). (Villani, 2008, Theorem 4.6, p.46) Under the same assumption of Kantorovich duality theorem, let π∗ be an optimal transport plan, and π̃ be a nonnegative measure on X × Y , such that π̃ ≤ π∗ and π̃[X × Y] > 0. Let π′ := π̃π̃[X×Y] , µ
′, ν′ be the maginals of π′. Then optimality is inherited to its restriction, i.e., π′ is an optimal transference plan between µ′ and ν′. Theorem 3 (Restriction for the Kantorovich duality). (Villani, 2008, Theorem 5.19, p.75-p.76) Under the same assumption of Kantorovich duality theorem, let π∗ be an optimal transport plan and ϕ∗ be an optimal potential that attain the maximum, and π̃ be a measure on X × Y , such that 0 ≤ π̃ ≤ π∗ and Z := π̃[X × Y] > 0. Let π′ = π̃/Z, and µ′, ν′ be the marginals of π′. Then there exists ϕ′ such that ϕ′ = ϕ∗ on projX (supp(π∗)), and ϕ′ solves the dual Kantorovich problem between (X , µ′) and (Y, ν′).
We are now ready to provide our duality result for the restricted measure. Lemma 2. Let
`OT ′ := min π ∫ X×Y\(A∪B) c(x, y)dπ(x, y)
`OT := min π ∫ X×Y c(x, y)dπ(x, y)
Then, we can replace `OT ′ with `OT in the sense that both of them share the same optimal Kantorovich potential. More specifically, if ϕ∗ is the optimal Kantorovich potential, which satisfies
`OT = ∫ X ϕ∗(x)dµ(x) + ∫ Y (ϕ∗)c(y)dν(y),
then
`OT ′ = ∫ X ϕ∗(x)[1− ρ(x)]dµ(x) + ∫ Y (ϕ∗)c(y)[1− σ(y)]dν(y).
Proof. This is a corollary of Theorem 2 and 3 by taking π̃ = π∗|X×Y\(Γ(GΘ)∪Γ(H)). Then, π̃ is nonnegative measure on X × Y . Use the same notation for Z, µ′, ν′ in Theorem 2 and 3. When Z = 0 , `OT ′ = 0, and µ′ = 0, ν′ = 0, so the statement trivially holds. Suppose Z > 0. Then by Theorem 3, there exists ϕ′ such that ϕ′ = ϕ∗ on projX (supp(π∗)), and ϕ′ solves the dual Kantorovich problem between (X , µ′) and (Y, ν′).
`OT ′/Z= ∫ X×Y c(x, y; Θ,H)dπ(x, y) Z
= ∫ X ϕ∗(x)dµ′(x) + ∫ Y (ϕ∗)c(y)dν′(y)
The remaining part of the proof is to show the relationship that µ′ = 1−ρZ µ and ν ′ = 1−σZ ν. From the proof of Proposition 1, the marginals of π∗|Γ(GΘ)∪Γ(H) are ρµ and σν. Therefore, the marginals of π̃ = π∗ − π∗|Γ(GΘ)∪Γ(H) are µ− ρµ and ν − σν so that Zµ′ = (1− ρ)µ and Zν′ = (1− σ)ν.
∴ `OT ′ = ∫ X ϕ∗(x)[1− ρ(x)]dµ(x) + ∫ Y (ϕ∗)c(y)[1− σ(y)]dν(y)

heading: B PROOF OF MAIN RESULTS

heading: B.1 PROOF OF LEMMA 1
text: Proof. From Theorem 2, we extract information about two separated integrals. For some optimal transportation plan π∗, the restrictions π∗|Γ(GΘ)∪Γ(H) and π∗|X×Y\(Γ(GΘ)∪Γ(H)) acquire optimality. Therefore,
T(Θ,H) = ∫ X×Y c(x, y; Θ,H)dπ∗(x, y)
= ∫ A∪B c(x, y; Θ,H)dπ∗(x, y) + ∫ X×Y\(A∪B) c(x, y; Θ,H)dπ∗(x, y)
= `cycle(Θ,H) + `OT ′(Θ,H)
This concludes the proof.

heading: B.2 PROOF OF PROPOSITION 1
text: Proof. By virtue of the definitions of the cost function and the sets Γ(GΘ) and Γ(H), we have
c(x, y; Θ,H) = ‖y −HGΘ(y)‖q on Γ(GΘ),
c(x, y; Θ,H) = ‖x−GΘ(Hx)‖p on Γ(H), c(x, y; Θ,H) = 0 on Γ(GΘ) ∩ Γ(H).
Hence,∫ Γ(GΘ)∪Γ(H) c(x, y; Θ,H)dπ∗(x, y) = ∫ Γ(H) ‖x−GΘ(Hx)‖pdπ∗(x, y)+ ∫ Γ(GΘ) ‖y−HGΘ(y)‖qdπ∗(x, y)
Furthermore, we use disintegration theorem (Simmons, 2012) to split the joint measure π∗ as follows: dπ∗(x, y) = π∗(dy|x)µ(dx) = π∗(dx|y)ν(dy), where π∗(Y|x) = 1 µ-a.e. and π∗(X|y) = 1 ν-a.e. Then, we have∫
Γ(H) ‖x−GΘ(Hx)‖pdπ∗(x, y) = ∫ Γ(H) ‖x−GΘ(Hx)‖pπ∗(dy|x)µ(dx)
= ∫ X ∫ y=Hx ‖x−GΘ(Hx)‖pπ∗(dy|x)µ(dx)
= ∫ X ‖x−GΘ(Hx)‖p ∫ 1y=Hxπ ∗(dy|x)µ(dx)
= ∫ X ‖x−GΘ(Hx)‖pρ(x)µ(dx)
where 1S denotes the indicator function for the set S and ρ(x) = ∫ 1y=Hxπ
∗(dy|x). In a similar fashion, with σ(y) = ∫ 1x=GΘ(y)π
∗(dx|y), we have∫ Γ(GΘ) c(x, y; Θ,H)dπ∗(x, y) = ∫ Y ‖y −HGΘ(y)‖qσ(y)ν(dy) .
This concludes the proof.

heading: B.3 PROOF OF PROPOSITION 2
text: Proof. This is a just simple corollary of the original Kantorovich’s duality formulation proof and the classical results of optimal transport on the restricted measure. We start with the first term of lOT (Θ,H) in (14). Specifically, we have∫ X×Y\(Γ(GΘ)∪Γ(H)) ‖GΘ(y)− x‖pπ∗(dx, dy) = max ϕ ∫ X ϕ(x)µ′(dx) + ∫ Y ϕc(GΘ(y))ν ′(dy)
where µ′ and ν′ are marginals of the restriction of optimal transportation plan on the restricted set X × Y \ (Γ(GΘ) ∪ Γ(H)). Now, using Lemma 2 in Appendix, we have∫
X×Y\(Γ(GΘ)∪Γ(H)) ‖GΘ(y)− x‖π∗(dx, dy)
= max ϕ ∫ X ϕ(x)[1− ρ(x)]µ(dx) + ∫ Y ϕc(GΘ(y))[1− σ(y)]ν(dy)
= max ϕ∈L1(µ) ∫ X ϕ(x)[1− ρ(x)]µ(dx)− ∫ Y ϕ(GΘ(y))[1− σ(y)]ν(dy)
for the non-restricted marginals µ and ν, where for the last equality we use ϕc = −ϕ for p = 1 when ϕ is 1-Lipschitz (Villani, 2008). Using the same technique, the second term in (14) can be computed as ∫
X×Y\(Γ(GΘ)∪Γ(H)) ‖y −Hx‖dπ(x, y)
= max ψ∈L1(ν) ∫ Y ψ(y)[1− σ(y)]dν(y)− ∫ X ψ(Hx)[1− ρ(x)]dµ(x).
By collecting terms, we conclude the proof.

heading: C OTHER VARIANTS OF UNSUPERVISED LEARNING
text: Here, we derive the dual optimal transport formulation when the PLS cost with deep learning prior as in (8) is used as the transportation cost for the OT problem. Specifically, the transportation cost of the OT problem with respect to the DIP cost in (8) becomes
min Θ
c(y, z; Θ,H) = ‖y −HQΘ(z)‖ (27)
where we consider y and z as random variable and we use non-squared norm for simplicity. Then, the corresponding OT problem can be dualized as follows:
min Θ min π ∫ X×Z ‖y −HGΘ(z)‖dπ(x, z) = min Θ max ψ∈L1(ν) ∫ Y ψ(y)dν(y)− ∫ X ψ(HGΘ(z))dη(z).
Although this is a nice way of pretraining deep image prior model using unmatched training data set, the final image estimate still comes from the following optimization problem:
x = GΘ(z ∗) where z∗ = arg min z ‖y −HGΘ∗(z)‖
where Θ∗ is the estimated network parameters from previous training step. This is equivalent to deep generator model (Bora et al., 2017). Therefore, this is not useful in designing a feed-forward neural network in an unsupervised manner.

heading: D EXPERIMENTAL DETAILS

heading: D.1 ACCELERATED MRI
text: We use dataset for fastMRI challenge (Zbontar et al., 2018) for our experiments. This dataset is composed of MR images of knees. We extracted 3500 MR images from fastMRI single coil validation set. Then, 3000 slices are used for training/evaluation, and 500 slices are used for test. These MR images are fully sampled images, so we make undersampled images by a randomly subsampling k-space lines. The acceleration factor is four, and autocalibration signal (ACS) region contains 4% of k-space lines. Each slice is normalized by standard deviation of the magnitude of each slice. To handle complex values of data, we concatenate real and imaginary values along the channel dimension. Each slice has different size, so we use only single batch. The images are center-cropped to the size of 320× 320, and then the peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM) values are calculated.
We use U-Net generator to reconstruct fully sampled MR images from undersampled MR images as shown in Fig. 6(a). Our generator consists of 3 × 3 convolution, Instance normalization, and leaky ReLU operation. Also, there are skip-connection and pooling layers. At the last convolution layer, we do not use any operation. Our discriminator is same as the discriminator of original CycleGAN. We use PatchGAN (Isola et al., 2017), so the discriminator classifies inputs at patch scales. The discriminator also consists of convolution layer, instance normalization, and leaky ReLU operation as shown in Fig. 6(b). We use Adam optimizer to train our network, with momentum parameters β1 = 0.5, β2 = 0.9, and learning rate of 0.0001. The discriminator is updated 5 times for every generator updates. We use batch size of 1, and trained our network during 100 epochs. Our code was implemented by TensorFlow.

heading: D.2 DECONVOLUTION MICROSCOPY
text: The network architecture of the high resolution image generator GΘ from the low-resolution image is a modified 3D-Unet (Çiçek et al., 2016) as shown in Fig 7. Our U-net structure consists of contracting and expanding paths. The contracting path consists of the repetition of the following blocks: 3D conv- Instance Normalization (Ulyanov et al., 2016)- ReLU. Here, the generator has symmetric configuration so that both encoder and decoder have the same number of layers, i.e. κ = 7. Throughout the network, the convolutional kernel dimension is 3 × 3 × 3. There exists a pooling layer and skipped connection for every other convolution operations. To enhance the image contrast, we add additional sigmoid layer at the end of U-Net. On the other hand, the low-resolution image generator from high resolution input is based on a single 3D convolution layer that models a 3D blurring kernel. The size of the 3D PSF modeling layer is chosen depending on the problem set by considering their approximate PSF sizes. In this paper, the size of the 3D PSF layer is set to 20× 20× 20. As for the discriminators, we follow the original CycleGAN that uses multi-PatchGANs (mPGANs) (Isola et al., 2017), where each discriminator has input patches with different sizes used. As shown in Fig 8, it consist of three independent discriminators. Each discriminator takes patches at different sizes: original, and half, and quarter size patches.

heading: D.2.1 REAL EXPERIMENT
text: A total 18 epi-fluorescent (EPF) microscopy images of tubulin with a size of 512 × 512 × 30 were used for training, and one for validation. As for unmatched sharp image volume, we used deblurred image generated by utilizing a commercial software AutoQuant X3 (Media Cybernetics, Rockville). The EPF volume depth was increased to 64 using the reflection boundary condition. Due to GPU memory limitations, the EPF volume was split into 64×64×64 size patches. For data augmentation, rotation, flip, translation, and scale were imposed on the input patches. We normalized the patches and set them to [0,1]. Adam optimizer (Kingma & Ba, 2014) was also used for training. The learning rate was initially set to 0.0001, which is decreased linearly after 40 epoch, and the total number of epoch was 200 epoch. For the optimizer, we used only a single batch.

heading: D.2.2 SIMULATION STUDY
text: For simulation studies with the ground-truth data, we used synthetic microtubule network data set (Sage et al., 2017) to train and validate our model. Specifically, from the ground-truth high resolution synthetic microtubule images, we generate blurred images by convolving them with a model PSF. In particular, the numerical PSF was computed using the Born and Wolf model (Kirshner et al., 2011), which is given by
h(rx, ry, rz) = ∣∣∣C ∫ 1
0
J0
[ k NA ni ρ √ r2x + r 2 y ] e − 12 jkρ 2rz ( NA ni )2 ρdρ ∣∣∣2 (28)
where C is a normalizing constant, k = 2π/λ is the wave number of emitted light, λ is the wavelength, NA is the numerical aperture, and ni is the refractive index of immersion layer. For all simulation, we use NA = 1.4, and ni = 1.5. Then, the simulation is performed with different λ for the PSF model (28). The convolution was performed using DeconvolutionLab2 (Sage et al., 2017).
The convolved data was then added with the mixture of zero mean Gaussian with the standard deviation σ = 3 and Poisson noises with the parameter m = 7. The size of the synthetic data was 256 × 512 × 128, and 18 samples of the synthetic data were used for training and the other one is used for validation. Due to memory limitations, the volume was split into 64×64×64 size patches, which are used as inputs to the neural network.
Fig. 9 shows qualitative comparisons of different reconstruction methods over datasets generated by the theoretical PSFs at different wavelengths. Here, the training was performed using the blurred data with λ = 500nm PSF kernnel, but the inference was performed using data with wavelength 400 nm, 500 nm, and 600 nm, respectively. The goal of this study is to investigate the generalization power of the proposed method with respect to different wavelengths.
The supervised method removed most blur and noise; however, it lost many fine details of the tubule structure and showed false continuity over the tubule structure. The conventional cycleGAN method showed better qualitative results than the supervised one, presenting finer details and the correct continuity of the tubule structure, but still some of the structures were not fully recovered. On the other hand, the proposed method showed the best qualitative results, recovering most of the microtubule structures. Moreover, as shown in Table 2, the proposed method showed the best f and SSIM scores. This confirmed that the proposed method generalizes well, even when the PSF model for the training data is not matched well with that of the inference phase.

heading: E COMPARISON OF THE CLASSICAL AND THE PROPOSED PENALTY FOR INVERSE PROBLEMS
text: Recall that the inverse problem y = Hx have many solutions due to the ill-poseness of the mapping H. The main motivation of the classical PLS presented in Section 2.2 is to introduce the penalty for choosing the solution x, which is most likely based on the prior distribution of data (see Fig. 10(a)). On the other hand, unsupervised learning with training data provides an another way to resolve the ambiguity of the feasible solutions. Specifically, if we define a single-valued function GΘ(y) and
impose the constraint x = GΘ∗(y) with the learned parameter Θ∗, many of the feasible solutions for y = Hx can be cut out as shown in Fig. 10(b). Now, the remaining question is how to determine the optimal parameter Θ∗. This is where the optimal transport theory comes into play. Specifically, our OT formulation tries to find the parameter Θ∗ by minimizing the average transportation cost c(x, y;H,Θ∗) so that it eventually approaches zero, the primal and the dual solution may become equivalent, and GΘ∗ can be an inverse of the forward operator H. This is the main motivation for using the new penalty function in our formulation.

Review summaries:
The paper presents an interesting connection between cycleGAN, penalized Least Squares (PLS) and optimal transport (OT). From the PLS with deep regularization formulation, the authors cast a general OT problem with a specific cost combining measurement and reconstruction errors (eq. 9). From this formulation, the problem is expressed as a combination of a cycle consistency loss and an OT loss. The formulation is more generic than the classical cycleGAN formulation. Qualitative experiments are conducted on two different scenarii: accelerated MRI and deconvolution microscopy, for which the proposed method achieves good performances.