Title: OPTIMAL TRANSPORT, CYCLEGAN, AND PENALIZED LS FOR UNSUPERVISED LEARNING IN INVERSE PROB-LEMS

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 cy-cleGAN 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.

1. INTRODUCTION:

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

EQUATION

where w is the measurement noise, and H : X → Y is the measurement operator. In inverse problems originating from physics, the measurement operator is usually represented by an integral equation:

EQUATION

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)

for q ≥ 1, where R(x) is a regularization (or penalty) function (l 1 , 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 min x 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; Peyré 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 cycle-GAN 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.

2. RELATED WORKS:

2.1 OPTIMAL TRANSPORT (OT)

OT compares two measures in a Lagrangian framework (Villani, 2008; Peyré et al., 2019) . Formally, we say that T :

EQUATION

which is often simply represented by ν = T # µ, where T # is often called the push-forward operator. Monge's original OT problem (Villani, 2008; Peyré 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; Peyré 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; Peyré 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)| ≤ c X (x) + c Y (y) for some c X ∈ L 1 (µ) and c Y ∈ L 1 (ν)

, where L 1 (µ) denotes the set of 1-Lipschitz functions with the measure µ. Then, there is duality:

min πΠ(µ,ν) X ×Y c(x, y)dπ(x, y) = max ϕ∈L 1 (µ) 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) = sup x (c(x, y) -ϕ(x)) is properly defined.

2.2. PLS WITH DEEP LEARNING PRIOR:

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:

EQUATION

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 Θ.

3. MAIN CONTRIBUTIONS:

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

EQUATION

If this holds for all y, then G Θ is an inverse operator of H. 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):

EQUATION

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 → X is single-valued, then the average transportation cost T(Θ, H) in (11) can be decomposed as

EQUATION

where

EQUATION

EQUATION

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

EQUATION

EQUATION

where ϕ and ψ are 1-Lipschitz Kantorovich potential functions. 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:

EQUATION

where

EQUATION

where cycle (Θ, H) denotes the cycle-consistency loss in (15) and GAN (Θ, h; Φ, Ξ) is the GAN loss given by:

EQUATION

+ 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:

EQUATION

Note that our formulation using ( 19) with ( 20) and ( 21) has many unique features compared to the standard cycleGAN (Zhu et al., 2017) .

• 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 (Peyré 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.

4. APPLICATIONS TO INVERSE PROBLEMS:

4.1 ACCELERATED MRI

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

EQUATION

)

where F is the 2-D Fourier transform and P Ω 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 -1 P Ω Fx + F -1 w (23)

This results in the following cost function for the PLS formulation:

EQUATION

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 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.

4.2. DECONVOLUTION MICROSCOPY:

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

EQUATION

)

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 Here, PSNR = 20 log 10 (n

x * ∞ / x -x * 2 ) and SSIM = (2µ x µ x * + c 1 )(2σ xx * + c 2 )/(µ 2 x + µ 2 x * + c 1 )(σ 2 x + σ 2 x * + c 2 )

, where x and x * denote the reconstructed images and ground truth, respectively, n is the number of pixels, µ x is an average of x, σ 2

