Forward Operator Estimation in Generative Models with Kernel Transfer Operators
Year: 2021

Abstract:
Generative models which use explicit density modeling (e.g., variational autoencoders, flow-based generative models) involve finding a mapping from a known distribution, e.g. Gaussian, to the unknown input distribution. This often requires searching over a class of non-linear functions (e.g., representable by a deep neural network). While effective in practice, the associated runtime/memory costs can increase rapidly, usually as a function of the performance desired in an application. We propose a much cheaper (and simpler) strategy to estimate this mapping based on adapting known results in kernel transfer operators. We show that our formulation enables highly efficient distribution approximation and sampling, and offers surprisingly good empirical performance that compares favorably with powerful baselines, but with significant runtime savings. We show that the algorithm also performs well in small sample size settings (in brain imaging).

heading: 1 Introduction
text: Generative modeling, in its unconditional form, refers to the problem of estimating the data generating distribution: given i.i.d. samples X with an unknown distribution PX , a generative model seeks to find a parametric distribution that closely resembles PX . In modern deep generative models, we often approach this problem via a latent variable – i.e., we assume that there is some variable Z ∈ Z associated with the observed data X ∈ X that follows a known distribution PZ (also referred to as the prior in generative models). Thus, we can learn a mapping f : Z → X such that the distribution after transformation, denoted by Pf(Z), aligns well with the data generating distribution PX . Therefore, sampling from PX becomes convenient since PZ can be efficiently sampled. Frequently, f is parameterized by deep neural networks and optimized with stochastic gradient descent (SGD).
Existing generative modeling methods variously optimize the transformation f , most commonly modeling it as a Maximum Likelihood Estimation (MLE) or distribution matching problem. For instance, given data X = {x1, . . . ,xn}, a variational autoencoder (VAE) (Kingma & Welling, 2013) first constructs Z through the approximate posterior qZ|X and maximizes a lower bound of likelihood pf(Z)(X). Generative adversarial networks (GANs) (Goodfellow et al., 2014) relies on a simultaneously learned discriminator such that samples of Pf(Z) are indistinguishable from X. Results in (Arjovsky et al., 2017; Li et al., 2017) suggest that GANs minimize the distributional discrepancies between Pf(Z) and PX . Flow-based generative models optimize pf(Z)(X) explicitly through the change of variable rule and efficiently calculating the Jacobian determinant of the inverse mapping f−1.
In all examples above, the architecture or objective notwithstanding, the common goal is to find a suitable function f that reduces the difference between Pf(Z) and PX . Thus, a key component in many deep generative models is to learn a forward operator as defined below.
Definition 1.1 (Forward operator). A forward operator f? ∈ C : Z → X is defined to be a mapping associated with some latent variable Z ∼ PZ such that f? = arg minf∈C d(Pf(Z), PX) for some function class C and a distance measure d(·, ·).
Motivation: The specifics of the forward operator may differ from case to case. But its properties and how it is estimated numerically greatly influences the empirical performance
of the model. For instance, mode collapse issues in GANs are well known and solutions continue to emerge (Srivastava et al., 2017). To learn the forward operator, VAEs use an approximate posterior qZ|X that may sometimes fail to align with the prior (Kingma et al., 2016; Dai & Wipf, 2019). Flow-based generative models enable direct access to the posterior likelihood, yet in order to tractably evaluate the Jacobian of the transformation during training, one must either restrict the expressiveness at each layer (Dinh et al., 2017; Kingma & Dhariwal, 2018) or use more involved solutions (Chen et al., 2018). Of course, solutions to mitigate these weaknesses (Ho et al., 2019) remains an active area of research.
The starting point of our work is to evaluate the extent to which we can radically simplify the forward operator in deep generative models. Consider some desirable properties of a hypothetical forward operator (in Def. (1.1)): (a) Upon convergence, the learned operator f? minimizes the distance between PX and Pf(Z) over all possible operators of a certain class. (b) The training directly learns the mapping from the prior distribution PZ , rather than a variational approximation. (c) The forward operator f? can be efficiently learned and sample generation is also efficient. It would appear that these criteria violate the “no free lunch rule”, and some compromise must be involved. Our goal is to investigate this trade-off: which design choices can make this approach work? Specifically, a well studied construct in dynamical systems, namely the Perron-Frobenius operator (Lemmens & Nussbaum, 2012), suggests an alternative linear route to model the forward operator. Here, we show that if we are willing to give up on a few features in existing models – this may be acceptable depending on the downstream use case – then, the forward operator in generative models can be efficiently approximated as the estimation of a closed-form linear operator in the reproducing kernel Hilbert space (RKHS). With simple adjustments of existing results, we identify a novel way to replace the expensive training for generative tasks with a simple principled kernel approach.
Contributions. Our results are largely based on results in kernel methods and dynamical systems, but we demonstrate their relevance in generative modeling and complement recent ideas that emphasize links between deep generative models and dynamical systems. Our contributions are (a) We propose a non-parametric method for transferring a known prior density linearly in RKHS to an unknown data density – equivalent to learning a nonlinear forward operator in the input space. When compared to its functionally-analogous module used in other deep generative methods, our method avoids multiple expensive training steps yielding significant efficiency gains; (b) We evaluate this idea in multiple scenarios and show competitive generation performance and efficiency benefits with pre-trained autoencoders on popular image datasets including MNIST, CIFAR-10, CelebA and FFHQ; (c) As a special use case, we demonstrate the advantages over other methods in limited data settings.
2 Preliminaries
We briefly introduce reproducing kernel Hilbert space (RKHS) and kernel embedding of probability distributions, concepts we will use frequently.
Definition 2.1 (RKHS (Aronszajn, 1950)). For a set X , let H be a set of functions g : X → R. Then, H is a reproducing kernel Hilbert space (RKHS) with a product 〈·, ·〉H if there exists a function k : X × X → R (called a reproducing kernel) such that (i) ∀x ∈ X , g ∈ H, g(x) = 〈g, k(x, ·)〉H; (ii) H = cl(span({k(x, ·), x ∈ X})), where cl(·) is the set closure.
The function φ(x) = k(x, ·) : X → H is referred to as the feature mapping of the induced RKHS H. A useful identity derived from feature mappings is the kernel mean embedding: it defines a mapping from a probablity measure in X to an element in the RKHS. Definition 2.2 (Kernel Mean Embedding (Smola et al., 2007)). Given a probability measure p on X with an associated RKHS H equipped with a reproducing kernel k such that supx∈X k(x, x) < ∞, the kernel mean embedding of p in RKHS H, denoted by µp ∈ H, is
defined as µp = Ep[φ(x)] = ∫ k(x, ·)p(x)dx, and the mean embedding operator E : L1(X )→ H is defined as µp = Ep.
Remark 1. For characteristic kernels, the operator E is injective. Thus, two distributions (p, q) in X are identical iff Ep = Eq.
This property allows using of Maximum Mean Discrepancy (MMD) for distribution matching (Gretton et al., 2012; Li et al., 2017) and is common, see (Muandet et al., 2017; Zhou et al., 2018). For a finite number of samples {xi}ni=1 drawn from the probability measure p, an unbiased empirical estimate of µH is µ̂H = 1 n ∑n i=1 k(xi, ·) such that limn→∞ 1 n ∑n i=1 k(xi, ·) = µH.
Next, we review the covariance/cross covariance operators, two widely-used identities in kernel methods (Fukumizu et al., 2013; Song et al., 2013) and building blocks of our approach.
Definition 2.3 (Covariance/Cross-covariance Operator). Let X,Z be random variables defined on X × Z with joint distribution PX,Z and marginal distributions PX , PZ . Let (l, φ,H) and (k, ψ,G) be two sets of (a) bounded kernel, (b) their corresponding feature map, and (c) their induced RKHS, respectively. The (uncentered) covariance operator CZZ : H → H and cross-covariance operator CXZ : H → G are defined as
CZZ , Ez∼PZ [φ(z)⊗ φ(z)] CXZ , E(x,z)∼PX,Z [ψ(x)⊗ φ(z)] (1)
where ⊗ is the outer product operator.