x is a variance of x and σ xx * is a covariance of x and x * , and c 1 , c 2 are two variables to stabilize the division. (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. 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, 

5. CONCLUSIONS:

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.

A ADDITIONAL THEOREMS AND LEMMA:

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 proj X (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 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 proj X (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)

B PROOF OF MAIN RESULTS:



B.1 PROOF OF LEMMA 1:

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.

B.2 PROOF OF PROPOSITION 1:

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) p dπ * (x, y)+ Γ(GΘ)

y-HG Θ (y) q dπ * (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) p dπ * (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 1 y=Hx π * (dy|x)µ(dx) = X x -G Θ (Hx) p ρ(x)µ(dx)

where 1 S denotes the indicator function for the set S and ρ(x) = 1 y=Hx π * (dy|x). In a similar fashion, with σ(y) = 1 x=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.

B.3 PROOF OF PROPOSITION 2:

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 l OT (Θ, 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 ϕ∈L 1 (µ) 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 ψ∈L 1 (ν) Y ψ(y)[1 -σ(y)]dν(y) - X ψ(Hx)[1 -ρ(x)]dµ(x).

By collecting terms, we conclude the proof.

C OTHER VARIANTS OF UNSUPERVISED LEARNING:

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

EQUATION

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 ψ∈L 1 (ν) 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. 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. 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) 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

EQUATION

where C is a normalizing constant, k = 2π/λ is the wave number of emitted light, λ is the wavelength, N A is the numerical aperture, and n i is the refractive index of immersion layer. For all simulation, we use N A = 1.4, and n i = 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.

E COMPARISON OF THE CLASSICAL AND THE PROPOSED PENALTY FOR INVERSE PROBLEMS:

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 . 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.

Figure 10 : Two strategies for resolving ambiguities in the feasible solutions in an ill-posed inverse problem. (a) Classical PLS approach using a close form prior distribution, and (b) our PLS approach using an inverse mapping to define a prior.

:

• 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.

Figures and Tables:

Figure 1: (a) Supervised learning with matched references, (b) OT with absolute continuous measures, and (c) OT having measures with singularities. CycleGAN corresponds to scenario (c).

Figure 2: Proposed cycleGAN architecture for accelerated MRI with 1-D downsampling patterns.

Figure 3: Unsupervised learning results for accelerated MRI using proposed cycleGAN. The values in the corners are PSNR/SSIM values for each image.Here, PSNR = 20 log 10 (nx * ∞ / x -x * 2 ) and SSIM = (2µ x µ x * + c 1 )(2σ xx * + c 2 )/(µ 2 x + µ 2 x * + c 1 )(σ 2 x + σ 2 x * + c 2 ), where x and x * denote the reconstructed images and ground truth, respectively, n is the number of pixels, µ x is an average of x, σ 2x is a variance of x and σ xx * is a covariance of x and x * , and c 1 , c 2 are two variables to stabilize the division.

formulation: c(x, y; Θ, h) = y -h * x + G Θ (y) -x .

Figure 4: Proposed cycleGAN architecture with a blur kernel for deconvolution microscopy. Here, G Θ denotes the CNN-based low-resolution to high-resolution generator. The blur generator is composed of a linear blur kernel h.

Fig.5shows 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,

Figure 5: Comparison of reconstruction results by various methods: (a) blurred image measurements, (b) classical deconvolution method, (c) supervised learning, (d) conventional cycleGAN, and (e) proposed cycleGAN. The regions of interest (marked yellow) show enlargements of the areas marked in yellow.

, y)dπ(x, y)

Figure 6: Proposed network architecture for (a) generator and (b) discriminator in accelerated MRI.

Figure 7: A modified 3D U-net architecture for our high-resolution image generator in the deconvolution microscopy problem.

Figure 8: Multi-PatchGANs discriminator architecture in the deconvolution microscopy problem.

Figure 9: Generalization performance comparison of various methods. The following parameters were used for the training and inference: N A = 1.4, n i = 1.5, and σ = 3, m = 7. For the training data, the wavelength was λ = 500nm. (a) Reconstruction result at the inference phase using the data with λ = 400nm. (b) Reconstruction result at the inference phase using the data with λ = 500nm. (c) Reconstruction result at the inference phase using the data with λ = 600nm. Two different views (YZ, XZ) were displayed along the corresponding lines. The ROIs (marked yellow) show the area for the enlarged parts.

Quantitative comparison for various algorithms on 500 test sets of fastMRI data.

Table content: <table><tr><td>Metric</td><td colspan="4">Input (×4) Supervised Learning Conventional CycleGAN Proposed Method</td></tr><tr><td>PSNR (dB)</td><td>26.8056</td><td>27.9609</td><td>25.0092</td><td>28.1704</td></tr><tr><td>SSIM</td><td>0.6419</td><td>0.6419</td><td>0.4689</td><td>0.6550</td></tr></table>

, 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. Performance comparison of various methods of Fig. 9 in terms of PSNR and SSIM.

Table content: <table><tr><td>D.2.1 REAL EXPERIMENT</td></tr><tr><td>A total 18 epi-fluorescent (EPF) microscopy images of tubulin with a size of 512 × 512 × 30 were</td></tr><tr><td>used for training, and one for validation. As for unmatched sharp image volume, we used deblurred</td></tr><tr><td>image generated by utilizing a commercial software AutoQuant X3 (Media Cybernetics, Rockville).</td></tr></table>

Peer review (iclr2020):
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. In general I like the paper and the corresponding idea. However, I had a hard time understanding some of the key elements of the proposed method. Notably, in lemma 1, it is hard to understand in which case mathcal{X} 	imes mathcal{X} / A union B is not an empty set. It is clearly the case whenever G or H are bijections, or whenever an exact reconstruction is achievable, but I guess a more open discussions could be conducted. Also, the definition of A and B are clearly dependent of G and H and this could be highlighted in the notation. As noted in proposition 2, however, the corresponding OT loss can be computed on the entire subset X and Y (since the Kantorovich potentials match those obtained on the restricted set). As a consequence, we can see the overall training procedure when H is fixed (as it is the case in both experiments), as learning for a generator in a WGAN way, that also enforces a consistency constraint (that can be seen as a regularization of the OT problem). While I like this idea, I think there might be some better justification for it (it all starts with Eq. 9 and the choice of this model. Why enforcing this minimal cost equation ?). Note that the model imposes that p=q=1, which is somehow limited. I guess other choices of p and q could be possible, eventually using regularized version of OT to remove the hard constraints on the kantorovich potentials I am willing to revise my note positively provided that a sufficient number of convincing explanations are given on my previous remarks. Minor remark Eq. 8 should include a weighting factor between the two terms since the dimensions are not the same. In the proof of Lemma 1, since you derive an equality for T, do you have to establish the first inequality ? Prop. 1 I guess there is an error on the first line (the integral is not over mathcal{X} 	imes mathcal{X} / A union B given the previous definition Missing related work DLOW: Domain Flow for Adaptation and Generalization Rui Gong, Wen Li, Yuhua Chen, Luc Van Gool; The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019, pp. 2477-2486