heading: 3 Simplifying the estimation of the forward operator
text: Forward operator as a dynamical system: The dynamical system view of generative models has been described by others (Chen et al., 2018; Grathwohl et al., 2019; Behrmann et al., 2019). These strategies model the evolution of latent variables in a residual neural network in terms of its dynamics over continuous or discrete time t, and consider the output function f as the evaluation function at a predetermined boundary condition t = t1. Specifically, given an input (i.e., initial condition) z(t0), f is defined as
f(z(t0)) = z(t0) + ∫ t1 t0 ∆t(z(t))dt (2)
where ∆t is a time-dependent neural network function and z(t) is the intermediate solution at t. This view of generative models is not limited to specific methods or model archetypes, but generally useful, for example, by viewing the outputs of each hidden layer as evaluations in discrete-time dynamics. After applying f on a random variable Z ∈ Z, the marginal density of the output over any subspace Λ ⊆ X can be expressed as∫
Λ
pf(Z)(x)dx = ∫ z∈f−1(Λ) pZ(z)dz (3)
If there exists some neural network instance ∆?t such that the corresponding output function f? satisfies PX = Pf?(Z), by Def. 1.1, f
? is a forward operator. Let X be a set of i.i.d. samples drawn from PX . In typical generative learning, either maximizing the likelihood 1|X| ∑ x∈X pf(Z)(x) or minimizing the distributional divergence d(Pf(Z), PX) requires evaluating and differentiating through f or f−1 many times.
Towards a one-step estimation of forward operator: Since f and f−1 in (3) will be highly nonlinear in practice, evaluating and computing the gradients can be expensive. Nevertheless, the dynamical systems literature suggests a linear extension of f∗, namely the Perron-Frobenius operator or transfer operator, that conveniently transfers pZ to pX .
Definition 3.1 (Perron-Frobenius operator (Mayer, 1980)). Given a dynamical system f : X → X , the Perron-Frobenius (PF) operator P : L1(X ) → L1(X ) is an infinitedimensional linear operator defined as ∫ Λ (PpZ)(x)dx = ∫ z∈f−1(Λ) pZ(z)dz for all Λ ⊆ X .
Although in Def. 3.1, the PF operator P is defined for self-maps, it is trivial to extend P to mappings f : Z → X by restricting the RHS integral ∫ z∈f−1(Λ) pZ(z)dz to Z.
It can be seen that, for the forward operator f∗, the corresponding PF operator P satisfies pX = PpZ . (4) If P can be efficiently computed, transferring the tractable density pZ to the target density pX can be accomplished simply by applying P. However, since P is an infinite-dimensional operator on L1(X ), it is impractical to instantiate it explicitly and exactly. Nonetheless, there exist several methods for estimating the Perron-Frobenius operator, including Ulam’s method (Ulam, 1960) and the Extended Dynamical Mode Decomposition (EDMD) (Williams et al., 2015a). Both strategies project P onto a finite number of hand-crafted basis functions – this may suffice in many settings but may fall short in modeling highly complex dynamics.
Kernel-embedded form of PF operator: A natural extension of PF operator is to represent P by an infinite set of functions (Klus et al., 2020), e.g., projecting it onto the bases of an RKHS via the kernel trick. There, for a characteristic kernel l, the kernel mean embedding uniquely identifies an element µX = ElpX ∈ G for any pX ∈ L1(X ). Thus, to approximate P , we may alternatively solve for the dynamics from pZ to pX in their embedded form. Using Tab. 1 notations, we have the following linear operator that defines the dynamics between two embedded densities.
Definition 3.2 (Kernel-embedded Perron-Frobenius operator (Klus et al., 2020)). Given pZ ∈ L1(X ) and pX ∈ L1(X ). Denote k as the input kernel and l as the output kernel. Let µX = ElpX and µZ = EkpZ be their corresponding mean kernel embeddings. The kernel-embedded Perron-Frobenius (kPF) operator, denoted by PE : H → G, is defined as PE = CXZC−1ZZ (5) Proposition 3.1 (Song et al. (2013)). With the above definition, PE satisfies µX = PEµZ (6) under the conditions: (i) CZZ is injective (ii) µt ∈ range(CZZ) (iii) E[g(X)|Z = ·] ∈ H for any g ∈ G.
The last two assumptions can sometimes be difficult to satisfy for certain RKHS (see Theorem 2 of Fukumizu et al. (2013)). In such cases, a relaxed solution can be constructed by replacing C−1ZZ by a regularized inverse (CZZ + λI)−1 or a Moore-Penrose pseudoinverse C † ZZ .
The following proposition shows commutativity between the (kernel-embedded) PF operator and the mean embedding operator, showing its equivalence to P when l is characteristic. Proposition 3.2 ((Klus et al., 2020)). With the above notations, El ◦ P = PE ◦ Ek.
Transferring embedded densities with the kPF operator: The kPF operator is a powerful tool that allows transferring embedded densities in RKHS. The main steps are:
(1) Use mean embedding operator El on pZ . Let us denote it by µZ . (2) Transfer µZ using kPF operator PE to get the mean embedded pX , given by µX .
Of course, in practice with finite data, {xi}i∈[n] ∼ PX and {zi}i∈[n] ∼ PX , PE must be estimated empirically (see Klus et al. (2020) for an error analysis).
P̂E = ĈXZ(ĈZZ)−1 ≈ Ψ(ΦTΦ + λnI)−1ΦT ≈ Ψ(ΦTΦ)†ΦT
where Φ = [k(z1, ·), · · · , k(zn, ·)], Ψ = [l(x1, ·), · · · , l(xn, ·)] are simply the corresponding feature matrices for samples of PX and PZ , and λ is a small penalty term.
Learning kPF for unconditional generative modeling: Some generative modeling methods such as VAEs and flow-based formulations explicitly model the latent variable Z as conditionally dependent on the data variable X. This allows deriving/optimizing the likelihood pf(Z)(X). This is desirable but may not be essential in all applications. To learn a kPF, however, X and Z can be independent RVs. While it may not be immediately obvious why we could assume this independence, we can observe the following property for the empirical kPF operator, assuming that the empirical covariance operator ĈZZ is non-singular:
P̂E µ̂Z = ĈXZ Ĉ−1ZZ µ̂Z = ΨΦ >︸ ︷︷ ︸
ĈXZ (ΦΦ>︸ ︷︷ ︸ ĈZZ )−1Φ1n = Ψ(Φ >Φ)−1Φ>Φ1n = Ψ1n = µ̂X (7)
Suppose that {xi}i∈[n] and {zj}i∈[n] are independently sampled from the marginals PX and PZ . It is easy to verify that (7) holds for any pairing {(xi, zj)}(i,j)∈[n]×[n]. However, instantiating the RVs in this way rules out the use of kPF for certain downstream tasks such as controlled generation or mode detection, since Z does not contain information regarding X. Nevertheless, if sampling is our only goal, then this instantiation of kPF will suffice.
Mapping Z to G: Now, since PE is a deterministic linear operator, we can easily set up a scheme to map samples of Z to elements of G where the expectation of the mapped samples equals µX
Define φ(z) = k(z, ·) and ψ(x) = l(x, ·) as feature maps of kernels k and l. We can rewrite µX as
µX = PEEkpZ = PEEZ [φ(Z)] = EZ [PE(φ(Z))] = EZ [ψ ( ψ−1 (PEφ (Z)) ) ] (8)
Here ψ−1 is the inverse or the preimage map of ψ. Such an inverse, in general, may not exist (Kwok & Tsang, 2004; Honeine & Richard, 2011). We will discuss a procedure to approximate ψ−1 in §4.1. In what follows, we will temporarily assume that an exact preimage map exists and is tractable to compute.
Define Ψ∗ = P̂Eφ(Z) as the transferred sample in G using the empirical embedded PF operator P̂E . Then the next result shows that asymptotically the transferred samples converge in distribution to the target distribution.
Proposition 3.3. As n → ∞, ψ−1 (Ψ∗) D→ PX . That is, the preimage of the transferred sample approximately conforms to PX under previous assumptions when n is large.
Proof. Since P̂E asymp.→ P, the proof immediately follows from (8).

heading: 4 Sample generation using the Kernel transfer operator
text: At this point, the transferred sample Ψ∗, obtained by the kPF operator, remains an element of RKHS G. To translate the samples back to the input space, we must find the preimage x∗ such that ψ(x∗) = Ψ∗.

heading: 4.1 Solving for an approximate preimage
text: Solving the preimage in kernel-based methods is known to be ill-posed (Mika et al., 1999) because the mapping ψ(·) is not necessarily surjective, i.e., a unique preimage x∗ = ψ−1(Ψ∗),Ψ∗ ∈ H may not exist. Often, an approximate preimage ψ†(X,Ψ∗) ≈ ψ−1(Ψ∗) is constructed instead based on relational properties among the training data in the RKHS. We consider two options in our framework (1) MDS-based method (Kwok & Tsang, 2004; Honeine & Richard, 2011),
ψ†MDS(X,Ψ ∗) = 12 (X ′X′>)−1X′(diag(X′>X′)− d>),where ∀i ∈ [γ],di = ‖l(x′i, ·)−Ψ∗‖G (9) which optimally preserves the distances in RKHS to the preimages in the input space, and (2) weighted Fréchet mean (Friedman et al., 2001), which in Euclidean space takes the form
ψ†wFM(X,Ψ ∗) = ψ†wFM(X ′; s) = X′s/‖s‖1,where ∀i ∈ [γ], si = 〈l(x′i, ·),Ψ∗〉 (10)
Algorithm 1 Sample Generation from kPF
1: Input: Training data X = {x1, . . . ,xn}, Optional autoencoder (E,D), input/output kernels (k, l), neighborhood size γ 2: Training 3: X = (E (x1) , . . . , E (xn)) if E is provided 4: Sample {zi}i∈[n] ∼ PZ independently 5: Construct L,K ∈ Rn×n s.t. Lij = l(xi,xj),Kij = k (zi, zj) 6: Kinv = (K + λnI) −1 or K† 7: Inference 8: Generate new prior sample z∗ ∼ PZ 9: s = L ·Kinv[k(z1, z∗) . . . k(zn, z∗)]>
10: ind = argsort(s)[−γ :] 11: x∗ = ψ†wFM (X[ind] ; s[ind] ). 12: Output D(x∗) if D is provided else x∗
where X′ a neighborhood of γ training samples based on pairwise distance or similarity in RKHS, following (Kwok & Tsang, 2004). The weighted Fréchet mean preimage uses the inner product weights 〈Ψ∗, ψ(xi)〉 as measures of similarities to interpolate training samples. On the toy data (as in Fig. 2), weighted Fréchet mean produces fewer samples that deviate from the true distribution and is easier to compute. Based on this observation, we use the weighted Fréchet mean as the preimage module for all experiments that requires samples, while acknowledging that other preimage methods can also be substituted in.
With all the ingredients in hand, we now present an algorithm for sample generation using the kPF operator in Alg. 1. The idea is simple yet powerful: at training time, we construct the empirical kPF operator P̂E using the training data {xi}i∈[s] and samples of the known prior {zi}ı∈[n]. At test time, we will transfer new points sampled from PZ to feature maps in H, and construct their preimages as the generated output samples.

heading: 4.2 Image generation
text: Image generation is a common application for generative models (Goodfellow et al., 2014; Dinh et al., 2017). While our proposal is not image specific, constructing sample preimages in a high dimensional space with limited training samples can be challenging, since the space of images is usually not dense in a reasonably sized neighborhood. However, empirically images often lie near a low dimensional manifold in the ambient space (Seung & Lee, 2000), and one may utilize an autoencoder (AE) (E,D) to embed the images onto a latent space that represents coordinates on a learned manifold. If the learned manifold lies close to the true manifold, we can learn densities on the manifold directly (Dai & Wipf, 2019).
Therefore, for image generation tasks, the training data is first projected onto the latent space of a pretrained AE. Then, the operator will be constructed using the projected latent representations, and samples will be mapped back to image space with the decoder of AE. Our setup can be viewed analogously to other generative methods based on so called “ex-post” density estimation of latent variables (Ghosh et al., 2020). We also restrict the AE latent space to a hypersphere Sn−1 to ensure that (a) k(·, ·) and l(·, ·) are bounded and (b) the space is geodesically convex and complete, which is required by the preimage computation. To compute the weighted Fréchet mean on a hypersphere, we adopt the recursive algorithm in Chakraborty & Vemuri (2015) (see appendix D for details).

heading: 5 Experimental Results
text: Goals. In our experiments, we seek to answer three questions: (a) With sufficient data, can our method generate new data with comparable performance with other state-of-the-art generative models? (b) If only limited data samples were given, can our method still estimate the density with reasonable accuracy? (c) What are the runtime benefits, if any?
Datasets/setup. To answer the first question, we evaluate our method on standard vision datasets, including MNIST, CIFAR10, and CelebA, where the number of data samples is much larger than the latent dimension. We compare our results with other VAE variants (Two-stage VAE Dai & Wipf (2019), WAE Arjovsky et al. (2017), CV-VAE Ghosh et al. (2020)) and flow-based generative models (Glow Kingma & Dhariwal (2018), CAGlow Liu et al. (2019)) The second question is due to the broad use of kernel methods in small sample size settings. For this more challenging case, we randomly choose 100 training samples (< 1% of the full dataset) from CelebA and evaluate the quality of generation compared to other density approximation schemes. We also use a dataset of T1 Magnetic Resonance (MR) images from the Alzheimer’s Disease Neuromaging Initiative (ADNI) study.
Distribution transfer with many data samples. We evaluate the quality by calculating the Fréchet Inception Distance (FID) (Heusel et al., 2017) with 10K generated images from each model. Here, we use a pretained regularized autoencoder (Ghosh et al., 2020) with a latent space restricted to the hypersphere (denoted by SRAE) to obtain smooth latent representations. We compare our kPF to competitive end-to-end deep generative baselines (i.e. flow and VAE variants) as well as other density estimation models over the same SRAE latent space. For the latent space models, we experimented with Glow (Kingma &
Dhariwal, 2018), VAE, Gaussian mixture model (GMM), and two proposed kPF operators with Gaussian kernel (RBF-kPF) and NTK (NTK-kPF) as the input kernel. The use of NTK is motivated by promising results at the interface of kernel methods and neural networks (Jacot et al., 2018; Arora et al., 2020). Implementation details are included in the appendix.
Comparative results are shown in Table 2. We see that for images with structured feature spaces, e.g., MNIST and CelebA, our method matches other non-adversarial generative models, which provides evidence in support of the premise that the forward operator can be simplified. Further, we present qualitative results on all datasets (in Fig. 3), where we compare our kPF operator based model with other density estimation techniques on the latent space. Observe that our model generates comparable visual results as SRAEGlow.
Since kPF learns the distribution on a pre-trained AE latent space for image generation, using a more powerful AE can offer improvements in generation quality. In Fig. 4, we present representative images by learning our kPF on NVAE Vahdat & Kautz (2020) latent space, pre-trained on the FFHQ dataset. NVAE builds a hierarchical prior and achieves state-of-the-art generation quality among VAEs. We see that kPF can indeed generate high-quality and diverse samples with the help of NVAE encoder/decoder. In fact, any AE/VAE may be substituted in, assuming that the latent space is smooth.
Summary: When a sufficient number of samples are available, our algorithm performs as well as the alternatives, which is attractive given the efficient training. In Fig. 5, we present comparative result of FIDs with respect to the training time. Since kPF can be computed in closed-form, it achieves significant training efficiency gain compared to other deep generative methods while delivering competitive generative quality.
Distribution transfer with limited data samples. Next, we present our evaluations when only a limited number of samples are available. Here, each of the density estimators was trained on latent representations of the same set of 100 randomly sampled CelebA images, and 10K images were generated to evaluate FID (see Table 3). Our method outperforms Glow and VAE, while offering competitive performance with GMM. Surprisingly, GMM remains a strong baseline for both tasks, which agrees with results in Ghosh et al. (2020). However, note that GMM is restricted by its parametric form and is less flexible than our method (as shown in Fig 2).
VAE Glow GMM RBF-kPF NTK-kPF
59.3 77.0 39.6 40.6 40.9
Table 3: FID values for few samples setting density approximation on CelebA.
Learning from few samples is common in biomedical applications where acquisition is costly. Motivated by interest in making synthetic but statistic preserving data (rather than the real patient records) pub-
licly available to researchers (see NIH N3C Data Overview), we present results on generating high-resolution (160 × 196 × 160) brain images: 183 samples from group AD (diagnosed as Alzheimer’s disease) and 291 samples from group CN (control normals). For n = 474 d = 5017600, using our kernel operator, we can generate high-quality samples that are in-distribution. We present comparative results with VAEs. The generated samples in Fig. 6 clearly show that our method generates sharper images. To check if the results are also scientifically meaningful, we test consistency between group difference testing (i.e., cases versus controls differential tests on each voxel) on the real images (groups were AD and CN) and the same test was performed on the generated samples (AD and CN groups), using a FWER corrected two-sample t-test Ashburner & Friston (2000). The results (see Fig 6) show that while there is a deterioration in regions identified to be affected by disease (different across groups), many statistically-significant regions from tests on the real images are preserved in voxel-wise tests on the generated images.
Summary: We achieve improvements in the small sample size setting compared to other generative methods. This is useful in many data-poor settings. For larger datasets, our method still compares competitively with alternatives, but with a smaller resource footprint.

heading: 6 Limitations
text: Our proposed simplifications can be variously useful, but deriving the density of the posterior given a mean embedding or providing an exact preimage for the generated sample in RKHS is unresolved at this time. While density estimation from kPF has been partially addressed in Schuster et al. (2020b), finding the pre-image is often ill-posed. The weighted Fréchet mean preimage only provides an approximate solution and evaluated empirically, and due to the interpolation-based sampling strategy, samples cannot be obtained beyond the convex hull of training examples. Making Z and X independent RVs also limits its use for certain downstream task such as representation learning or semantic clustering. Finally, like other kernel-based methods, the sup-quadratic memory/compute cost can be a bottleneck on large datasets and kernel approximation (e.g. (Rahimi & Recht, 2009)) may have to be applied; we discuss this in appendix F.

heading: 7 Conclusions
text: We show that using recent developments in regularized autoencoders, a linear kernel transfer operator can potentially be an efficient substitute for the forward operator in some generative models, if some compromise in capabilities/performance is acceptable. Our proposal, despite its simplicity, shows comparable empirical results to other generative models, while offering efficiency benefits. Results on brain images also show promise for applications to highresolution 3D imaging data generation, which is being pursued actively in the community.

heading: A Choice of Kernel is relevant yet flexible
text: In some cases, one would focus on identifying (finite) eigenfuntions and modes of the underlying operator (Williams et al., 2015b; Brunton et al., 2016). But rather than finding certain modes that best characterize the dynamics, we care most about minimizing the error of the transferred density and therefore whether the span of functions is rich/expressive enough. In particular, condition (iii) in Proposition 3.2 requires an input RKHS spanned by sufficiently rich bases (Fukumizu et al., 2013). For this reason, the choice of kernel here cannot be completely ignored since it determines the family of functions contained in the induced RKHS.
To explore the appropriate kernel setup for our application, we empirically evaluate the effect of using several different kernels via a simple experiment on MNIST. We first train an autoencoder to embed MNIST digits on to a hypersphere S2, then generate samples from kPF by the procedure described by Alg. 1 using the respective kernel function as the input kernel k. Subplot (b) and (c) in Fig. 7 show the generated samples using Radial Basis Function (RBF) kernel and arc-cosine kernel, respectively. Observe that the choice of kernel has an influence on the sample population, and a kernel function with superior empirical behavior is desirable.
Motivated by this observation, we evaluated the Neural Tangent Kernel (NTK) (Jacot et al., 2018), a well-studied neural kernel in recent works. We use it for a few reasons, (a) NTK, in theory, corresponds to a trained infinitely-wide neural network, which spans a rich set of functions that satisfies the assumption. (b) For well-conditioned inputs (i.e., no duplicates) on hypersphere, the positive-definiteness of NTK is proved in (Jacot et al., 2018). Therefore, invertibility of the Gram matrix KZZ = Φ
TΦ is almost guaranteed if the prior distribution pZ is restricted on a hypersphere (c) NTK can be non-asymptotically approximated (Arora et al., 2019). (d) Unlike other parametric kernels such as RBF kernels, NTK is less sensitive to hyperparameters, as long as the number of units used is large enough (Arora et al., 2019). Subplot (d) of Fig. 7 shows that kPF learned with NTK as input kernel is able to generate samples that are more consistent with the data distribution. However, we should note that NTK is merely a convenient choice of kernel that requires less tuning, and is not otherwise central to our work. In fact, as shown in our experiment in Tab. 2, a well-tuned RBF kernel may also achieve a similar performance. Indeed, in practice, any suitable choice of kernel may be conveniently adopted into the proposed framework without major modifications.

heading: B Density Estimation with kPF Operator
text: The displayed transferred density with the kPF operator on toy data in Fig. 1 is approximated using the empirical kernel conditional density operator (CDO) (Schuster et al., 2020a), since there is currently no known methods that can exactly reconstruct density from the transferred mean embeddings. The marginalized transferred density pPE has the following form
pPE = C −1 ρ PEµZ = C−1ρ CXZC−1ZZµZ , (11)
where Cρ = Ey∼ρ[l(y, ·)⊗ l(y, ·)] is the covariance operator of a independent reference density ρ in G. The above density function is also an element of RKHS G, and therefore we can evaluate the density at x∗ by using the reproducing property pPE (x
∗) = 〈pPE , l(x∗, ·)〉. The results in Schuster et al. (2020a) show that the empirical estimate of pPE may be constructed from m samples of the reference density {yi}i∈[m] and n training samples {xi}i∈[n] and {zi}i∈[n] as
p̂PE = (Ĉρ + α ′I)−1ĈXZ(ĈZZ + αI) −1µ̂Z = m∑ i=1 βil(yi, ·) (12)
where β = m−2(LY + α ′I)−2LY X(KZ + αI) −1Φ>µ̂Z (13)
and
LY = [l(yi, yj)]ij ∈ Rm×m, LY X = [l(yi, xj)]ij ∈ Rm×n,KZ = [k(zi, zj)]ij ∈ Rn×n (14)
In Fig. 8, we use a uniform density in the square (±4,±4) as the reference density ρ and constructed p̂PE using m = 10000 samples {yi}i∈[m] from ρ. Due to the form of the empirical kernel CDO, where the estimated density function p̂PE is a linear combination of {l(yi, ·)}i∈[m] (as in Eq. 12), the approximation can be inaccurate if reference samples are relatively sparse and the densities are ‘sharp’. In those cases, to obtain a better density estimate, we may either increase the number of reference samples used to construct the empirical CDO (which can be computationally difficult due to the need to compute (LY + α
′I)−2), or, with some prior knowledge to the true distribution, choose a reference density which is localized around the ground truth density.
Therefore, to show a more faithful density estimate of the transferred distribution for visualization purpose, we use a composite of the uniform density and the true density with weight 4 : 1 as the reference density ρ in Fig. 2. In this case, approximately 20% of the reference samples concentrates around the high-density areas of the true density, which helps to form a better basis for p̂PE . Note that the choice of ρ does not affect the transferred density embedding P̂E µ̂Z since it is independent of pX and pZ . After this modification, the reconstructed density more accurately reflects the true density compared to GMM and GLOW, indicating the transferred distribution by kPF in RKHS matches better to the ground truth distribution. This is also reflected in the generated samples in Fig. 9, where samples generated by the kPF operator are clearly more aligned with the ground truth distribution.

heading: C Effect of γ on Sample Quality
text: In the sampling stage, our proposed method finds the approximate preimage of the transferred kernel embeddings by taking the weighted Fréchet mean of the top γ neighbors among the training samples. The choice of γ therefore influences the quality of generation.
From Figure 10, we can observe that, in general, FID worsens as γ increases. This observation aligns with our intuition of preserving only the local similarities represented by the kernel, and similar ideas have been previously used in the literature (Hastie et al., 2001; Kwok & Tsang, 2004). However, significantly decreasing γ leads to the undesirable result where the generator merely generates the training samples (in the extreme case where γ = 1, generated samples will just be reconstructions of training samples). Therefore, in our experiments, we choose γ = 5 to achieve a balance between generation quality and the distance to training samples.

heading: D Weighted Fréchet Mean on the Hypersphere
text: While the weighted Fréchet Mean in Euclidean space can be computed in closed-form as a weighted arithmetic mean (as in Eq. 10), on the hypersphere there is no known closed-form solution. Thus, we adopt the iterative algorithm in (Chakraborty & Vemuri, 2015) for an approximate solution given data points X = {x1 . . .xγ} and weight vector s:
M1 = x1
Mi+1 := cos(‖si+1v‖)Mi + sin(‖si+1v‖) v
‖v‖
where, v = θsin(θ) (Xi+1 −Mi cos(θ)) , θ = arccos(M t iXi+1). This algorithm iterates through the data points once, yielding a complexity of only O(γd), where d is the dimension of X . Under the prescribed iteration, Mn converges asymptotically to the true Fréchet mean for finite data points. We refer the readers to (Chakraborty & Vemuri, 2015) for further details.

heading: E Fast approximation of Moore-Penrose inverse
text: When computing the inverted kernel matrix Kinv in Algo. 1, conventional approaches typically performs SVD or Cholesky decomposition. Both procedures are hard to parallelize, and therefore, can be slow when K is large. Alternatively, we can utilize an iterative procedure proposed in Razavi et al. (2014) to approximate the Moore-Penrose inverse.
Z1 = K/(‖K‖1‖K‖∞) (15) Zi+1 := Zi(13I −KZi(15I −KZi(7I −KZi))) (16)
where
‖K‖1 = max j n∑ i=0 Kij , ‖K‖∞ = max i n∑ j=0 Kij (17)
Since this iterative procedure mostly involves matrix multiplications, it can be efficiently parallelized and implemented on GPU. The same procedure has also seen success in approximating large self-attention matrices in language modeling (Xiong et al., 2021). For the NVAE experiment, we run this iteration for 10 steps and use Kinv = Z10.

heading: F Nystrom Approximation of kPF
text: Due to the need to store and compute a kernel matrix inverse (K + λnI)−1 or K†, the memory and computational complexity of kPF is at least O(n2) and O(n3), respectively. The sup-quadratic complexity hinders the use of kPF on extremely large datasets. In our experiments, we already adopted a simple subsampling strategy which randomly select 10k training samples from each dataset (≤ 50k samples) to fit our hardware configuration which works well. But for larger datasets with potentially more modes, a larger set of subsamples must be considered, and in those cases kPF may not be suitable for commodity/affordable hardware. In order to overcome this problem, we can combine kPF with conventional kernel approximation methods such as the Nyström method (Williams & Seeger, 2001).
Let (X?,Z?) be a size v subset of the training set (which we refer to as the landmark points) and (Ψ?,Φ?) be their corresponding kernel feature maps. The weighting coefficients s for each prior sample z∗ ∼ Z derived in Alg. 1 can be approximated by
s = L(K + λnI)†Φ>k(z?, ·) ≈ LΨ?W † Ψ? L>Ψ?(KΦ?W † Φ? K>Φ? + λnI) †Φ>k(z?, ·)
= LΨ?W † Ψ? L>Ψ?(λn) −1(I −K>Φ?(λnW † Φ? +KΦ?K > Φ?) †KΦ?)Φ >k(z?, ·) (18)
where LΨ? = Ψ >Ψ? ∈ Rn×v, WΨ? = Ψ>? Ψ? ∈ Rv×v, KΦ? = Φ>Φ? ∈ Rn×v, WΦ? = Φ>? Φ? ∈ Rv×v, and the last identity is due to applying the Woodbury formula on (KΦ?W † Φ? K>Φ?+λnI)
†. Assuming v n, the memory complexity is reduced toO(nv) and the computation complexity to O(nv2).
We empirically evaluated the Nyström-approximated kPF on the CelebA experiment and present the result in Tab. 4. It can be observed that when v is sufficiently large, the performance of Nyström approximated kPFs is as good as the ones using the full kernel matrices.

heading: G Does kPF Memorize Training Data?
text: Since in kPF, samples are generated by linearly interpolating between training samples, it is natural to wonder whether it ‘fools’ the metrics by simply replicating training samples. For comparison, we consider an alternative scheme that generates data through direct manipulation of the training data, namely Kernel Density Estimation (KDE).
We fit KDEs by varying noise levels σ and compare their FIDs and nearest samples in the latent space to kPF in Fig. 11. We observe that, although KDE can reach very low FIDs when σ is small, almost all new samples closely resemble some instance in the training set, which is a clear indication of memorization. In contrast, kPF can generate diverse samples that do not simply replicate the observed data.

heading: H Assessing Sample Diversity
text: Although FID is one of the most frequently used measures for assessing sample quality of generative models, certain diversity considerations, such as mode collapse, may not be conveniently deduced from it (Sajjadi et al., 2018). To enable explicit examination of generative models with respect to both accuracy (i.e., generating samples within the support of the data distribution) and diversity (i.e., covering the support of the data distribution as much as possible), Sajjadi et al. (2018) proposed an approach to evaluate generative models with generalized definitions of precision and recall between distributions. Quality of generation can then be assessed by evaluating the PRD curve, which depicts the trade-offs between accuracy (precision) and diversity (recall). We present the PRD curves in Fig. 12. The observations align with our results in Tab. 2 and kPF performs competitively in both accuracy and sample diversity.

heading: I Exploring Kernel Configurations
text: To investigate the implication of kernel choices on generation quality, we tested 25 different kernel configurations for CelebA generation (results are presented in Tab. 5). For RBF kernels used in the CelebA experiments of the main text, we use a bandwidth of σin = √
2|Z|/8 ≈ 0.71 when used as input kernel and σout ≈ 0.34, and we adopt the same notation here.
It can be seen that kernel configurations and parameters indeed has a non-trivial impact on the generation quality, with NTK-kPF being the most robust to the choice of parameters. This aligns with our previous observations and offers some support for using NTK as an input kernel despite the additional compute cost.

heading: J Experimental Details
text: In this section, we provide the detailed specifications for all of our experiments. We have also provided our code in the supplemental material.
J.1 Density Estimation on Toy Densities
We generated 10000 samples from each of the toy densities to learn the kPF operator. The input kernel k is a ReLU-activated NTK corresponding to a fully-connected network with depth L = 4 and width w = 10000 at each layer, and the output kernel l is a Gaussian kernel. Unless specified otherwise, we always uses a Gaussian kernel as the output kernel for the remainder of this appendix. The bandwidth of the output kernel was adjusted separately for density estimation and sampling for the purpose of demonstration. For comparison, we also fit/estimate a 10-component GMM and a Glow model with 50 coupling layers, where each of them were trained until convergence.
J.2 Image Generation with Computer Vision Datasets
To generate results in Tab. 2 and Tab. 3, we first trained an autoencoder for each dataset following the model setup in (Ghosh et al., 2020), which uses a modified Wasserstein autoencoder (Tolstikhin et al., 2018) architecture with batch normalization. Additionally, we applied spectral normalization on both the encoder and the decoder, following (Ghosh et al., 2020), to obtain a regularized autoencoder. The latent representations were projected onto a hypersphere before decoding to image space. We trained the models on two NVIDIA GTX 1080TI GPUs. A detailed model specification is provided below in Table 6.
We used an NTK with L = 8 and w = 10000 as the input kernel k (i.e. the embedding kernel of pZ) for NTK-kPF, and a Gaussian kernel with bandwidth σin = √
2|Z|/8 for RBF-kPF. The bandwidth for the output Gaussian kernel is selected by grid search over {2−i ∗ σdata|i ∈ [8]}, where σdata is the empirical data standard deviation, based on crossvalidation of a degree 3 polynomial kernel MMD between the sampled and the ground-truth latent points. Further, to mitigate the deterioration of performance of kernel methods in a high-dimensional setting due to the curse of dimensionality (Evangelista et al., 2006), in practice, we model Z as a space with fewer dimensions than the input space X . As a rule of the thumb, we choose Z such that |Z| = |X |/4. To generate images from kPF learned on NVAE latent space, we used the pre-trained checkpoints provided in (Vahdat & Kautz, 2020) to obtain the latent embeddings for 2000 FFHQ images. We then construct the kPF from the concatenated latent space of the lowest resolution (8× 8). During sampling, prior samples at those resolutions are replace by the kPF samples, while for other resolutions samples remain generated from inferred Gaussian distributions. The batchnorm statistics were readjusted for 500 iterations following (Vahdat & Kautz, 2020). We use rbf kernels as input and output kernels, with bandwidths σk, σl chosen by the median heuristic (∼ 100 for input and ∼ 70 for output in our experiments).
J.3 Image Generation for Brain Images
For the high-resolution brain imaging dataset, we used a custom version of ResNet (He et al., 2016) with 3D convolutions. The detailed architecture is shown in Fig. 7. Due to the large size of the data, we trained the model on 4 NVIDIA Tesla V100 GPUs.
Mandatory ADNI statement regarding data use. Data used in preparation of this article were obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found in the ADNI Acknowledgement List.
The T1 MR brain dataset we utilize consists of images from 184 subjects diagnosed with Alzheimers’s disease and 292 healthy controls/ normal subjects. Images were first coregistered
Encoder
5x5 conv, stride 4
ResBlock64 × 2 ResBlock64 × 2 ResBlock128 × 2 ResBlock256 × 2
Decoder
ResBlock256 × 2 ResBlock128 × 2 ResBlock64 × 2 ResBlock64 × 2
5x5 conv, stride 4
Table 7: Model architecture for experiments for image generation on brain imaging dataset. Subscript denotes the number of output channels. Upsampling and downsampling are performed using strided convolutions.
to a MNI template and segmented to preserve only the white matter and grey matter. Then, all images were resliced and resized to 160× 196× 160 and rescaled to the range of [−1, 1]. Voxel-based morphometry (VBM) was used to obtain the p-value map of data and generated images.

heading: K More Samples
text: In this section we present additional uncurated set of samples on MNIST, CIFAR-10, CelebA based on pre-trained SRAE and FFHQ based on NVAE. From the figures, it can be seen that kPF produces consistent and diverse samples, often better in quality than the alternatives.

Review summaries:
The authors propose a new type of generative model. The new scheme is based on a kernel transfer operator that leads to a cheap method for distribution matching. The authors rely on rigorous theory on RKHS and propose a framework for transferring a prior distribution linearly (in RKHS) to the data distribution. The authors demonstrate that the proposed approach leads to improved approximations of observed distributions. Specifically, the new approach can lead to the generation of new images form a given distribution and requires less training time compared to existing baselines. The paper is mostly well written, and the method relies on solid justifications. The authors demonstrate the usefulness of the method and compare it to several existing methods. 

The paper presents a novel (or an unusual type of) generative modeling approach that is kernel-based and non-parametric. The basic idea is using an operator that maps from the RKHS of Z to the RKHS of X, so that data generation can be done by mapping the prior distribution $p_{prior}(z)$ to the RKHS of Z, applying the operator, and project down the result in the RKHS of X to the data space. The operator is constructed as the RKHS analogy of the desired conditional distribution p(x|z) (which makes $\int p(x|z) p_{prior}(z) dz = p_{data}(x)$), so that it transforms the RHKS embedding of $p_{prior}(z)$ to that of $p_{data}(x)$. However, to guarantee a desired $p(x|z)$, the operator is chosen by assuming $x$ and $z$ are independent. The operator can be estimated by samples of the joint distribution, which amounts to samples of $p_{data}(x)$ (i.e., training dataset) and samples from $p_{prior}(z)$ under the independent assumption, and the maps between the sample space and the RKHS can also be estimated. Experiments show the utility of the method for data generation for densely supported distributions.

The paper proposes to leverage the Perron-Frobenius operator to simplify the mapping from Z -> X in latent variable generative models. Specifically, this “forward operator” can be approximated by a closed-form linear operator in the RKHS. They evaluate their method on both synthetic settings and high-dimensional image datasets such as MNIST, CIFAR-10, CelebA, and FFHQ.