Quantum-inspired probability metrics define a complete, universal space for statistical learning

TL;DR

This paper introduces quantum probability metrics (QPMs), a new class of tools for comparing probability distributions that are more complete and sensitive, especially in high-dimensional, non-compact spaces, with practical applications in generative modeling.

Contribution

The paper develops quantum probability metrics by embedding probability measures into quantum state space, extending kernel methods, and overcoming MMD's limitations on non-compact spaces.

Findings

QPMs are complete and universal for bounded, uniformly continuous functions.

QPMs outperform MMD in high-dimensional generative modeling tasks.

Eigenvalue-based computation allows practical implementation of QPMs.

Abstract

Comparing probability distributions is a core challenge across the natural, social, and computational sciences. Existing methods, such as Maximum Mean Discrepancy (MMD), struggle in high-dimensional and non-compact domains. Here we introduce quantum probability metrics (QPMs), derived by embedding probability measures in the space of quantum states: positive, unit-trace operators on a Hilbert space. This construction extends kernel-based methods and overcomes the incompleteness of MMD on non-compact spaces. Viewed as an integral probability metric (IPM), QPMs have dual functions that uniformly approximate all bounded, uniformly continuous functions on $\mathbb{R}^n$, offering enhanced sensitivity to subtle distributional differences in high dimensions. For empirical distributions, QPMs are readily calculated using eigenvalue methods, with analytic gradients suited for learning and optimization. Although computationally more intensive for large sample sizes ($O(n^3)$ vs. $O(n^2)$), QPMs can significantly improve performance as a drop-in replacement for MMD, as demonstrated in a classic generative modeling task. By combining the rich mathematical framework of quantum mechanics with classical probability theory, this approach lays the foundation for powerful tools to analyze and manipulate probability measures.

Tables

Purpose	Variable	Conversion formulas
Originally defined on $ℝ$,	$x$	$x=\left(z+\overline{w}\right)/2$
but extended to $ℂ$	$y$	$y=\left(z-\overline{w}\right)/\left(2i\right)$
Originally defined on $ℂ$	$z$	$z=x+iy$
	$w$	$w=\overline{x}+i\overline{y}$

Equations

$$T(\mu )={\mathbb{E}}_{x\sim \mu }[\varphi (x)]\equiv {\int }_X\varphi (x)d\mu (x)$$

$$d(\mu,\nu)=\sup_{\begin{subarray}{c}\lambda\in E^{*}\\ \|\lambda\|\leq 1\end{subarray}}\big{|}\lambda(T(\mu))-\lambda(T(\nu))\big{|}$$

$$d(\mu,\nu)=\sup_{\begin{subarray}{c}\lambda\in E^{*}\\ \|\lambda\|\leq 1\end{subarray}}\big{|}\,\mathbb{E}_{x\sim\mu}[\lambda(\phi(x)]-\mathbb{E}_{x\sim\nu}[\lambda(\phi(x)]\,\big{|}$$

$$T(\mu )={\int }_{\mathbb{C}}\varphi (z)d\mu (z)={\int }_{\mathbb{C}}{\hat{\rho }}_{z}d\mu (z)={\int }_{\mathbb{C}}|z⟩⟨z|d\mu (z)\equiv \hat{\mu }$$

$$\Vert T({\mu }_n)-T({\delta }_{x_1})\Vert$$

$$=\frac{1}{n}\Vert T({\delta }_{x_n})-T({\delta }_{x_1})\Vert$$

$$\ge \frac{1}{n}\Vert T({\delta }_{x_n})\Vert -\frac{1}{n}\Vert T({\delta }_{x_1})\Vert$$

$$>1-\frac{1}{n}\Vert T({\delta }_{x_1})\Vert \underset{n\to \infty }{\to }1$$

$$T(m)[\lambda ]={\int }_X\lambda (T({\delta }_x))dm(x)$$

$$\displaystyle\big{|}T(m)[\lambda]\big{|}$$

$$\le \Vert \lambda \Vert {\int }_X\Vert T({\delta }_x)\Vert d|m|(x)\le \Vert \lambda \Vert \Vert m\Vert \underset{x\in X}{\sup }\Vert T({\delta }_x)\Vert$$

$$\underset{\alpha }{\lim }T(m_{\alpha })[\lambda ]$$

$$={\int }_X\lambda (T({\delta }_x))dm(x)=T(m)[\lambda ]$$

$$\displaystyle\lim_{\alpha}\,\big{<}T(\mu_{n_{\alpha}}),f\big{>}$$

$$\displaystyle=\lim_{\alpha}\,\int_{X}\big{<}K_{x},f\big{>}\,d\mu_{n_{\alpha}}(x)$$

$$=\underset{\alpha }{\lim }{\int }_{X}f(x)d{\mu }_{n_{\alpha }}(x)$$

$$\left|{\int }_{C}f(x)d{\mu }_n(x)\right|\le \frac{NM}{2n+1}$$

$$|\mathbf{z}⟩=e^{-|\mathbf{z}{|}^2/2}\sum _{\mathbf{n}=0}^{\infty }(\mathbf{n}!{)}^{-1/2}{\mathbf{z}}^{\mathbf{n}}|\mathbf{n}⟩$$

$$⟨\mathbf{w}|\mathbf{z}⟩=e^{-\frac{1}{2}(|\mathbf{w}{|}^2+|\mathbf{z}{|}^2-2\overline{\mathbf{w}}\mathbf{z})}$$

$$\psi (\mathbf{z})\equiv ⟨\mathbf{z}|\psi ⟩=e^{-|\mathbf{z}{|}^2/2}\sum _{\mathbf{n}}(\mathbf{n}!{)}^{-1/2}⟨\mathbf{n}|\psi ⟩{\overline{\mathbf{z}}}^{\mathbf{n}}$$

$$\int e^{-|z{|}^2}{z}^n{\overline{z}}^m{d}^{2}z=\pi n!{\delta }_{n,m}$$

$$⟨\varphi |\psi ⟩=\frac{1}{{\pi }^n}\int \overline{\varphi (\mathbf{z})}\psi (\mathbf{z})d^{2n}\mathbf{z}=\frac{1}{{\pi }^n}\int ⟨\varphi |\mathbf{z}⟩⟨\mathbf{z}|\psi ⟩d^{2n}\mathbf{z}$$

$$\hat{I}\simeq \frac{1}{{\pi }^n}\int |\mathbf{z}⟩⟨\mathbf{z}|d^{2n}\mathbf{z}$$

$$K(\mathbf{w},\mathbf{z})=|k(\mathbf{w},\mathbf{z}){|}^2=\overline{k(\mathbf{w},\mathbf{z})}k(\mathbf{w},\mathbf{z})=⟨\mathbf{z}|\mathbf{w}⟩⟨\mathbf{w}|\mathbf{z}⟩$$

$$({\hat{\rho }}_{\mathbf{w}},{\hat{\rho }}_{\mathbf{z}})=\mathrm{tr}({\hat{\rho }}_{\mathbf{w}}{\hat{\rho }}_{\mathbf{z}})=\mathrm{tr}(|\mathbf{w}⟩⟨\mathbf{w}|\mathbf{z}⟩⟨\mathbf{z}|)=⟨\mathbf{z}|\mathbf{w}⟩⟨\mathbf{w}|\mathbf{z}⟩$$

$$⟨K_{\mathbf{w}}(\cdot ),K_{\mathbf{z}}(\cdot )⟩=K(\mathbf{w},\mathbf{z})=⟨\mathbf{z}|\mathbf{w}⟩⟨\mathbf{w}|\mathbf{z}⟩.$$

$$f(\cdot )=\sum _i{b}_i{K}_{{\mathbf{z}}_i}(\cdot )=\sum _i{b}_{i}K(\cdot ,{\mathbf{z}}_i)$$

$$⟨f(\cdot ),g(\cdot )⟩$$

$$=\sum _{i,j}\overline{b_i}{c}_j⟨K_{{\mathbf{z}}_i}(\cdot ),K_{{\mathbf{z}}_j}(\cdot )⟩=\sum _{i,j}\overline{b_i}{c}_{j}K({\mathbf{z}}_i,{\mathbf{z}}_j).$$

$$f(\mathbf{w})$$

$$=\sum _i{b}_i⟨K_{\mathbf{w}}(\cdot ),K_{{\mathbf{z}}_i}(\cdot )⟩=\sum _i{b}_{i}K(\mathbf{w},{\mathbf{z}}_i)=f(\mathbf{w})$$

$$f(\cdot )=\sum _i{b}_i{K}_{{\mathbf{z}}_i}(\cdot )\to \hat{f}=\sum _i{b}_i{\hat{\rho }}_{{\mathbf{z}}_{\mathbf{i}}}$$

$$f(\mathbf{w})$$

$$=\sum _i{b}_i⟨K_{\mathbf{w}}(\cdot ),K_{{\mathbf{z}}_i}(\cdot )⟩$$

$$=\sum _i{b}_i\mathrm{tr}({\hat{\rho }}_{\mathbf{w}}{\hat{\rho }}_{{\mathbf{z}}_i})=\mathrm{tr}\left({\hat{\rho }}_{\mathbf{w}}\sum _i{b}_i{\hat{\rho }}_{{\mathbf{z}}_i}\right)$$

$$=\mathrm{tr}({\hat{\rho }}_{\mathbf{w}}\hat{f})=\mathrm{tr}(|\mathbf{w}⟩⟨\mathbf{w}|\hat{f})=⟨\mathbf{w}|\hat{f}|\mathbf{w}⟩.$$

Full text

Given two probability measures on a set $X$, how can we measure their similarity? On this space $\mathcal{P}(X)$, one often chooses the weak topology, which is the coarsest topology ensuring convergence of expected values for bounded continuous functions $f:X\to\mathbb{R}$ [1]. If there is a metric $d(x_{1},x_{2})$ defining the distance between points in $X$, we can construct a metric on $\mathcal{P}_{\mathrm{w}}(X)$ that metrizes the weak topology (indicated by the subscript $\mathrm{w}$; we assume the Borel $\sigma$-algebra throughout). Many such metrics are integral probability metrics, or IPMs. Given a family of functions $\mathcal{F}$, the IPM between two probability measures $\mu$ and $\nu$ is the supremum of the difference in expected values for all $f\in\mathcal{F}$: d(\mu,\nu)=\sup_{f\in\mathcal{F}}\big{|}\mathbb{E}_{x\sim\mu}[f(x)]-\mathbb{E}_{x\sim\nu}[f(x)]\big{|}.

The functions in $\mathcal{F}$ are witness functions for the IPM: they detect differences between probability measures. For example, when $\mathcal{F}$ is the set of functions with Lipschitz constant $L\leq 1$, the resulting metric is the 1-Wasserstein metric, or earth mover’s metric. As this definition is not very useful for practical calculations, IPMs are often constructed by embedding $\mathcal{P}_{\mathrm{w}}(X)$ into a normed linear space $E$, such as a Banach or Hilbert space [2, 3]. Given a continuous, injective function $\phi:X\to E$, define the affine barycenter map $T:\mathcal{P}_{\mathrm{w}}(X)\to E$ as the expected value of $\phi$:

[TABLE]

For appropriate $\phi$, this integral is well-defined (as a Bochner or Pettis integral) and provides a topological embedding that respects the convex (affine) structure of $\mathcal{P}_{\mathrm{w}}(X)$. The norm on $E$ induces a metric, $d(\mu,\nu)=\left\|T(\mu)-T(\nu)\right\|$, equivalently defined over the dual space $E^{*}$ of continuous linear functionals $\lambda$ that act on $E$:

[TABLE]

The barycenter map $T$ plus linearity of the expected value gives:

[TABLE]

This is an IPM: $\lambda\in E^{*}$ defines a dual function $\lambda(\phi(x))$ on $X$; the witness functions have $\|\lambda\|\leq 1$. Typically, calculating the norm of embedded probability measures is simpler than finding the supremum over an unwieldy class of functions. For instance, if we embed $\mathcal{P}_{\mathrm{w}}(X)$ into a Hilbert space $\mathcal{H}$ with inner product $\langle\cdot,\cdot\rangle$, then $\phi$ is the feature map for a reproducing kernel Hilbert space (RKHS) with kernel $k(x,y)=\langle\phi(x),\phi(y)\rangle$ [2]. With an appropriate kernel this embedding defines a metric on $\mathcal{P}_{\mathrm{w}}(X)$ known as the Maximum Mean Discrepancy (MMD). The dual functions are precisely the functions in the RKHS, with IPM witness functions contained in the unit ball.

Limitations of MMD

The limited space of RKHS functions imposes limits on MMD. For the widely used Gaussian, Laplacian, or inverse-multiquadric kernels, all RKHS functions vanish at infinity; even constants are excluded. In addition, when the space $X$ is not compact, the MMD metric is incomplete. Consider generating a sequence of probability measures $\mu_{n}$ by iterative sampling or modeling. Without an a priori limit, we can only check if the sequence is Cauchy: do the terms get closer to one another as $n$ grows? In a complete metric space, every Cauchy sequence converges to a valid probability measure. But in an incomplete space, a seemingly convergent sequence may diverge to something invalid. This problem arises with embeddings into reflexive Banach spaces, including all Hilbert spaces and most reproducing kernel Banach spaces [3]. We prove in SI Appendix, A.1:

Theorem 1**.**

For noncompact $X$, an affine embedding of $\mathcal{P}_{\mathrm{w}}(X)$ into a reflexive Banach space $E$ is never closed (in the norm topology). Thus, the induced metric is incomplete.

Quantum coherent states

A complete metric needs a non-reflexive Banach space, such as the quantum state space. Elementary quantum states $|\psi\rangle$ are vectors in a Hilbert space $\mathcal{H}$, but a proper treatment uses the full state space $\mathcal{S}$ of positive, unit-trace operators on $\mathcal{H}$. The extreme points of this convex space are rank-1 projection operators $\hat{\rho}=|\psi\rangle\langle\psi|$ called pure states. Their convex combinations form mixed states such as $\hat{\rho}=\sum_{i}p_{i}|\psi_{i}\rangle\langle\psi_{i}|$ for a set of states $|\psi_{i}\rangle$ with probabilities $p_{i}$. This space $\mathcal{S}$ is contained in the non-reflexive Banach space of trace-class operators $\mathcal{B}_{1}(\mathcal{H})$, with norm $\|\hat{A}\|_{1}=\mathrm{tr}(|\hat{A}|)$, whose dual is the space of bounded operators $\mathcal{B}(\mathcal{H})$ with the operator norm. (We use Dirac notation and “hats” for operators; see SI Appendix, B, for details.)

A special class of quantum states are the coherent states of the quantum harmonic oscillator [4]. (In the math literature, this construction is known as the Fock space over $\mathbb{C}$.) Given a pair of complex numbers $w,z$, the corresponding coherent states satisfy $\langle w|z\rangle=\exp[-\frac{1}{2}(|w|^{2}+|z|^{2}-2w^{*}z)]$. These states form a closed embedding of $\mathbb{C}$ in $\mathcal{H}$ with the map $\phi(z)=|z\rangle$, and define an RKHS with kernel $k(w,z)=\langle w|z\rangle$.

Although the coherent state kernel is not typically used for RKHS embeddings, its norm-square is precisely the Gaussian kernel on $\mathbb{C}$: $|k(w,z)|^{2}=\exp(-|w-z|^{2})$. This kernel arises naturally if we embed $z\in\mathbb{C}$ as a pure quantum state in $\mathcal{S}$, with embedding $\phi(z)=|z\rangle\langle z|\equiv\hat{\rho}_{z}$. The Hilbert-Schmidt inner product between two embedded states is $\mathrm{tr}(\hat{\rho}_{\mathrm{w}}\hat{\rho}_{z})=\exp(-|w-z|^{2})$. Extending to $\mathbb{C}^{n}$ gives (SI Appendix, A.2):

Theorem 2**.**

Coherent states of an $n$-mode quantum harmonic oscillator, with the Hilbert-Schmidt inner product, are isometric to kernel functions in the Gaussian RKHS over $\mathbb{C}^{n}$.

Quantum probability metrics

For pure quantum states, the “natural” topology is induced equivalently by the trace norm and the Hilbert-Schmidt norm. For mixed quantum states, however, we will find the trace distance to be advantageous. It is conventionally normalized as $d(\hat{\rho},\hat{\sigma})=\frac{1}{2}\mathrm{tr}(|\hat{\rho}-\hat{\sigma}|)$ to take values in $[0,1]$. Given eigenvalues $\lambda_{i}$ of the compact self-adjoint operator $(\hat{\rho}-\hat{\sigma})$, the trace distance is $\frac{1}{2}\sum_{i}|\lambda_{i}|$.

We first build a Quantum Probability Metric (QPM) for measures on $\mathbb{C}$. Embed $z\in\mathbb{C}$ using the coherent-state map $\phi(z)=|z\rangle\langle z|$. The barycenter map $T:\mathcal{P}_{\mathrm{w}}(\mathbb{C})\to\mathcal{S}$ gives:

[TABLE]

where the “hat” signals the operator embedding of the measure. In physics this is known as the $P$-representation for the measure $d\mu(z)=P(z)d^{2}z$ on $\mathbb{C}$; this is not a contour integral, but rather over the 2D Lebesgue measure $d\Re z\;d\Im z$. From Theorem 2, this is precisely the Gaussian kernel embedding. The MMD is the Hilbert-Schmidt distance $d(\mu,\nu)=\mathrm{tr}(|\hat{\mu}-\hat{\nu}|^{2})^{1/2}$, which metrizes $\mathcal{P}_{\mathrm{w}}(\mathbb{C})$, albeit with an incomplete metric. The Gaussian kernel is characteristic, which means that the barycenter map is injective [2]. This construction, which we call the Fock embedding of $\mathcal{P}_{\mathrm{w}}(\mathbb{C})$, extends to probability measures on $\mathbb{C}^{n}$, and even to measures on a separable infinite-dimensional Hilbert space [5].

The QPM between two measures is the trace distance $d(\mu,\nu)=\frac{1}{2}\mathrm{tr}(|\hat{\mu}-\hat{\nu}|)$, which completely metrizes $\mathcal{P}_{\mathrm{w}}(\mathbb{C})$ (i.e. the metric is complete). Completeness relies on some remarkable properties of the barycenter map on a closed set of pure states $\mathcal{C}\subset\mathcal{S}$: it is a continuous, closed surjection from $\mathcal{P}_{\mathrm{w}}(\mathcal{C})$ onto the norm-closed convex hull of $\mathcal{C}$ [6]; if it is injective, then it is a homeomorphism onto its image. Thus, given a characteristic kernel, the barycenter map yields a closed topological embedding of $\mathcal{P}_{\mathrm{w}}(\mathcal{C})$ into $\mathcal{S}$, and the trace distance completely metrizes the topology on $\mathcal{P}_{\mathrm{w}}(\mathcal{C})$.

We can now construct a general QPM. Choose $\phi:X\to\mathcal{S}$ to map points $x\in X$ to pure states $|x\rangle\langle x|\equiv\hat{\rho}_{x}$ on a Hilbert space. We require that $\phi$ is a closed embedding, so its image $\mathcal{C}$ is closed and homeomorphic to $X$, and also that its kernel $k(x,y)=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$ is characteristic. Embedding $\mathcal{P}_{\mathrm{w}}(X)$ into $\mathcal{S}$ using $\phi$ and the (injective) barycenter map, the QPM is the trace distance between the embedded measures. Adding some technical details in SI Appendix, A.3, we have:

Theorem 3**.**

Every QPM completely metrizes $\mathcal{P}_{\mathrm{w}}(X)$.

All QPMs embed into the same state space $\mathcal{S}$, which is a universal embedding space: using [7, 5] in SI Appendix, A.4, if $X$ is any separable, completely metrizable space—a Polish space—there is a map $\phi$ that meets the requirements of Theorem 3, so there is a QPM embedding into $\mathcal{S}$:

Theorem 4**.**

Every Polish space has a QPM.

Dual functions for a QPM

Since every QPM is also an IPM, we can identify its dual functions using the duality between real Banach spaces of self-adjoint operators: $\mathcal{B}_{1,\mathrm{sa}}(\mathcal{H})$ (trace class) and $\mathcal{B}_{\mathrm{sa}}(\mathcal{H})$ (bounded operators). In quantum mechanics, an operator $\hat{F}\in\mathcal{B}_{\mathrm{sa}}(\mathcal{H})$ is called an observable; it is a linear functional on $\hat{\rho}\in\mathcal{B}_{1,\mathrm{sa}}(\mathcal{H})$ via the trace map, $\mathrm{tr}(\hat{\rho}\hat{F})$. Given a QPM that maps $x\in X$ to a pure state $\hat{\rho}_{x}=|x\rangle\langle x|$, every observable defines a real function $F(x)=\mathrm{tr}(\hat{\rho}_{x}\hat{F})=\langle x|\hat{F}|x\rangle$. These are dual functions for the QPM; since $\hat{F}$ is bounded they are bounded and continuous. Embedding $\mu$ as a state $\hat{\mu}\in\mathcal{S}$, the (quantum) expectation value $\mathrm{tr}(\hat{\mu}\hat{F})$ equals the (classical) expected value $\mathbb{E}_{x\sim\mu}[F(x)]$. The IPM witness functions for a QPM have operator norm $\|\hat{F}\|\leq 1$. Even on compact domains, QPMs have a larger class of witness functions than MMD, which requires Hilbert-Schmidt norm $\|\hat{F}\|_{2}\leq 1$.

On $\mathbb{R}^{n}$, the Fock embedding QPM has a rich class of dual functions that are not required to vanish at infinity, unlike those of MMD. In SI Appendix, A.5, using the Schur test [8], we show that all bounded band-limited functions on $\mathbb{R}^{n}$ are included as dual functions; combined with [9], we prove:

Theorem 5**.**

Dual functions for the Fock (coherent state) QPM are dense in $BUC(\mathbb{R}^{n})$: they can uniformly approximate any bounded, uniformly continuous function on $\mathbb{R}^{n}$.

Replacing MMD with a QPM

In most cases, MMD is calculated for atomic measures (convex combinations of Dirac measures). Although the underlying RKHS is usually infinite-dimensional, the MMD between two such measures can be calculated using kernel evaluations $k(x,y)$ [2]. If their joint support contains $n$ points, MMD requires the Gram matrix for every pair of embedded vectors, with complexity $O(n^{2})$.

The analogous QPM calculation needs eigenvalues of the rank-$n$ operator $(\hat{\mu}-\hat{\nu})$, which has the form $\sum_{i}c_{i}|x_{i}\rangle\langle x_{i}|$. Its eigenvalues are equal to those of the Hermitian matrix $M=H^{\dagger}CH$, where we factor the Gram matrix as $G=HH^{\dagger}$ and define $C=\mathrm{diag}(c_{i})$. Although eigenvalue computation is $O(n^{3})$, this is a simple “drop-in” replacement for MMD, perhaps with a modified kernel (SI Appendix, C). Anyone who uses MMD should try this swap; the improved metric performance is likely worth the computational cost. Eigenvalue calculations are numerically stable with analytic gradients for optimization. By comparison, the widely-used $p$-Wasserstein metrics are also $O(n^{3})$ with the Hungarian algorithm, but they are difficult to differentiate, need finite $p$th moments, and have the opposite problem of MMD: sequences that diverge in the Wasserstein metric might converge in $\mathcal{P}_{\mathrm{w}}(X)$.

As a proof of concept, we swap MMD with QPM in a classic model, the Generative Moment Matching Network [10]; see SI Appendix, D, for details. These networks are trained to minimize the MMD (in the image space) between batches of generated images and batches drawn from the dataset. This is a great idea that gives poor results, as shown in Figure 1 where we reproduce the original work on the MNIST dataset [11] (784 dimensions) and extend it to the CelebA-64 dataset [12] (12,288 dimensions). The CelebA result illustrates why this approach fails: MMD believes that it has met its objective, as shown by an MMD two-sample test [2] that finds no statistically significant difference between the generated images and the data ($p\approx 0.23$). The MMD witness functions, which vanish at infinity, simply fail in this high-dimensional case; a different architecture or training method would be unlikely to help. Replacing MMD with a QPM, keeping everything else constant, yields much better images. And QPM knows that it has not met its objective: its two-sample test gives $p<10^{-3}$. With its richer class of witness functions on this compact, high-dimensional space, QPM easily distinguishes the generated images from the data. This experiment highlights QPM as a target for generative modeling and a discerning metric for probability distributions.

Future work

Using QPMs, we can study and manipulate probability measures with powerful tools from quantum mechanics such as unitary transformations, group representations, tensor products, entropy, information theory, and more. In addition, this construction offers an interesting independent motivation for the mathematical structure of quantum mechanics.

Acknowledgments

The AI model “GPT-5 Thinking” from OpenAI [13] assisted with checking proofs in the SI Appendix.

SI Appendix, Materials and Methods

A Proofs

A.1 Theorem 1

For noncompact $X$, an affine embedding of $\mathcal{P}_{\mathrm{w}}(X)$ into a reflexive Banach space $E$ is never closed (in the norm topology). Thus, the induced metric is incomplete.

Note that if $\mathcal{P}_{\mathrm{w}}(X)$ can be embedded in a Banach space, then $\mathcal{P}_{\mathrm{w}}(X)$ is metrizable, so $X$ is metrizable via the canonical, continuous embedding that maps $x\in X$ to its point mass $\delta_{x}$. We begin with two preparatory lemmas.

Lemma 1**.**

$\mathcal{P}_{\mathrm{w}}(X)$* is not closed for noncompact, metrizable $X$.*

Proof.

We review results from [1], which proves this lemma slightly differently. For noncompact $X$, the dual of $C_{b}(X)$ (the Banach space of bounded continuous real functions) is the space of signed, finite, regular, finitely additive measures on $X$; we denote these generalized measures as $\mathcal{M}(X)$. The probability measures $\mathcal{P}(X)$ are a proper subset of $\mathcal{M}(X)$: they are positive, countably additive, and normalized to one. The “weak topology” on $\mathcal{P}_{\mathrm{w}}(X)$ is the subspace topology inherited from the weak-* topology on $\mathcal{M}(X)$.

Using sequential noncompactness, the metrizable space $X$ has a countable closed subset $S$ with no accumulation points; we index it as $S=\{x_{n}:n\in\mathbb{Z}\}$. Since $C_{b}(S)\cong\ell^{\infty}(\mathbb{Z})$, its dual can be identified with the space of generalized measures $\mathcal{M}(S)\cong\mathcal{M}(\mathbb{Z})$, or the ba space over $\mathbb{Z}$. The probability measures $\mathcal{P}_{\mathrm{w}}(S)$ correspond to the unit simplex in $\ell^{1}(\mathbb{Z})$, with the weak-* topology inherited from $\mathcal{M}(S)$. Note that $\mathcal{M}(S)$ is a weak--closed subset of $\mathcal{M}(X)$: define $I_{S}=\{f\in C_{b}(X):f|_{S}\equiv 0\}$, then $\mathcal{M}(S)=\bigcap_{f\in I_{S}}\{\mu\in\mathcal{M}(X):\mu(f)=0\}$, an intersection of weak--closed sets.

Define the left shift operator $L$ on $\mathcal{M}(S)$: for any index set $I\subseteq\mathbb{Z}$, $L\mu(\{x_{k}:k\in I\})=\mu(\{x_{k}:k+1\in I\})$. Build a shift-invariant measure as follows: For each $n\in\mathbb{N}$, define $S_{n}=\{x_{k}:-n\leq k\leq n\}$ and let $\mu_{n}$ be the uniform probability measure on $S_{n}$. By Banach-Alaoglu, some subnet $\mu_{n_{\alpha}}$ converges weak-* to a measure $\mu\in\mathcal{M}(S)$. A short calculation using the Følner property of $S_{n}$ shows that $\mu$ is weakly shift-invariant (i.e. $L\mu(f)=\mu(f)$ for $f\in C_{b}(S)$). By applying the shift $m$ times we see that the singleton measure $\mu(\{x_{k}\})=\mu(\{x_{k+m}\})$ for any $k$, which implies $\mu(\{x\})$ is constant for all $x\in S$.

Yet a probability measure on a countable set cannot be shift-invariant. Countable additivity forces either $\mu(S)=0$, if $\mu(\{x\})=0$, or $\mu(S)=\infty$, if $\mu(\{x\})>0$, contradicting $\mu(S)=1$. So $\mu$ is not countably additive (just like a Banach limit) and $\mu\notin\mathcal{P}_{\mathrm{w}}(S)$. This shows $\mathcal{P}_{\mathrm{w}}(S)$ is not closed. Now $\mathcal{P}_{\mathrm{w}}(S)$ = $\mathcal{P}_{\mathrm{w}}(X)\cap\mathcal{M}(S)$. Since $\mathcal{M}(S)$ is weak-* closed, but $\mathcal{P}_{\mathrm{w}}(S)$ is not, $\mathcal{P}_{\mathrm{w}}(X)$ is not closed. ∎

Lemma 2**.**

Any continuous affine map $T:\mathcal{P}_{\mathrm{w}}(X)\to E$ (norm topology on the Banach space $E$) has a canonical linear extension to $T:\mathcal{M}(X)\to E^{**}$ that is continuous from the weak- topology on $\mathcal{M}(X)$ to the weak-* topology on the double dual $E^{**}$.*

Proof.

The central idea is that $T$ is determined by its action on the point measures $\delta_{x}$, whose convex combinations are dense in $\mathcal{P}_{\mathrm{w}}(X)$ (Theorem 4.3 in [1]). Their images $T(\delta_{x})$ are bounded in $E$: suppose, for contradiction, that they are unbounded. Then there is a sequence $x_{n}$ with $\|T(\delta_{x_{n}})\|>n,\forall n\in\mathbb{N}$. Construct a sequence of probability measures $\mu_{n}=\tfrac{1}{n}\delta_{x_{n}}+(1-\tfrac{1}{n})\delta_{x_{1}}$, which converges in $\mathcal{P}_{\mathrm{w}}(X)$ to $\delta_{x_{1}}$. Yet the image of this sequence under the affine map $T$ does not converge to $T(\delta_{x_{1}})$:

[TABLE]

which contradicts continuity of $T$. In addition, $T(\delta_{x})$, viewed as a function of $x$, is continuous from $X$ to $E$, as it is the composition of the continuous embedding $x\to\delta_{x}$ with the continuous map $T$.

Now extend $T(m)$ for a generalized measure $m\in\mathcal{M}(X)$. Given a linear functional $\lambda\in E^{*}$, define $T(m)\in E^{**}$ by its action on $\lambda$:

[TABLE]

Note that this “integral” of the real-valued function $\lambda(T(\delta_{x}))$ over the generalized measure $m$ arises from viewing $m\in\mathcal{M}(X)$ as an element of the continuous dual of $C_{b}(X)$; it is not a Lebesgue integral since $m$ may not be countably additive. This expression defines a linear functional of $\lambda\in E^{*}$ that is also linear in $m$. We show that it is continuous in $\lambda$ (making it an element of $E^{**}$), and that it is continuous in $m$ (weak-* to weak-* topology).

Continuity in $\lambda$ follows from boundedness via the following calculation. Note that $|m|$ is the absolute value $(m_{+}-m_{-})$ of the signed measure $m$, while $\|m\|$ is its total variation norm—the dual norm for $\mathcal{M}(X)\cong C_{b}(X)^{*}$:

[TABLE]

where the final bound follows from boundedness of $T(\delta_{x})\in E$ and finite total variation of $m\in\mathcal{M}(X)$. Hence $T(m)\in E^{**}$.

Continuity in $m$ follows from the observation that $\lambda(T(\delta_{x}))\in C_{b}(X)$: it is the composition of continuous maps, and it is bounded since $T(\delta_{x})$ is bounded. Thus, given a weak-* convergent net of generalized measures $m_{\alpha}\to m$, we find:

[TABLE]

by weak-* convergence in $\mathcal{M}(X)$ for $\lambda(T(\delta_{x}))\in C_{b}(X)$.

Using the affine property of $T$, linearity of $\lambda$, and linearity of the integral, this extended $T$ agrees with the original $T(\mu)$ for convex combinations of point measures, which are dense in $\mathcal{P}_{\mathrm{w}}(X)$. By weak-* continuity, therefore, it agrees on all of $\mathcal{P}_{\mathrm{w}}(X)$.

Note that this construction can also be viewed as the Dunford integral of the $E$-valued map $T(\delta_{x})$ over the generalized measure $m$; such “integrals” generally lie in the double dual $E^{**}$. ∎

We now return to the main theorem.

Proof.

(Theorem 1, SI Appendix, A.1) Assume we have an affine closed embedding map $T$ (a homeomorphism from the weak topology on $\mathcal{P}_{\mathrm{w}}(X)$ to its norm-closed image in $E$); its inverse $T^{-1}$ is also affine. Since $\mathcal{P}_{\mathrm{w}}(X)$ is convex, its image is convex. Using Lemma 1, choose a net $\mu_{\alpha}\in\mathcal{P}_{\mathrm{w}}(X)$ that is weak-* convergent to a shift-invariant generalized measure $\mu$. By Lemma 2, the image of this net $T(\mu_{\alpha})$ converges weak-* to $T(\mu)\in E^{**}$. Since $E$ is reflexive, weak-* convergence is equivalent to weak convergence, and $E^{**}\cong E$, so $T(\mu_{\alpha})$ converges weakly to $T(\mu)\in E$. Define $v_{\alpha}\equiv T(\mu_{\alpha})$ as its image in $E$, and $v\equiv T(\mu)$ as its weak limit.

Using Mazur’s Lemma for nets, construct $w_{\alpha,n}$ on the directed set defined by $(\alpha^{\prime},n^{\prime})\geq(\alpha,n)$ when $\alpha^{\prime}\geq\alpha$ and $n^{\prime}\geq n$, such that (i) each $w_{\alpha,n}$ is a finite convex combination of tails of $v_{\alpha}$, i.e. from $\{v_{\beta}:\beta\geq\alpha\}$, and (ii) $\|w_{\alpha,n}-v\|\leq 1/n$, so $w_{\alpha,n}$ converges strongly to $v=T(\mu)$ in $E$. By assumption, this net and its (strong) limit lie within the norm-closed convex image of $\mathcal{P}_{\mathrm{w}}(X)$.

The affine inverse $T^{-1}$ defines a net $\nu_{\alpha,n}\equiv T^{-1}(w_{\alpha,n})$; these are probability measures since they are convex combinations of tails of $\mu_{\alpha}$. Now $T^{-1}$ is continuous from the norm topology on $E$ to the weak topology on $\mathcal{P}_{\mathrm{w}}(X)$, so the net $\nu_{\alpha,n}$ converges in $\mathcal{P}_{\mathrm{w}}(X)$ to $\nu\equiv T^{-1}(v)$; it also converges in the ambient weak-* topology on $\mathcal{M}(X)$. Since each $\nu_{\alpha,n}$ is built from tails of $\mu_{\alpha}$, both nets must converge to the same limit in this topology. Yet the limit of $\mu_{\alpha}$ is the shift-invariant generalized measure $\mu$, which cannot be a probability measure by Lemma 1. This contradicts the assumption that the image of $\mathcal{P}_{\mathrm{w}}(X)$ in $E$ is norm-closed.

Note that reflexivity of $E$ is essential. Otherwise, the image $T(\mu)$ of the shift-invariant generalized measure $\mu$ might not lie within $E$; the extension of $T$ only ensures it lies in $E^{**}$, and it is merely a weak-* limit, not a weak limit. Moreover, Mazur’s Lemma can “upgrade” weakly-convergent nets to norm-convergent nets using convex combinations, but this does not work for weak-*-convergent nets, so the key contradiction does not hold for nonreflexive $E$. ∎

As a corollary, we note that an embedding of $\mathcal{P}_{\mathrm{w}}(X)$ that is not closed can lead to pathological features such as including the zero vector in the closure of the embedded image. Moreover, since the image is convex, including zero means that images of subnormalized positive measures with $\mu(X)<1$ are also included. This is the case for the usual Kernel Mean Embedding of $\mathcal{P}_{\mathrm{w}}(X)$ for many widely-used kernels, such as the Gaussian, Laplacian, or inverse-multiquadric kernels on $\mathbb{R}^{n}$ [2], because these kernels vanish at infinity:

Corollary 1**.**

For noncompact, locally-compact $X$, the zero vector is contained in the closure of the image of $\mathcal{P}_{\mathrm{w}}(X)$ under the Kernel Mean Embedding into any RKHS of bounded functions that vanish at infinity.

Proof.

Recall that the Kernel Mean Embedding is the affine map $T:\mathcal{P}_{\mathrm{w}}(X)\to\mathcal{H}$ defined by $T(\delta_{x})=K_{x}$ for the kernel function $K_{x}(\,\cdot\,)\equiv K(\,\cdot\,,x)\in\mathcal{H}$ [2]. By Lemma 2, this map can be extended to all generalized measures, and by the argument in Theorem 1 (SI Appendix, A.1), given any weak-* convergent net of probability measures $\mu_{\alpha}\to\mu$, the image of the limit $T(\mu)$ is contained in the norm closure of the image of $\mathcal{P}_{\mathrm{w}}(X)$.

Using the construction in Lemma 1, choose a weak-* convergent subnet $\mu_{n_{\alpha}}$ that converges to a shift-invariant generalized measure $\mu$ supported on the discrete set $S$. Choose any function $f\in\mathcal{H}$ and evaluate the inner product with $T(\mu_{n_{\alpha}})$ using the reproducing property of the kernel:

[TABLE]

We show that this expression converges to zero for any $f\in\mathcal{H}$, which implies that the limit $T(\mu)$ must vanish. Choose $\epsilon>0$. Since $f$ vanishes at infinity, there is a compact set $C\subset X$ such that $|f(x)|<\epsilon/2$ for any $x\in X\backslash C$. Now recall that $\mu_{n_{\alpha}}$ is a subnet of the sequence $\mu_{n}$, and each $\mu_{n}$ is the uniform probability measure on the set $S_{n}=\{x_{k}:-n\leq k\leq n\}$. Let $N$ be the number of points in $S$ contained in $C$; it must be finite since $S$ has no accumulation points. The probability mass inside $C$ vanishes at the rate $\mu_{n}(C)\leq N/(2n+1)$, which gives the bound:

[TABLE]

for $M=\sup|f(x)|$, which is bounded for $f\in\mathcal{H}$. Choose $n$ so this integral is less than $\epsilon/2$; we already know that the integral of $|f|\,d\mu_{n}$ over $X\backslash C$ is less than $\epsilon/2$. So the sequence $\langle T(\mu_{n}),f\rangle$ converges to zero, and thus any subnet $\langle T(\mu_{n_{\alpha}}),f\rangle$ also converges to zero. Hence the weak limit (zero vector) is contained in the weak closure, and hence also in the norm closure, of the (convex) image of $\mathcal{P}_{\mathrm{w}}(X)$ (recall Hahn-Banach implies that weak- and norm-closures coincide). This reflects a failure of tightness—in the sense of Prokhorov’s theorem—since for noncompact X probability measures can “escape to infinity,” which causes the expected value to vanish for any function that vanishes at infinity. ∎

A.2 Theorem 2

Coherent states of an $n$-mode quantum harmonic oscillator, with the Hilbert-Schmidt inner product, are isometric to kernel functions in the Gaussian RKHS over $\mathbb{C}^{n}$.

We first review the well-known Segal-Bargmann-Fock correspondence. Unlike the usual Bargmann representation, which uses holomorphic functions, we follow Segal and choose anti-holomorphic functions [4]:

Lemma 3**.**

The Hilbert space of an $n$-mode quantum harmonic oscillator is isometrically isomorphic to the RKHS of complex square-integrable functions over $\mathbb{C}^{n}$ of the form $f(\mathbf{z})=\exp(-|\mathbf{z}|^{2}/2)\overline{F(\mathbf{z})}$, where $F$ is an entire function of $\mathbf{z}$.

Proof.

The Hamiltonian of the $n$-mode quantum harmonic oscillator has a complete set of orthonormal eigenstates known as Fock states indexed by $n$-dimensional vectors of nonnegative integers; we denote these states as $|\mathbf{n}\rangle\equiv|n_{1},n_{2},\ldots,n_{n}\rangle$. Each complex vector $\mathbf{z}\in\mathbb{C}^{n}$ defines a normalized coherent state [4]:

[TABLE]

where $\mathbf{n}!=n_{1}!n_{2}!\ldots n_{n}!$ and $\mathbf{z}^{\mathbf{n}}=z_{1}^{n_{1}}z_{2}^{n_{2}}\ldots z_{n}^{n_{n}}$ are standard multi-index notation, and we sum over all $\mathbf{n}\in\mathbb{N}_{0}^{n}$. A calculation using the orthonormality of Fock states gives the inner product of two coherent states (and confirms $\langle\mathbf{z}|\mathbf{z}\rangle=1$):

[TABLE]

where $\overline{\mathbf{w}}\mathbf{z}$ is the inner product in $\mathbb{C}^{n}$. We denote the Hilbert space of the $n$-mode harmonic oscillator as $\mathcal{H}_{F}$; the coherent states are a total set in this space. This construction of the Fock space over $\mathbb{C}^{n}$ can be performed over any separable Hilbert space [4].

Every vector $|\psi\rangle\in\mathcal{H}_{F}$ defines a unique function $\psi(\mathbf{z})$ on $\mathbb{C}^{n}$:

[TABLE]

where the power series converges for all $\mathbf{z}$ by boundedness of $|\langle\mathbf{n}|\psi\rangle|\leq\|\psi\|$ and rapid decrease of the prefactor $(\mathbf{n}!)^{-1/2}$. This function has the desired form: a Gaussian weight $\exp(-|\mathbf{z}|^{2}/2)$ times the complex conjugate of an entire function of $\mathbf{z}$. Uniqueness follows from uniqueness of the power series coefficients, which give the expansion of $|\psi\rangle$ in the basis of Fock states. Convergence of this power series ensures weak convergence in $\mathcal{H}_{F}$ of the power-series definition of the coherent state $|\mathbf{z}\rangle$, which in turn ensures strong convergence since these states are normalized to 1.

The following identity (readily derived in polar coordinates) shows that the monomials $z^{n}$ are orthogonal when integrated over $\mathbb{C}$ with Gaussian weight:

[TABLE]

This identity links the orthogonality of single-mode Fock states with orthogonality of the corresponding monomials. In $\mathbb{C}^{n}$ the integral factorizes, yielding $\pi^{n}\prod_{k}n_{k}!\,\delta_{\mathbf{n},\mathbf{m}}$, so one can easily show that the inner product in $\mathcal{H}_{F}$ is given by an integral over the corresponding functions:

[TABLE]

where $d^{2n}\mathbf{z}$ is the $2n$-dimensional Lebesgue measure over the real and imaginary components of $\mathbf{z}$. This inner product defines an RKHS of functions on $\mathbb{C}^{n}$. Conversely, any function of the form $\exp(-|\mathbf{z}|^{2}/2)\overline{F(\mathbf{z})}$ is uniquely determined by the power-series coefficients of the entire function $F(\mathbf{z})$. Using the coherent-state expansion, one can extract the unique Fock-basis coefficients of the corresponding vector in $\mathcal{H}_{F}$. Thus we have an isometric bijection between $\mathcal{H}_{F}$ and the desired coherent state RKHS.

The kernel for this RKHS is simply $k(\mathbf{w},\mathbf{z})=\langle\mathbf{w}|\mathbf{z}\rangle$. This kernel is strictly positive definite (s.p.d.), which means that for any distinct $\mathbf{z}_{i}$, the matrix $M_{ij}=k(\mathbf{z}_{i},\mathbf{z}_{j})$ is an s.p.d. matrix [2]. This s.p.d. Gram matrix $\langle\mathbf{z}_{i}|\mathbf{z}_{j}\rangle$ implies linear independence for any distinct, finite set of coherent states.

Although this RKHS may look like $L^{2}(\mathbb{C}^{n})$, it only contains functions of the form $\exp(-|\mathbf{z}|^{2}/2)\overline{F(\mathbf{z})}$, with $F$ entire. This is a closed subspace of $L^{2}(\mathbb{C}^{n})$ [4]: briefly, $L^{2}$ convergence of Gaussian-weighted entire functions implies uniform convergence on compact sets, which yields a Gaussian-weighted entire limit by Morera’s theorem. The expansion of the inner product in coherent states suggests the well-known resolution of the identity operator:

[TABLE]

where the notation $\simeq$ means equality in the sense of the weak operator topology, i.e., it is an equality when “sandwiched” between any pair of states $\langle\phi|\cdot|\psi\rangle$. ∎

We now return to the main theorem.

Proof.

(Theorem 2, SI Appendix, A.2) Lemma 3 connects the space $\mathcal{H}_{F}$ with the coherent state RKHS of square-integrable functions over $\mathbb{C}^{n}$. We now connect this coherent state RKHS with the Gaussian RKHS over $\mathbb{C}^{n}$.

Every RKHS is uniquely defined by its kernel (Theorem 2.5 in [2]); the Gaussian kernel is $K(\mathbf{w},\mathbf{z})=\exp(-|\mathbf{w}-\mathbf{z}|^{2})$. A short calculation using the expression for the inner product of coherent states (eq. 2 of Lemma 3) shows that the Gaussian kernel is given by the square of the coherent state kernel:

[TABLE]

As the coherent state kernel is strictly positive definite, so is the Gaussian kernel: its matrix is the Hadamard product of the individual kernel matrices, so strict positive definiteness follows by the Schur product theorem.

Define the pure coherent states $\hat{\rho}_{\mathbf{z}}=|\mathbf{z}\rangle\langle\mathbf{z}|$. These rank-one operators belong to the Hilbert-Schmidt space $\mathcal{B}_{2}(\mathcal{H}_{F})$ of operators with inner product $(\hat{A},\hat{B})\equiv\mathrm{tr}(\hat{A}^{\dagger}\hat{B})$. The inner product between two coherent states is:

[TABLE]

where the final equality is derived by taking the trace in an orthonormal basis that includes $|\mathbf{z}\rangle$ as one of its vectors.

Now recall that the Gaussian RKHS includes Gaussian kernel functions $K_{\mathbf{z}}(\,\cdot\,)=K(\,\cdot\,,\mathbf{z})$. The reproducing property means that the inner product between two kernel functions is defined as:

[TABLE]

This provides an isometric bijection between the kernel functions $K_{\mathbf{z}}(\,\cdot\,)=\exp(-|\,\cdot\,-\mathbf{z}|^{2})$ (Gaussians centered at $\mathbf{z}$) and the pure coherent states $\hat{\rho}_{\mathbf{z}}=|\mathbf{z}\rangle\langle\mathbf{z}|$ of the quantum harmonic oscillator. Note that the kernel functions and the pure coherent states are nonlinear images of $\mathbb{C}^{n}$ that sit as submanifolds of their respective (infinite-dimensional) Hilbert spaces. In addition, the Gaussian RKHS is often defined with a scale parameter $\lambda$, i.e. $K(\mathbf{w},\mathbf{z})=\exp(-|\mathbf{w}-\mathbf{z}|^{2}/\lambda)$. In the coherent state picture, this parameter can be absorbed by restoring physical dimensions to the phase space coordinates $\mathbf{z}$ and adjusting Planck’s constant $\hbar$. Equivalently, one can merely rescale the coordinates as $\mathbf{z}\mapsto\mathbf{z}/\sqrt{\lambda}$ to achieve the same effect. ∎

This isometry offers a precise characterization of the functions contained in the Gaussian RKHS over $\mathbb{C}^{n}$:

Corollary 2**.**

Any function in the Gaussian RKHS over $\mathbb{C}^{n}$ can be written uniquely as $f(\mathbf{z})=\langle\mathbf{z}|\hat{f}|\mathbf{z}\rangle$ for a Hilbert-Schmidt operator $\hat{f}$ on the state space of the $n$-mode quantum harmonic oscillator.

Proof.

In any RKHS, the linear span of the kernel functions is a dense subspace of the full Hilbert space [2]. Indeed, the Gaussian RKHS can be constructed from the kernel functions $K_{\mathbf{z}}(\,\cdot\,)$ as follows. Any function in the linear span of $K_{\mathbf{z}}$ can be written as:

[TABLE]

for coefficients $b_{i}$ and vectors $\mathbf{z}_{i}\in\mathbb{C}^{n}$. This expansion is unique because the kernel functions are linearly independent (due to strict positive definiteness of the Gaussian kernel). As the inner product is a sesquilinear form, the inner product between such functions is:

[TABLE]

The reproducing property of the kernel means that:

[TABLE]

as desired. The full Gaussian RKHS is defined as the Cauchy completion of these functions with respect to this inner product.

Now we apply the isometric map $K_{\mathbf{z}}(\,\cdot\,)\to\hat{\rho}_{\mathbf{z}}$. By linear independence, every function $f$ in the linear span of $K_{\mathbf{z}}(\,\cdot\,)$ is mapped uniquely to an operator $\hat{f}$ in the linear span of $\hat{\rho}_{\mathbf{z}}$:

[TABLE]

Apply the reproducing property of the kernel, and use the isometry of the inner product under the map $K_{\mathbf{z}}(\,\cdot\,)\to\hat{\rho}_{\mathbf{z}}$:

[TABLE]

Thus, for functions $f$ in the linear span of $K_{\mathbf{z}}(\,\cdot\,)$, we have a linear isometric map $f\to\hat{f}$ to a Hilbert-Schmidt operator $\hat{f}\in\mathcal{B}_{2}(\mathcal{H}_{F})$. Since the inner product is the same in both spaces by Theorem 2 (SI Appendix, A.2), the full Gaussian RKHS, which is the Cauchy completion of the span of $K_{\mathbf{z}}(\,\cdot\,)$, is isometrically mapped to the Cauchy completion of the span of $\hat{\rho}_{\mathbf{z}}$ in $\mathcal{B}_{2}(\mathcal{H}_{F})$. ∎

It turns out that this map $f\to\hat{f}$ is bijective, which establishes the equivalence of these Hilbert spaces:

Corollary 3**.**

There is an isometric isomorphism between the Hilbert-Schmidt space $\mathcal{B}_{2}(\mathcal{H}_{F})$ and the Gaussian RKHS over $\mathbb{C}^{n}$.

Proof.

The isometric equivalence between $K_{\mathbf{z}}(\,\cdot\,)$ and $\hat{\rho}_{\mathbf{z}}$ ensures that the $\hat{\rho}_{\mathbf{z}}$ are linearly independent: the Gram matrix is identical, so it must be strictly positive definite. So any operator $\hat{f}$ in the span of $\hat{\rho}_{\mathbf{z}}$ has a unique preimage $f$ in the span of $K_{\mathbf{z}}(\,\cdot\,)$. Likewise, since the inner product is the same in both spaces, any $\hat{f}$ in the Cauchy completion of the span of $\hat{\rho}_{\mathbf{z}}$ has a unique preimage $f$ in the Gaussian RKHS. Thus the map $f\to\hat{f}$ is injective.

Now we prove surjectivity: every Hilbert-Schmidt operator $\hat{A}\in\mathcal{B}_{2}(\mathcal{H}_{F})$ is the image of some function $A(\,\cdot\,)$ in the Gaussian RKHS. Equivalently, the linear span of the coherent states $\hat{\rho}_{\mathbf{z}}$ is dense in the entire Hilbert-Schmidt space $\mathcal{B}_{2}(\mathcal{H}_{F})$. Denote the closure of this linear span as $\mathcal{S}$. If $\mathcal{S}$ is a proper closed subspace of $\mathcal{B}_{2}(\mathcal{H}_{F})$, there must be a nonzero operator $\hat{A}$ that is orthogonal to every operator in $\mathcal{S}$; i.e., $\mathrm{tr}(\hat{S}\hat{A})=0$ for $\hat{S}\in\mathcal{{S}}$. Define the function $Q_{A}(\mathbf{z})=\langle\mathbf{z}|\hat{A}|\mathbf{z}\rangle$, which is the Hilbert-Schmidt inner product $\mathrm{tr}(\hat{\rho}_{\mathbf{z}}\hat{A})$ between $\hat{A}$ and the pure coherent state $\hat{\rho}_{\mathbf{z}}$. This is known as the $Q$-representation of the operator $\hat{A}$, which is unique for any $\hat{A}$ (see Lemma 4 below). If $\hat{A}\perp\mathcal{{S}}$, then its $Q$-representation is zero since it is orthogonal to every coherent state in $\mathcal{S}$. By uniqueness of the $Q$-representation, however, this requires that $\hat{A}=0$. So the span of $\hat{\rho}_{\mathbf{z}}$ is dense in the Hilbert-Schmidt space. We conclude that the Gaussian RKHS over $\mathbb{C}^{n}$ is isometrically isomorphic to the Hilbert-Schmidt operators on $\mathcal{{H}_{F}}$. ∎

Lemma 4**.**

For any bounded operator $\hat{A}$, the $Q$-representation $Q_{A}(\mathbf{z})=\langle\mathbf{z}|\hat{A}|\mathbf{z}\rangle$ over the coherent states $|\mathbf{z}\rangle$ is a unique function of $\hat{A}$. Moreover, for self-adjoint $\hat{A}$, $\langle\mathbf{z}|\hat{A}|\mathbf{z}\rangle$ is a bounded real entire function of both variables $\overline{\mathbf{z}}$ and $\mathbf{z}$.

Proof.

Define the function $\tilde{A}(\overline{\mathbf{w}},\mathbf{z})$ as:

[TABLE]

so the Q-representation is $Q_{A}(\mathbf{z})=\tilde{A}(\overline{\mathbf{z}},\mathbf{z})$. Expanding in the Fock basis using multi-index notation (eq. 1 in Lemma 3) gives:

[TABLE]

which converges for all $\overline{\mathbf{w}},\mathbf{z}$ by boundedness of $\langle\mathbf{n}|\hat{A}|\mathbf{m}\rangle\leq\|\hat{A}\|$ and rapid decrease of $(\mathbf{n}!\mathbf{m}!)^{-1/2}$. Thus $\tilde{A}(\overline{\mathbf{w}},\mathbf{z})$ is a sesquiholomorphic entire function (holomorphic in $\mathbf{z}$, anti-holomorphic in $\mathbf{w}$) whose series coefficients are defined by the Fock-basis matrix elements of $\hat{A}$. Since the series expansion of an entire function is unique, and a bounded operator is uniquely determined by its matrix elements in a complete basis, the map $\hat{A}\to\tilde{A}(\overline{\mathbf{w}},\mathbf{z})$ is injective.

The $Q$-representation is:

[TABLE]

which is an entire function of $\overline{\mathbf{z}},\mathbf{z}$ (treated as independent variables), and bounded since $\hat{A}$ is bounded. The map $\hat{A}\to\tilde{A}(\overline{\mathbf{z}},\mathbf{z})$ is injective by the same argument as above: uniqueness of series coefficients for an entire function, and uniqueness of matrix elements in a complete basis. Thus the $Q$-representation is unique. Finally, if $\hat{A}$ is self-adjoint, its Q-representation $\langle\mathbf{z}|\hat{A}|\mathbf{z}\rangle$ must be real by the Hermitian symmetry of diagonal matrix elements. ∎

A.3 Theorem 3

Every QPM completely metrizes $\mathcal{P}_{\mathrm{w}}(X)$.

Proof.

We assume that we have a QPM on $\mathcal{P}_{\mathrm{w}}(X)$. As defined in the main text, this requires a function $\phi:X\to\mathcal{S}$ that maps points $x\in X$ to pure states $|x\rangle\langle x|\equiv\hat{\rho}_{x}$ on a Hilbert space. We require that $\phi$ is a closed embedding, so its image $\mathcal{C}$ is closed and homeomorphic to $X$, and also that its kernel $k(x,y)=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$ is characteristic. All of these properties are preserved if we normalize the states $\langle x|x\rangle=\mathrm{tr}(\hat{\rho}_{x})=1$, so we assume normalized states without loss of generality; this gives a normalized kernel $k(x,x)=1$.

The barycenter map is:

[TABLE]

which is well-defined as a strong (Bochner) integral since the map $x\to\hat{\rho}_{x}$ is continuous from $X$ to the separable space $\mathcal{B}_{1}(\mathcal{H})$, and the trace norm $\|\hat{\rho}_{x}\|_{1}=1$. Since $\mu$ is a probability measure, it is straightforward to show that $\hat{\mu}$ is an element of the quantum state space $\mathcal{S}$: it is a positive operator with unit trace. By definition, the barycenter map is injective for a characteristic kernel.

By assumption, the image $\phi(X)\equiv\mathcal{C}$ is a (norm-)closed subset of $\mathcal{S}$ that is homeomorphic to $X$. Labeling operators in $\mathcal{C}$ by their corresponding points in $X$ and viewing $\mu$ as the pushforward measure on $\mathcal{C}$, we can write the barycenter map as an integral on $\mathcal{C}$:

[TABLE]

which is now viewed as a measure on the closed set $\mathcal{C}\subset\mathrm{extr}\,\mathcal{S}$ (recall that pure states are the extreme points of the convex quantum state space $\mathcal{S}$). Measures on $\mathcal{C}$ inherit the weak topology from measures on $X$, and vice versa, via the homeomorphism $\phi$. The barycenter map $T(\mu)$ is continuous from $\mathcal{P}_{\mathrm{w}}(\mathcal{C})$ to $\mathcal{S}$ with the trace norm [6], it is convex (affine), and its image is contained in the closed convex hull of $\mathcal{C}$, denoted $\overline{\mathrm{co}}(\mathcal{C})$.

Proposition 4 from [6] states that the quantum state space $\mathcal{S}$ is a $\mu$-compact set, which means that the barycenter map on $\mathcal{S}$ is a proper map (the preimage of a compact set is compact). This property extends to closed subsets such as $\mathcal{C}$: the preimage of any compact subset of $\overline{\mathrm{co}}(\mathcal{C})$ is a compact subset of $\mathcal{P}_{\mathrm{w}}(\mathcal{C})$. The key results of Choquet theory, which usually apply only to compact sets, extend to $\mu$-compact sets. For instance, Proposition 2 of [6] shows that the barycenter map $T(\mu)$ is surjective onto $\overline{\mathrm{co}}(\mathcal{C})$: every element of $\overline{\mathrm{co}}(\mathcal{C})$ is the barycenter of some measure on $\mathcal{C}$. Since the barycenter map is injective (by assumption that the kernel $k(x,y)=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$ is characteristic), these measures are unique. Thus we have a continuous bijection from measures on $\mathcal{C}$ to quantum states in the closed convex set $\overline{\mathrm{co}}(\mathcal{C})$.

To show that the barycenter map has a continuous inverse, note that any continuous proper map $f:X\to Y$ between metrizable spaces is closed: it maps closed sets to closed sets. (If $U\subseteq X$ is closed but $f(U)$ is not, pick $y_{n}\in f(U)$ converging to $y\notin f(U)$. The set $\{y_{n}\}\cup\{y\}$ is compact, so by properness, $\{x_{n}:f(x_{n})=y_{n}\}$ has a convergent subsequence $x_{n_{k}}\to x\in U$. Continuity gives $f(x)=y$, contradicting $y\notin f(U)$. Thus $f(U)$ must be closed.) Since any continuous closed bijection is a homeomorphism, the barycenter map is a homeomorphism between the metrizable spaces $\mathcal{P}_{\mathrm{w}}(\mathcal{C})$ and $\overline{\mathrm{co}}(\mathcal{C})$. Note that this construction works for any characteristic kernel $k(x,y)=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$; homeomorphism follows from the topological properties of the barycenter map on $\mathcal{S}$.

Finally, since the embedding map $\phi$ is a homeomorphism from $X$ onto its image $\mathcal{C}$, the spaces $\mathcal{P}_{\mathrm{w}}(\mathcal{C})$ and $\mathcal{P}_{\mathrm{w}}(X)$ are likewise homeomorphic. Thus the barycenter map establishes a homeomorphism between $\mathcal{P}_{\mathrm{w}}(X)$ and $\overline{\mathrm{co}}(\mathcal{C})$. Since this is a closed set in the complete metric space $\mathcal{B}_{1}(\mathcal{H})$, its metric (the trace distance) completely metrizes $\mathcal{P}_{\mathrm{w}}(X)$. ∎

As a characteristic kernel is essential for this construction, the following two lemmas highlight important features of such kernels:

Lemma 5**.**

For any normalized kernel $k(x,y)=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$, $k$ is characteristic if and only if it is integrally strictly positive definite, i.e. for all finite complex measures $m\neq 0$ we have:

[TABLE]

Proof.

First we show that being integrally s.p.d. is equivalent to having an injective barycenter map for all complex measures. Suppose $\iint k(x,y)\,d\overline{m}(x)\,dm(y)=0$ for some nonzero $m$. Choose an arbitrary pair of distinct measures $m_{1},m_{2}$ with $m_{1}-m_{2}=m$ and calculate the Hilbert-Schmidt distance between the barycenter embeddings of these measures (Bochner integrals):

[TABLE]

so the barycenter map is not injective for $m_{1},m_{2}$. (The trace and integral can be interchanged using Fubini-Tonelli by boundedness of $|\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})|$ and finite total variation of $|m|$). Conversely, if the barycenter map is not injective for some pair $m_{1},m_{2}$, then the same calculation shows that for $m=m_{1}-m_{2}$ we must have $\iint k(x,y)\,d\overline{m}(x)\,dm(y)=0$, so the kernel is not integrally s.p.d.

Now we prove the lemma. Since the barycenter map with an integrally s.p.d. kernel is injective for all complex measures, it is injective for probability measures. Thus integrally s.p.d. implies characteristic. Although in general the converse does not hold [2], it does hold for kernels of the form $k(x,y)=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$. Suppose we have two distinct complex measures $m_{1},m_{2}$. Since the states are normalized, $\langle x|x\rangle=\mathrm{tr}(\hat{\rho}_{x})=1$, the trace of the barycenter embedding of any measure is equal to its total measure $m(X)$:

[TABLE]

So the embedded measures will be distinct (different traces) unless they have the same total measure, $m_{1}(X)=m_{2}(X)$. Note that this argument is not available in the usual RKHS context: the trace of the barycenter embedding map is equivalent to integration against the constant witness function $\mathbf{1}$, which may not belong to the RKHS. As pointed out in the main text, for many commonly-used kernels such as the Gaussian or Laplacian kernels, all RKHS functions vanish at infinity so constant functions are not included. Using the trace embedding, on the other hand, the identity operator maps to the constant function $\mathbf{1}$, which is always included as a dual function (witness function).

If $m_{1}(X)=m_{2}(X)$, define $m=m_{1}-m_{2}$ so that $m(X)=0$. Assume that the barycenter map is not injective for $m_{1},m_{2}$, i.e. $T(m_{1})-T(m_{2})=T(m)=0$. The complex measure $m$ has Jordan decomposition $m=(\mu_{1}-\mu_{2})+i(\mu_{3}-\mu_{4})$, and since $m(X)=0$, these positive measures must have total measure $\mu_{1}(X)=\mu_{2}(X)$ and $\mu_{3}(X)=\mu_{4}(X)$. Unless they vanish identically, each pair can be normalized to yield probability measures for which the barycenter map is injective. Thus, $T(\mu_{1})-T(\mu_{2})=T(\mu_{1}-\mu_{2})=0$ forces $\mu_{1}=\mu_{2}$, and likewise for $\mu_{3}=\mu_{4}$, which gives $m=m_{1}-m_{2}=0$, contradicting the assumption that $m_{1},m_{2}$ are distinct. So the barycenter map is indeed injective for all complex measures, which means that the kernel $k(x,y)=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$ is integrally s.p.d. ∎

The next lemma gives a simple construction for kernels of the desired form starting with any integrally s.p.d. kernel $k^{\prime}(x,y)=\langle x|y\rangle$ on the underlying Hilbert space: simply use the norm-squared kernel $k=|k^{\prime}|^{2}=\langle y|x\rangle\langle x|y\rangle=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$. We show that the resulting kernel is integrally strictly positive definite:

Lemma 6**.**

If the kernel $k^{\prime}(x,y)=\langle x|y\rangle$ is integrally strictly positive definite, then so is the norm-squared kernel $k=|k^{\prime}|^{2}=\langle y|x\rangle\langle x|y\rangle=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$.

Proof.

On the dual $\mathcal{H}^{*}$ of the underlying Hilbert space $\mathcal{H}$, the barycenter map $B(m)$ for a (generally complex) measure $m$ is:

[TABLE]

where $\langle m|$ represents the embedded dual vector corresponding to the measure $m$. Embedding as dual vector $\langle m|\in\mathcal{H}^{*}$ instead of as a vector $|m\rangle$ is conventional because it allows the expected value of a function $|f\rangle$ in the RKHS to be calculated as $\langle m|f\rangle$.

Following the same argument in Lemma 5, since $k^{\prime}(x,y)=\langle x|y\rangle$ is integrally s.p.d., $B$ is injective [2]. Thus, for any two distinct measures $m_{1},m_{2}$ there is a witness vector $|w\rangle$ such that $\langle m_{1}|w\rangle\neq\langle m_{2}|w\rangle$. Since $\langle x|w\rangle$ is a continuous function of $x$, define new complex measures by setting $dM_{1}\equiv\langle x|w\rangle dm_{1}$ and $dM_{2}\equiv\langle x|w\rangle dm_{2}$ and note that they have different total measures:

[TABLE]

so $M_{1}\neq M_{2}$ (and both total measures are finite).

Now assume for contradiction that the norm-squared kernel $k=|k^{\prime}|^{2}=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$ is not integrally s.p.d. By Lemma 5, this means that its barycenter map

[TABLE]

fails to be injective, i.e. for some distinct measures $m_{1},m_{2}$ we have $T(m_{1})=T(m_{2})$. Following the argument above and applying the operator $T(m_{1})-T(m_{2})$ to the witness vector $|w\rangle$ gives:

[TABLE]

Taking the dual, we have:

[TABLE]

Since $B$ is injective, this requires $M_{1}=M_{2}$, which contradicts our earlier conclusion that $M_{1}\neq M_{2}$. Hence $k$ must be integrally strictly positive definite. ∎

Note that although the norm-squared kernel inherits the property of being integrally s.p.d. from its base kernel, the same is not true for being characteristic. As a counterexample, consider a two-element space $X=\{a,b\}$ and define the kernel

[TABLE]

For a probability measure $P=(P_{a},P_{b})$ with $P_{a}+P_{b}=1$, the embedding induced by $k^{\prime}$ is:

[TABLE]

In particular,

[TABLE]

[TABLE]

Thus, the embedding can be written as $B(P)=(P_{a}-P_{b},P_{b}-P_{a})$. Since the difference $P_{a}-P_{b}$ uniquely characterizes the measure, this embedding is injective and hence $k^{\prime}$ is characteristic.

Now, consider the squared kernel $k(x,y)=|k^{\prime}(x,y)|^{2}$ which has all entries equal to 1. Its embedding maps any probability measure $P=(P_{a},P_{b})$ to:

[TABLE]

Since all measures are mapped to the same point, the squared kernel is not characteristic. In this case, neither kernel is injective on the entire space of finite measures (as opposed to probability measures). For instance, the kernel $k^{\prime}$ is not injective—it gives the same embedding for measures $(1,1)$ and $(2,2)$.

Finally, note that Lemma 6 can be viewed as a continuous, integral generalization of the Schur product theorem. Recall that this theorem states that the entrywise Hadamard product of two (strictly) positive definite matrices is (strictly) positive definite. As a special case, if a matrix $A_{ij}$ is s.p.d. ($\langle v|A|v\rangle>0$ for all nonzero $v$), then its Hadamard product with its own conjugate, $B_{ij}=|A_{ij}|^{2}$, is also s.p.d. In our setting, integral strict positive definiteness corresponds to positivity under integration against nonzero complex measures, rather than positivity under summation. Thus, the claim presented here is a generalization of this special case of the Schur product theorem, replacing finite sums with integrals and discrete vectors with measures. Indeed, one can extend this proof to the general case of two integrally s.p.d. kernels $k_{1},k_{2}$ and show that their product $k_{1}k_{2}$ is also integrally s.p.d.

A.4 Theorem 4

Every Polish space has a QPM.

Proof.

Given a Polish space $X$—a separable, completely metrizable space—we seek a map $\phi:X\to\mathcal{S}$ that maps points $x\in X$ to pure states $|x\rangle\langle x|\equiv\hat{\rho}_{x}$ on a Hilbert space. We require $\phi$ to be a closed embedding with a characteristic kernel $k(x,y)=\mathrm{tr}(\hat{\rho}_{x}\hat{\rho}_{y})$.

First, choose a metric $d$ that completely metrizes $X$. Embed $(X,d)$ isometrically as a closed subset of the Urysohn universal metric space, which is homeomorphic to a separable Hilbert space $\mathcal{H}$ [7]. This closed embedding maps $x\in X$ to vectors $\mathbf{x}\in\mathcal{H}$.

Following the construction in Lemma 3, extended to infinite dimensions as in Section 10 of [4], build the Fock space $\mathcal{H}_{F}$ over $\mathcal{H}$. Each vector $\mathbf{x}$ in the original Hilbert space is mapped to a coherent state vector $|\mathbf{x}\rangle$ in $\mathcal{H}_{F}$, with inner product:

[TABLE]

where $(\mathbf{x},\mathbf{y})$ is the inner product in $\mathcal{H}$.

Define the map $\phi$ in terms of these coherent states as:

[TABLE]

where $\mathcal{S}$ is the quantum state space over $\mathcal{H}_{F}$; this gives a closed embedding of $X$ as pure states in $\mathcal{S}$. The kernel is:

[TABLE]

which is the Gaussian kernel for embedded vectors $\mathbf{x},\mathbf{y}\in\mathcal{H}$. The Gaussian kernel on a separable Hilbert space is characteristic (Theorem 3.1 of [5]). Thus the map $\phi$ has all of the properties required for Theorem 3 (SI Appendix A.3): it defines a QPM that completely metrizes the Polish space $X$. If one already has a closed embedding of $X$ into a vector space (as is often the case in machine learning), the Gaussian kernel—or, indeed, any integrally s.p.d. kernel—immediately gives a QPM for this space. ∎

A.5 Theorem 5

Dual functions for the Fock (coherent state) QPM are dense in $BUC(\mathbb{R}^{n})$: they can uniformly approximate any bounded, uniformly continuous function on $\mathbb{R}^{n}$.

This proof is somewhat long and technical, so we first offer a short overview. For reasons that will be clear later, we consider functions on $\mathbb{R}^{2n}$, represented as $F(\mathbf{x},\mathbf{y})$ for $\mathbf{x},\mathbf{y}\in\mathbb{R}^{n}$; any function on $\mathbb{R}^{m}$ can trivially be embedded in $\mathbb{R}^{2n}$ for $2n\geq m$. The proof proceeds as follows.

We start by considering real functions $F(\mathbf{x},\mathbf{y})$ that are bounded on $\mathbb{R}^{2n}$ and band-limited (i.e. distributional Fourier transform has compact support). Using a Paley-Wiener type result (Lemma 7 below), such a function extends to an entire function of exponential type on $\mathbb{C}^{2n}$, with exponential growth depending only on the imaginary parts of its arguments. 2. 2.

Using the Q-representation (eq. 7 in Lemma 4), this (real) entire function $F(\mathbf{x},\mathbf{y})$ determines the unique Fock-basis matrix elements of a Hermitian operator $\hat{F}$ on the Fock space $\mathcal{H}_{F}$, with $F(\mathbf{x},\mathbf{y})=Q_{F}(\mathbf{x}+i\mathbf{y})=\langle\mathbf{x}+i\mathbf{y}|\hat{F}|\mathbf{x}+i\mathbf{y}\rangle$ (these are coherent-state vectors for complex $\mathbf{z}=\mathbf{x}+i\mathbf{y}$, with $\mathbf{x},\mathbf{y}\in\mathbb{R}^{n}$). 3. 3.

Using the Segal-Bargmann-Fock correspondence (Lemma 3), the operator $\hat{F}$ on $\mathcal{H}_{F}$ is isometrically isomorphic to an integral operator on $L^{2}(\mathbb{C}^{n})$ whose integral kernel is bounded by the exponential growth bounds for $F$. 4. 4.

Using the Schur test [8], this integral kernel defines a bounded operator on $L^{2}(\mathbb{C}^{n})$, which implies that the operator $\hat{F}$ is indeed a bounded, self-adjoint operator on the Fock space $\mathcal{H}_{F}$. Proving boundedness is the key step in this theorem. 5. 5.

Since the bounded self-adjoint operators are dual to the self-adjoint trace-class operators, every bounded, band-limited function is a dual function, with $F(\mathbf{x},\mathbf{y})=\langle\mathbf{x}+i\mathbf{y}|\hat{F}|\mathbf{x}+i\mathbf{y}\rangle$. 6. 6.

Finally, using the results of [9], these dual functions can uniformly approximate any bounded, uniformly continuous function on $\mathbb{R}^{2n}$.

We begin with a technical lemma of Paley-Wiener type:

Lemma 7**.**

Any band-limited real function $f(\mathbf{x})$ with $|f(\mathbf{x})|\leq M$ on $\mathbb{R}^{n}$ extends to an entire function $f(\mathbf{z})$ on $\mathbb{C}^{n}$ with:

[TABLE]

for positive constants $\tau_{j}$.

Proof.

Begin with the one-dimensional case. By the Paley-Wiener-Schwartz theorem, any real function $f(x)$ that is band-limited (distributional Fourier transform supported in $[-\tau,\tau]$) extends to an entire function of exponential type $\tau$. Then, by Proposition 1.1 in [9], if in addition $|f(x)|\leq M$ on the real axis, then for all $z\in\mathbb{C}$ we have $|f(z)|\leq Me^{\tau|\Im z|}$. (We adopt the Fourier convention $\widehat{f}(\xi)=\int e^{-ix\xi}f(x)\,dx$.)

For the multivariate case, if a function $f(\mathbf{x})$ on $\mathbb{R}^{n}$ has a compactly-supported distributional Fourier transform, it extends to an entire function $f(\mathbf{z})$ of exponential type on $\mathbb{C}^{n}$, i.e., there are positive constants $A,B$ with $|f(\mathbf{z})|\leq Ae^{B|\mathbf{z}|}$. Such a function is of exponential type in each complex coordinate alone, treating the other coordinates as fixed. That is, if $f(\mathbf{z})$ is of exponential type, then for $z_{k}$ alone we have:

[TABLE]

[TABLE]

[TABLE]

for a new constant $A^{\prime}$ that depends on all coordinates except $z_{k}$, and these coordinates are held fixed. Thus $f(z_{1},\ldots,z_{k},\ldots,z_{n})$, treated as a single-variable function of $z_{k}$, is of exponential type in that variable alone.

Now, a function $f(\mathbf{x})$ as described in the lemma extends to an entire function of exponential type $f(\mathbf{z})$. Write each complex coordinate as $z_{j}=x_{j}+iy_{j}$. Since $f$ is bounded on $\mathbb{R}^{n}$, we have:

[TABLE]

Treating $f$ as a function of the first coordinate alone, with $x_{2},\ldots,x_{n}$ fixed, we have a function of exponential type that is bounded by $M$ on the real axis (i.e. for all $x_{1}\in\mathbb{R}$). Applying the single-variable result above, and extending the first coordinate to $x_{1}+iy_{1}\in\mathbb{C}$, we have, for some constant $\tau_{1}$,

[TABLE]

Now treat $f$ as a function of the second coordinate, with $x_{1}+iy_{1}$ and $x_{3},\ldots,x_{n}$ fixed; we have a function of exponential type that is bounded for all $x_{2}\in\mathbb{R}$ by $M\exp(\tau_{1}|y_{1}|)$. Extending the second coordinate to $\mathbb{C}$ and applying the single-variable result gives:

[TABLE]

for some constant $\tau_{2}$. Iterating in this fashion proves the lemma. With additional effort (not needed for this proof), one can show that the constants $\tau_{j}$ define the support of the distributional Fourier transform as the box $[-\tau_{1},\tau_{1}]\times[-\tau_{2},\tau_{2}]\times\ldots\times[-\tau_{n},\tau_{n}]$. ∎

We need one more lemma that connects the analytic continuation of $F(\mathbf{x},\mathbf{y})$ on $\mathbb{C}^{2n}$ with the normalized off-diagonal coherent-state matrix elements $\langle\mathbf{w}|\hat{F}|\mathbf{z}\rangle$ of the corresponding operator $\hat{F}$:

Lemma 8**.**

Given a real entire function $F$ of exponential type, with $F(\mathbf{x},\mathbf{y})=Q_{F}(\mathbf{x}+i\mathbf{y})=\langle\mathbf{x}+i\mathbf{y}|\hat{F}|\mathbf{x}+i\mathbf{y}\rangle$ for a Hermitian operator $\hat{F}$, the extension of $F$ to complex-valued $\mathbf{x},\mathbf{y}$ is given by:

[TABLE]

where $\overline{\mathbf{x}},\overline{\mathbf{y}}$ represent the complex conjugates of $\mathbf{x},\mathbf{y}\in\mathbb{C}^{n}$.

Proof.

Recall from eq. 7 of Lemma 4 that the Q-representation of the operator $\hat{F}$ has the absolutely convergent power series:

[TABLE]

From the usual perspective of complex analyticity, this is a peculiar power series because it is entire for both $\mathbf{z}$ and $\overline{\mathbf{z}}$ independently. In some ways, this expression looks more natural if we express it in terms of the independent real and imaginary parts of $\mathbf{z}=\mathbf{x}+i\mathbf{y}$:

[TABLE]

(we notate $\mathbf{x}^{2}=\mathbf{x}\cdot\mathbf{x}=x_{1}^{2}+\ldots+x_{n}^{2}$, and same for $\mathbf{y}^{2}$; these are holomorphic squares not norm-squares.) Written in this form, it is evident that the right-hand side is a real analytic function of $\mathbf{x},\mathbf{y}$ that extends to an entire function for complex $\mathbf{x},\mathbf{y}$. But the left-hand side cannot be extended to complex arguments, because the expression $\mathbf{x}+i\mathbf{y}$ for complex $\mathbf{x},\mathbf{y}$ does not capture the $4n$ real degrees of freedom contained in the two complex vectors $\mathbf{x},\mathbf{y}\in\mathbb{C}^{n}$.

Consider instead the normalized off-diagonal matrix elements defined in eq. 5 of Lemma 4 as:

[TABLE]

where we define the invertible map $\mathbf{w}=\overline{\mathbf{x}}+i\overline{\mathbf{y}},\;\;\mathbf{z}=\mathbf{x}+i\mathbf{y}$. Using eq. A.2 from Lemma 4, we obtain (note $\overline{\mathbf{w}}=\mathbf{x}-i\mathbf{y}$; we use multi-index notation $\mathbf{n}!=n_{1}!n_{2}!\ldots n_{n}!$ and $\mathbf{z}^{\mathbf{n}}=z_{1}^{n_{1}}z_{2}^{n_{2}}\ldots z_{n}^{n_{n}}$):

[TABLE]

Note that the values of the entire function $F(\mathbf{x},\mathbf{y})$ for $\mathbf{x},\mathbf{y}\in\mathbb{C}^{n}$ are determined by its values on the real subspace. Likewise, the matrix elements $\langle\mathbf{w}|\hat{F}|\mathbf{z}\rangle$ for any pair of coherent states are determined by their values on the diagonal subspace $\mathbf{w}=\mathbf{z}$. Moreover, since the coherent states are a total set that span the Fock space $\mathcal{H}_{F}$, these matrix elements completely define the operator $\hat{F}$. ∎

We now prove the main theorem.

Proof.

(Theorem 5, SI Appendix, A.5) First, we show that, for the Fock embedding, the dual functions include any function $F(\mathbf{x},\mathbf{y})$ that is bounded on $\mathbb{R}^{2n}$ and band-limited (i.e. its distributional Fourier transform has compact support). We start with $n=1$ for simplicity.

Suppose we are given a real, band-limited function that is bounded, $F(x,y)\leq M$ for $x,y\in\mathbb{R}$. We seek an operator $\hat{F}$ on the Fock space $\mathcal{H}_{F}$ such that, for $z=x+iy$, we have:

[TABLE]

If such an operator is bounded and self-adjoint, then it defines a dual function as desired. Using Lemma 7, we know that $F$ is an entire function with an everywhere-convergent power series in $x,y$. The product of two entire functions is entire, so the product $F(x,y)e^{(x^{2}+y^{2})}$ has an everywhere-convergent power series with unique coefficients $A_{nm}$:

[TABLE]

Defining $z=x+iy$ and substituting $x=(z+\overline{z})/2$, $y=(z-\overline{z})/(2i)$ we get:

[TABLE]

By expanding the binomials, rearranging terms, and adding a convenient combinatorial factor $(n!m!)^{-1/2}$, we can rewrite the power series in terms of $z,\overline{z}$ with new coefficients $F_{nm}$:

[TABLE]

If we treat $z,\overline{z}$ as independent variables, this expression involves the composition of entire functions, i.e. $F[x(z,\overline{z}),y(z,\overline{z})]e^{\overline{z}z}$, so the resulting power series is entire in $z,\overline{z}$ independently. Solving for $F(x,y)$ and comparing with the power-series expansion of $Q_{F}(z)=\langle z|\hat{F}|z\rangle$ in eq. 9 of Lemma 8, we see:

[TABLE]

if and only if the Fock-basis matrix elements of $\hat{F}$ are equal to the power-series coefficients, $\langle n|\hat{F}|m\rangle=F_{nm}$; since $F(x,y)$ is real, these coefficients must be Hermitian, $F_{nm}=\overline{F_{mn}}$. Thus, we define the operator $\hat{F}$ to have these matrix elements. At this stage, it is not at all clear that they define a bounded operator. Nonetheless, $\hat{F}$ with these specified Fock-basis matrix elements will satisfy:

[TABLE]

We now extend $F(x,y)$ to complex arguments. From Lemma 7, this extension is an entire function of exponential type, with growth bounded by the imaginary parts of its arguments:

[TABLE]

From eq. 8 of Lemma 8, this complex extension is given by the normalized off-diagonal coherent-state matrix elements of the same operator $\hat{F}$, defining $w=\overline{x}+i\overline{y}$ and $z=x+iy$:

[TABLE]

It is essential here that the proof of Lemma 8 does not depend on $\hat{F}$ being a bounded operator, only that its Fock-basis matrix elements match the power-series coefficients $F_{nm}$ as required.

To prove boundedness, we use the Schur test [8], in the following form. Given a Hilbert space $L^{2}(X,\mu)$, we define a linear operator $T$ that acts on $f\in L^{2}$ via the integral kernel $K(x,y)$ as:

[TABLE]

The Schur test states that this operator will be bounded (from $L^{2}(X,\mu)$ to itself) if there exist two constants $A,B$ such that:

[TABLE]

Furthermore, the operator norm will satisfy $\|T\|\leq\sqrt{AB}$.

In what follows, we will often convert between two sets of complex variables, so it may be helpful to have a table:

Note that these definitions are consistent with functions that are holomorphic in $(x,y)$ and sesqui-holomorphic in $(w,z)$, i.e., holomorphic in $\overline{w}$ and $z$. In particular, we have $\overline{w}=x-iy$.

To use the Schur test, we apply the unitary Segal-Bargmann-Fock correspondence (eq. 3 of Lemma 3). Consider the action of an operator $\hat{A}$ on a vector $|\psi\rangle$ in the Fock space $\mathcal{H}_{F}$. Each such vector corresponds to a complex function $\psi(z)=\langle z|\psi\rangle$ belonging to $L^{2}(\mathbb{C})$ with the following squared norm:

[TABLE]

where we recognize the resolution of the identity operator (in the weak operator topology, see eq. 4 of Lemma 3):

[TABLE]

The action of $\hat{A}$ on $|\psi\rangle$ is given by $\hat{A}|\psi\rangle$. Taking the inner product with the coherent-state vector $\langle w|$ and inserting the resolution of the identity, we get:

[TABLE]

On the left, we have a complex function of $w$ that corresponds to the Segal-Bargmann-Fock representation of the vector $\hat{A}|\psi\rangle$; let’s label this function as $T_{A}[\psi](w)$. Inside the integral, we recognize that $\langle z|\psi\rangle$ is a complex function of $z$ that corresponds to the Segal-Bargmann-Fock representation of the vector $|\psi\rangle$; we already identified this function with $\psi(z)$. Furthermore, we identify the integral kernel $K_{A}(w,z)=\langle w|\hat{A}|z\rangle$. These substitutions give:

[TABLE]

This derivation shows that the coherent-state matrix elements $\langle w|\hat{A}|z\rangle$ define the values of the integral kernel $K_{A}(w,z)$ for the Segal-Bargmann-Fock representation of the operator $T_{A}$ corresponding to the linear operator $\hat{A}$ on the Fock space.

Thus, as expected, we have an isometric isomorphism between vectors $|\psi\rangle$ in the Fock space and functions $\psi(z)$ in $L^{2}(\mathbb{C})$. Likewise, we have an isometric isomorphism between linear operators $\hat{A}$ on the Fock space and operators $T_{A}$ defined by integral kernels $K_{A}(w,z)=\langle w|\hat{A}|z\rangle$ on $L^{2}(\mathbb{C})$. Moreover, from eq. A.5 of Lemma 8, we know that:

[TABLE]

With this isomorphism, we can apply the Schur test to show that the operator with integral kernel $K_{F}(w,z)=\langle w|\hat{F}|z\rangle$ is bounded on $L^{2}(\mathbb{C})$; this implies that the isometric operator $\hat{F}$ is bounded on the Fock space. The Schur test requires the magnitude of the kernel, $|K_{F}(w,z)|$, which we calculate as:

[TABLE]

Now note that $z-w=(x+iy)-(\overline{x}+i\overline{y})=2i(\Im x)-2(\Im y)$, so $|z-w|^{2}=4(\Im x)^{2}+4(\Im y)^{2}$. Putting it all together, we have:

[TABLE]

To apply the Schur test, we first seek a constant $A$ that bounds the following integral $I(w)$ for all $w\in\mathbb{C}$:

[TABLE]

Since the complex Lebesgue measure $d^{2}z=d\Re z\,d\Im z$ is translation-invariant, and we are integrating over all of $\mathbb{C}$, for fixed $w$ we can shift the measure and write this integral as:

[TABLE]

where we use $z-w=2i(\Im x)-2(\Im y)$, and the Jacobian gives $|\det J|=4$. Substitute our earlier expression for $|K_{F}(w,z)|$:

[TABLE]

and apply Lemma 7, $|F(x,y)|\leq M\exp(\tau_{x}|\Im x|+\tau_{y}|\Im y|)$, to obtain:

[TABLE]

for constants $\tau_{x}$ and $\tau_{y}$ that represent half of the Fourier bandwidth in each coordinate. We now separate variables and rewrite these integrals in standard form over two half-lines:

[TABLE]

These Gaussian integrals can be evaluated in closed form by completing the square to obtain:

[TABLE]

These calculations show convergence for any $\tau_{j}$, so we have achieved our goal of bounding the integral $I(w)$ by a constant $A$. Since $\mathrm{erf}(x)<1,\forall x$, we can simplify the final expression at the expense of a slightly looser bound by simply setting $\text{erf}(x)=1$:

[TABLE]

For the other part of the Schur test, we seek a constant B that bounds the following integral for all $z\in\mathbb{C}$:

[TABLE]

Since the kernel $K_{F}(w,z)=\langle w|\hat{F}|z\rangle$ is Hermitian, we have $|K_{F}(w,z)|=|K_{F}(z,w)|$ and thus the bound $B$ is the same as the earlier bound $A$; note the Schur test gives $\|T_{F}\|\leq\sqrt{AB}$.

Putting this all together, we bound the operator norm of the integral operator $T_{F}$ with kernel $K_{F}(w,z)$ as:

[TABLE]

Since the Segal-Bargmann-Fock transform is isometric, this is also the operator norm of the operator $\hat{F}$ defined on the Fock space through the matrix elements in its power-series expansion. Thus $\hat{F}$ is a bounded, self-adjoint operator as desired, and the function $F(x,y)=\langle x+iy|\hat{F}|x+iy\rangle$ defines a dual function in the QPM.

Extension to higher dimensions is straightforward. Given a real, bounded function $F(\mathbf{x},\mathbf{y})\leq M$ on $\mathbb{R}^{2n}$ that is band-limited, the entire argument follows by replacing the complex variables $x,y,w,z\in\mathbb{C}$ with vectors $\mathbf{x},\mathbf{y},\mathbf{w},\mathbf{z}\in\mathbb{C}^{n}$, and indices $n,m$ with multi-indices $\mathbf{n},\mathbf{m}$. From eq. 4 of Lemma 3, the resolution of the identity in $\mathbb{C}^{n}$ has a numerical factor of $1/\pi^{n}$. In calculating the bound for the Schur test on $L^{2}(\mathbb{C}^{n})$, the Jacobian contributes a factor $4^{n}$, and the product of $2n$ Gaussian integrals gives a factor of $(2\pi)^{n}$ (using the looser bound that sets $\mathrm{erf}(x)=1$). The final bound for the norm of $\hat{F}$, for $F(\mathbf{x},\mathbf{y})=\langle\mathbf{x}+i\mathbf{y}|\hat{F}|\mathbf{x}+i\mathbf{y}\rangle$ on $\mathbb{R}^{2n}$, is:

[TABLE]

where we recall from Lemma 7 that the $\tau_{j}$’s define the support of the Fourier transform as the box $[-\tau_{1},\tau_{1}]\times[-\tau_{2},\tau_{2}]\times\ldots\times[-\tau_{2n},\tau_{2n}]$.

This is a loose bound: the constant function $F(\mathbf{x},\mathbf{y})=M$, with Fourier transform supported at the origin, is represented by a multiple of the identity operator, $\hat{F}=M\hat{I}$, which has operator norm $|M|$, not $8^{n}|M|$. Nonetheless, any bound is sufficient to prove that $\hat{F}$ is a bounded operator. And this bound may yield some intuition about the IPM witness functions for the Fock embedding; these functions must have operator representations with norm $\|\hat{F}\|\leq 1$. Using the Schur bound, we see that these witness functions can have arbitrary support on $\mathbb{R}^{2n}$, but their Fourier bandwidth may be limited as the operator norm could grow exponentially with the square of the bandwidth.

We complete the proof by using Theorem 3.1 from [9]: any bounded, uniformly continuous function on $\mathbb{R}$ can be uniformly approximated by a bounded real entire function of exponential type. The cited proof is a straightforward exercise in analytic interpolation, and extends readily to functions on $\mathbb{R}^{n}$. And the class of bounded real entire functions of exponential type are precisely the bounded, band-limited functions—which are all contained as dual functions in the Fock embedding QPM. Thus the dual functions are dense in $BUC(\mathbb{R}^{n})$: they can uniformly approximate any bounded, uniformly continuous function on $\mathbb{R}^{n}$. ∎

It is instructive to compare the Schur bound for the QPM, which uses the operator norm, with the bound for MMD using the Gaussian RKHS. Corollary 2 shows that the RKHS norm is simply the Hilbert-Schmidt norm of the operator embedding in the Fock space $\mathcal{H}_{F}$. The space of Hilbert-Schmidt operators is much smaller than the space of bounded operators. Every Hilbert-Schmidt operator must be compact, so even the identity operator is not Hilbert-Schmidt. For the Fock space on $\mathbb{C}$, if we have a self-adjoint operator $\hat{f}$ that corresponds to a function $f(z)=\langle z|\hat{f}|z\rangle$, the squared Hilbert-Schmidt norm is given by:

[TABLE]

where the first line expresses the trace in terms of coherent states, the second inserts the resolution of the identity, and the third identifies the kernel $K_{f}(w,z)=\langle w|\hat{f}|z\rangle$. The operator $\hat{f}$ is Hilbert-Schmidt (and defines a function in the Gaussian RKHS) if and only if this kernel is square-integrable over all $w$ and $z$. In contrast, the Schur test requires only that the kernel be absolutely integrable over $z$ for each fixed $w$, and vice versa—a much weaker requirement. As a result, even on a compact space, the IPM witness functions for the QPM are a much broader class of functions compared with the RKHS functions used in MMD. For instance, while all functions in the Gaussian RKHS must vanish at infinity, we have shown that all real, bounded, band-limited functions on $\mathbb{R}^{n}$ are dual functions in the QPM defined by the Fock embedding.

B Notational conventions

The “Dirac notation” for linear algebra on a Hilbert space is quite useful for the mathematical manipulations that arise in quantum mechanics. Vectors in $\mathcal{H}$ are represented as kets, notated $|\cdot\rangle$. The label inside the ket could represent an arbitrary vector $|\psi\rangle$, an index like $|3\rangle$ or $|n\rangle$, a coherent-state label $|z\rangle$ or $|\alpha\rangle$, or even an experimental outcome like $|\text{up}\rangle$ or $|\text{alive}\rangle$. For instance, instead of denoting an orthonormal basis as $e_{n}$ for indices $n\in\mathbb{N}$, one might simply write $|n\rangle$.

Dual vectors, the linear functionals in the dual space $\mathcal{H}^{*}$, are represented by bras, notated $\langle\cdot|$. Every ket has a dual bra, and vice versa, so we can write $\langle\psi|$, $\langle 3|$, $\langle n|$, $\langle z|$, $\langle\alpha|$, $\langle\text{up}|$ or $\langle\text{alive}|$. The dual relationship between vectors and linear functionals is emphasized in this notation. When a vector (ket) is multiplied by a scalar, the dual of $\lambda|\psi\rangle$ is its Hermitian adjoint $\overline{\lambda}\langle\psi|$. Scalar multiplication with bras or kets is commutative: $\overline{\lambda}\langle\psi|=\langle\psi|\overline{\lambda}$.

The inner product between two vectors $(\phi,\psi)$ is written as the linear functional (bra) $\langle\phi|$ acting on the vector (ket) $|\psi\rangle$ giving the bracket (bra–ket) $\langle\phi|\psi\rangle$. The inner-product notation should be consistent with scalar multiplication of the kets, i.e. $\langle\phi|\lambda|\psi\rangle=\lambda\langle\phi|\psi\rangle$ (since $\lambda$ multiplies the object to its right). This requires the inner product to be linear in the second argument and anti-linear in the first. That is, we would like the inner product notation to mirror the bra–ket notation:

[TABLE]

This is different from the usual math convention that places the dual object in the second slot, giving an inner product that is linear in the first argument. If you have ever wondered why physicists adopt the opposite convention, it is precisely to retain consistency with Dirac notation in this fashion. In practice, writing inner products as bra–kets avoids any ambiguity about which slot is linear: scalars written next to any bra or ket simply follow the usual (commutative and associative) rules for scalar multiplication.

Linear operators are often written with “hats” like $\hat{A}$ or $\hat{n}$. They act on vectors to the right: $\hat{A}|\psi\rangle$ represents the action of $\hat{A}$ on the vector $|\psi\rangle$. Of course, $\hat{A}|\psi\rangle$ is also a vector, so one can act on it using a linear functional $\langle\phi|$ to give $\langle\phi|\hat{A}|\psi\rangle$. A virtue of this notation is that one can equally view this expression as the linear functional $\langle\phi|\hat{A}$ acting on the vector $|\psi\rangle$. The linear functional $\langle\phi|\hat{A}$ is the bra corresponding to the ket $\hat{A}^{\dagger}|\phi\rangle$, where $\hat{A}^{\dagger}$ represents the operator adjoint (Hermitian conjugate) of $\hat{A}$. Scalar multiplication is consistent with the action of scalar multiples of the identity operator, e.g. $\langle\phi|\lambda|\psi\rangle=\langle\phi|(\lambda\hat{I})|\psi\rangle=\lambda\langle\phi|\psi\rangle$.

Dirac notation is especially useful for representing rank-1 projection operators, which represent pure states in quantum mechanics. Given a normalized vector $|\psi\rangle$, the projection operator $\hat{P}_{\psi}$ onto its one-dimensional subspace is simply $|\psi\rangle\langle\psi|$. This dyadic expression can be inserted anywhere you would insert the projector $\hat{P}_{\psi}$: for instance, the projection of a vector $|\phi\rangle$ onto this subspace is simply $\hat{P}_{\psi}|\phi\rangle=|\psi\rangle\langle\psi|\phi\rangle$. This notation eliminates any confusion about whether the projection of a vector $v$ along $w$ is given by $(v,w)w$ or $(w,v)w$.

This notation also shines in representing operators as sums over these dyads. For instance, given an orthonormal basis $|n\rangle$ for a separable Hilbert space, the identity operator can be written as:

[TABLE]

Likewise, a compact self-adjoint operator $\hat{T}$ with eigenvectors $|\phi_{n}\rangle$ and real eigenvalues $\lambda_{n}$ can be written in spectral form as:

[TABLE]

These dyadic sums yield a well-defined operator $\hat{T}$ as long as they converge between any pair of states $\langle\phi|\hat{T}|\psi\rangle$ (that is, in the weak operator topology) and are bounded as a sesquilinear form, i.e. $\langle\phi|\hat{T}|\psi\rangle\leq M\|\phi\|\|\psi\|$ for all $|\phi\rangle$, $|\psi\rangle$, and some finite $M$.

Given a trace-class operator $\hat{S}$, the trace can be written using Dirac notation as:

[TABLE]

for any orthonormal basis $|n\rangle$. The trace of a dyad $|\phi\rangle\langle\psi|$ is simply $\mathrm{tr}(|\phi\rangle\langle\psi|)=\langle\psi|\phi\rangle$, as can be derived using Dirac notation:

[TABLE]

where we recognize the expansion of the identity operator as a sum over orthonormal projectors; this is a lovely example of the utility of Dirac notation. We use this trace identity extensively. For instance, in defining the dual functions for a QPM, given an embedding that maps $x$ to the projector $\hat{\rho}_{x}=|x\rangle\langle x|$, and a bounded operator $\hat{F}$, the dual function is $F(x)=\mathrm{tr}(\hat{\rho}_{x}\hat{F})=\mathrm{tr}(|x\rangle\langle x|\hat{F})=\langle x|\hat{F}|x\rangle$. In quantum mechanics, we would describe this as the expected value of the operator (observable) $\hat{F}$ in the pure state $|x\rangle$.

Indeed, the QPM embedding is at its heart an embedding of probability measures as expected values of rank-1 projection operators. For instance, the standard Fock embedding over $\mathbb{C}$ using pure coherent states $\hat{\rho}_{z}=|z\rangle\langle z|$ is given by the barycenter map:

[TABLE]

which is well-defined as a strong (Bochner) integral on the Banach space $\mathcal{B}_{1}(\mathcal{H})$ of trace-class operators since the map $z\to\hat{\rho}_{z}$ is continuous from $\mathbb{C}$ to the separable space $\mathcal{B}_{1}(\mathcal{H})$, and the trace norm of the integrand is $\|\hat{\rho}_{z}\|_{1}=1$.

If we have a dual function $F(z)$ represented by a bounded operator $\hat{F}$, the expected value of $F(z)=\langle z|\hat{F}|z\rangle$ is given by:

[TABLE]

Let’s work this out carefully as an exercise in Dirac notation. The Banach dual of $\mathcal{B}_{1}(\mathcal{H})$ is the space of bounded operators, $\mathcal{B}(\mathcal{H})$. The dual pairing is given by the trace map, $(\hat{F},\hat{\mu})=\mathrm{tr}(\hat{F}\hat{\mu})$. Since the definition of $\hat{\mu}$ is given by a Bochner integral, it is also valid as a Pettis (weak) integral, which means that the dual pairing of $\hat{\mu}\in\mathcal{B}_{1}(\mathcal{H})$ with $\hat{F}\in\mathcal{B}(\mathcal{H})$ via the trace map can be brought inside the integral, i.e.:

[TABLE]

Finally, just as the expansion of the identity operator in terms of an orthonormal basis shows the strength of Dirac notation, we use a similar expansion for the identity operator in terms of coherent states on $\mathbb{C}^{n}$ (eq. 4 in Lemma 3):

[TABLE]

Here we are careful not to write this as an operator equality; it holds only in the sense of the weak operator topology when “sandwiched” between any pair of states $\langle\phi|\cdot|\psi\rangle$, i.e.:

[TABLE]

by virtue of the Segal-Bargmann-Fock correspondence. This is in contrast to the infinite sums found in operator expansions in terms of rank-1 dyads, which represent a convergent series of operators in various topologies (including the weak operator topology). Here the integral representation of the identity in terms of coherent states does not converge as an operator-valued Dunford or Pettis integral. Thus it is only well-defined when inserted between a pair of states—and this is precisely how we use it in our proofs.

C Replacing MMD with a QPM

Typical applications of MMD involve choosing $n$ i.i.d. samples from a probability measure $\mu$, and $m$ samples from another measure $\nu$, and calculating the MMD between the resulting empirical measures. More generally, let

[TABLE]

be probability measures with finite support $\{x_{k}\}_{k=1}^{n}$ and $\{y_{\ell}\}_{\ell=1}^{m}$. We wish to calculate the MMD between these discrete measures.

Let $\phi:X\to\mathcal{H}$ be the feature map of a characteristic kernel $k^{\prime}$, written $\phi(x)=|x\rangle$ with inner product $\langle x|y\rangle=k^{\prime}(x,y)$. We assume without loss of generality that the kernel is normalized, $\langle x|x\rangle=k^{\prime}(x,x)=1$; if not, multiply $k^{\prime}(x,y)$ by $(k^{\prime}(x,x)k^{\prime}(y,y))^{-1/2}$. Embed the measures by their barycenters,

[TABLE]

so MMD is the Hilbert space distance between these vectors:

[TABLE]

(If $(\mathbb{P},\mathbb{Q})$ are sampled i.i.d. from $(\mu,\nu)$, this is a biased estimate of $\text{MMD}(\mu,\nu)$, i.e. V-statistic.) Writing the difference vector $|D\rangle=|\mathbb{P}\rangle-|\mathbb{Q}\rangle$ as:

[TABLE]

where we define $N=n+m$ and:

[TABLE]

Now MMD is the norm of this difference vector, which we can expand using Dirac notation:

[TABLE]

Calculating this MMD requires evaluating the kernel between every pair of points in the joint support of the measures $\mathbb{P}$, $\mathbb{Q}$, for $O(N^{2})$ kernel evaluations. The matrix $G_{ij}=\langle x_{i}|x_{j}\rangle$ is simply the Gram matrix for those embedded vectors.

Now suppose that we want the norm-squared kernel instead, i.e., we seek the MMD for the kernel $k(x,y)=|k^{\prime}(x,y)|^{2}$. With a different kernel, we obtain a different RKHS. As discussed in the proof of Theorem 2, SI Appendix, A.2, the RKHS for the norm-squared kernel is the Hilbert-Schmidt space $\mathcal{B}_{2}(\mathcal{H})$, where we define the inner product between two Hilbert-Schmidt operators as $(\hat{A},\hat{B})=\mathrm{tr}(\hat{A}^{\dagger}\hat{B})$. In this space, the kernel evaluates to:

[TABLE]

with the self-adjoint pure states (rank-1 projection operators) defined as $\hat{\rho}_{x}=|x\rangle\langle x|$ and $\hat{\rho}_{y}=|y\rangle\langle y|$.

In this space, instead of embedding probability measures as barycenters of embedded vectors in $\mathcal{H}$, we embed as barycenters of embedded operators (pure states) in $\mathcal{B}_{2}(\mathcal{H})$:

[TABLE]

and define the MMD using the Hilbert-Schmidt norm:

[TABLE]

Following eq. 12, we define the difference operator $\hat{D}$ as:

[TABLE]

using the same notation as eq. 13-14; note $\mathrm{tr}(\hat{D})=\sum_{i}c_{i}=0$. An MMD calculation using the Hilbert-Schmidt inner product yields:

[TABLE]

which is identical to the MMD calculated using the norm-squared kernel $k(x,y)=|k^{\prime}(x,y)|^{2}$. Equivalently, if the eigenvalues of the self-adjoint operator $\hat{D}$ are given by $\lambda_{i}$, the MMD is given by:

[TABLE]

We calculate the QPM using the same embedding and the same difference operator $\hat{D}$ (eq. 15) but with trace distance instead of Hilbert-Schmidt distance. With the conventional normalization factor of $1/2$ (giving distances in [0,1]), we have:

[TABLE]

where $|\hat{D}|\equiv(\hat{D}^{2})^{1/2}$ represents the (operator) square root of the nonnegative operator $\hat{D}^{2}$. Here we see the central computational challenge in using QPM. Unlike MMD, which can be calculated using $O(N^{2})$ kernel evaluations, the operator square root required for QPM needs eigenvalues of the operator $\hat{D}$, with typical complexity $O(N^{3})$ for an $N\times N$ matrix. Recall that the eigenvalues of $|\hat{D}|$ are the absolute values of the eigenvalues of $\hat{D}$, and the trace of an operator is the sum of its eigenvalues. Thus, if the eigenvalues of $\hat{D}$ are given by $\lambda_{i}$, the QPM is:

[TABLE]

Now we need the eigenvalues of the difference operator $\hat{D}=\sum_{i=1}^{N}c_{i}|x_{i}\rangle\langle x_{i}|$. This is not directly amenable to matrix calculation, as the embedded vectors $|x_{i}\rangle$ do not define an orthonormal basis, and the only quantities we can easily calculate are kernel evaluations $k^{\prime}(x,y)=\langle x|y\rangle$. However, we can use a bit of linear algebra to help. Choose a set of $N$ orthonormal basis vectors $|e_{i}\rangle\in\mathcal{H}$, whose linear span includes the span of the $|x_{i}\rangle$, and define the operators $\hat{V}=\sum_{i=1}^{N}|x_{i}\rangle\langle e_{i}|$ and $\hat{C}=\sum_{j=1}^{N}c_{j}|e_{j}\rangle\langle e_{j}|$ so that the difference operator can be expressed as:

[TABLE]

using $\langle e_{j}|e_{i}\rangle=\delta_{ij}$. As operators on $\mathrm{span}\{|e_{i}\rangle\}$, their eigenvalues are preserved under cyclic permutations, so the eigenvalues of $\hat{D}=\hat{V}\hat{C}\hat{V}^{\dagger}$ are the same as those of $\hat{C}\hat{V}^{\dagger}\hat{V}$. Now observe that the operator $\hat{G}=\hat{V}^{\dagger}\hat{V}$ defines the Gram matrix in the $|e_{i}\rangle$ basis:

[TABLE]

i.e. $G_{ij}=\langle e_{i}|\hat{G}|e_{j}\rangle=\langle x_{i}|x_{j}\rangle$. So we can equivalently find the eigenvalues of the operator $\hat{C}\hat{G}$, whose matrix elements are given by the matrix product $M=CG$, i.e. $M_{ij}=c_{i}G_{ij}=c_{i}\langle x_{i}|x_{j}\rangle$. Although this is a non-Hermitian matrix, its eigenvalues will be real since they are identical to those of the Hermitian operator $\hat{D}$. They can be calculated using a general eigenvalue solver; any imaginary parts that arise due to rounding errors should be discarded to preserve the essential relationships $\mathrm{tr}(\hat{D})=\sum_{i}\lambda_{i}=\sum_{i}c_{i}=0$.

Nonetheless, calculating eigenvalues of non-Hermitian matrices can be computationally costly and numerically unstable. For better numerics, we factor the Gram matrix as $G=HH^{\dagger}$, e.g. using Cholesky factorization or the matrix square root. Using invariance under cyclic permutation, the eigenvalues of $CG=CHH^{\dagger}$ are the same as those of $H^{\dagger}CH$, which is a Hermitian matrix with numerically stable (real) eigenvalues. These are the eigenvalues we want—they are the eigenvalues of the difference operator $\hat{D}$.

We now have a simple drop-in “recipe” for replacing MMD with QPM in any calculation that uses finite discrete measures (e.g. empirical measures). MMD already requires the Gram matrix $G_{ij}=\langle x_{i}|x_{j}\rangle=k^{\prime}(x_{i},x_{j})$. Instead of calculating MMD using the quadratic form $c_{i}c_{j}G_{ij}$, factor the Gram matrix into $G=HH^{\dagger}$, define the coefficient matrix $C=\mathrm{diag}(c_{i})$, and find the eigenvalues $\lambda_{i}$ of the Hermitian matrix $H^{\dagger}CH$. The QPM will be $\tfrac{1}{2}\sum_{i=1}^{N}|\lambda_{i}|$.

There is one subtle point, however. The original MMD used the kernel $k^{\prime}(x,y)=\langle x|y\rangle$, while the QPM calculation effectively uses the norm-squared kernel $k(x,y)=|k^{\prime}(x,y)|^{2}=\langle y|x\rangle\langle x|y\rangle$. The difference is clear when we consider the distance between two Dirac measures (point masses), $\mathbb{P}=\delta_{x}$ and $\mathbb{Q}=\delta_{y}$, embedded as vectors in a Hilbert space $\mathcal{H}$: $\delta_{x}\to|x\rangle$ and $\delta_{y}\to|y\rangle$. Assuming (as we have throughout) that all embedded vectors are normalized, $\langle x|x\rangle=\langle y|y\rangle=1$, the MMD in $\mathcal{H}$ is given by:

[TABLE]

In contrast, we obtain the norm-squared kernel if we embed the Dirac measures as operators in the Hilbert-Schmidt space $\mathcal{B}_{2}(\mathcal{H})$: $\delta_{x}\to\hat{\delta}_{x}=|x\rangle\langle x|$ and $\delta_{y}\to\hat{\delta}_{y}=|y\rangle\langle y|$. The MMD is given by the norm of the difference operator $\hat{D}=\hat{\delta}_{x}-\hat{\delta}_{y}=|x\rangle\langle x|-|y\rangle\langle y|$:

[TABLE]

The QPM for these measures can be calculated by diagonalizing the difference operator $\hat{D}=|x\rangle\langle x|-|y\rangle\langle y|$. The eigenvalues of this operator are: $\lambda_{\pm}=\pm\sqrt{1-|\langle x|y\rangle|^{2}}$, so we calculate:

[TABLE]

Up to a factor of $\sqrt{2}$, therefore, the QPM between point masses is identical to the MMD calculated using the Hilbert-Schmidt metric, i.e. with the norm-squared kernel. Conversely, MMD using a kernel $K$ will give the same point-mass metric as QPM using the square-root kernel $\sqrt{K}$. This is an important detail when replacing MMD with QPM.

Thus, given a calculation that uses MMD, one should first ask if the MMD kernel can be written as the norm-square of another kernel. For many common kernels, including the Gaussian, Laplacian, and inverse-multiquadric kernels, the square root of each of these kernels is also a strictly positive definite kernel. In that case, the best recipe for replacing it with QPM is as follows:

We are given an MMD calculation between two probability measures with finite support, i.e.,

[TABLE] 2. 2.

For convenience, use a single set of indices from $1$ to $N=n+m$ with coefficients and points of support defined as:

[TABLE]

[TABLE] 3. 3.

Where MMD was calculated from the Gram matrix $\langle x_{i}|x_{j}\rangle=K(x_{i},x_{j})$, calculate a new Gram matrix using the square root of the kernel, i.e. $\langle x_{i}|x_{j}\rangle=\sqrt{K(x_{i},x_{j})}$. One must ensure here that $\sqrt{K}$ is also a valid (positive definite) kernel. 4. 4.

Factor this new Gram matrix $G=HH^{\dagger}$ using Cholesky factorization. Define $M=H^{\dagger}CH$ where $C=\mathrm{diag}(c_{i})$, and find the eigenvalues $\lambda_{i}$ of this Hermitian $M$. If some points are coincident, $x_{i}=x_{j}$ for $i\neq j$, Cholesky will fail; one could eliminate these coincidences to restore a positive-definite Gram matrix or use a general eigensolver and extract the (real) eigenvalues $\lambda_{i}$ of the non-Hermitian matrix $CG$. 5. 5.

Then $\text{QPM}(\mathbb{P},\mathbb{Q})=\tfrac{1}{2}\sum_{i=1}^{N}|\lambda_{i}|$ with conventional normalization. If isometry between the MMD and QPM metrics for point masses is desired instead, multiply the QPM by $\sqrt{2}$.

In case the kernel $K$ used in the existing MMD calculation cannot be written as the norm-square of another kernel, there are two reasonable options. One could proceed using the kernel $K$ directly, in which case the QPM base metric effectively uses $|K|^{2}$ as its kernel. This could have undesirable effects: for instance, the effective bandwidth of $|K|^{2}$ will be smaller than that of $K$.

Alternately, one could find a kernel whose norm-square closely approximates the original kernel over the relevant range. We use the second approach in the GMMN example below: the original kernel was defined as a sum of Gaussians with a range of widths, giving a “fat-tailed” kernel that was effective in that particular application. This sum of Gaussians cannot be written as the square of any kernel, but the fat-tailed kernel can be matched with an inverse-multiquadric kernel, which does have a valid positive-definite kernel square root. In practice, this second approach will tend to follow the original geometry more closely.

D QPMs in Generative Moment Matching Networks (GMMNs)

As an example of the “drop-in” replacement of MMD with QPM, we turn to a classic generative modeling task that relies on MMD to measure the distance between probability distributions: the Generative Moment Matching Network (GMMN) [10]. We reproduce this result and then replace the MMD metric with QPM, keeping all else constant (architecture, optimizer, training schedule, etc.).

The original GMMN implementation used a 5-layer fully-connected feed-forward network to map a random latent vector $z$ into a data space, such as the space of MNIST images [11]. This transforms a fixed probability distribution $P_{z}$ on the latent space into a distribution on the data space, and training seeks to match that distribution to that of the target dataset, using the MMD metric as its objective. This is a conceptually simple and appealing idea: the distribution of generated images should match as closely as possible to the actual distribution. In the specific example from [10], they begin with a uniform distribution on $[-1,1]^{10}$ in a 10-dimensional latent space and pass it through 5 layers using ReLU activation with 64, 256, 256, and 1024 hidden dimensions respectively, ending with sigmoid activation on 784 dimensions to match the $28\times 28$ grayscale MNIST images with pixel values in $[0,1]$. The net was trained by choosing random batches of 1000 images from the dataset, generating 1000 images from the network, and backpropagating to minimize MMD. The code is available from the original author at https://github.com/yujiali/gmmn.

We first re-implemented this network architecture using PyTorch and reproduced the results from the paper. Then, from this baseline, we swapped QPM in place of MMD as the objective function to run a side-by-side comparison of these two probability metrics. This presented two key challenges. First, the kernel used for MMD in the original paper was an equal mixture of 6 Gaussians with a range of values for $\sigma^{2}=\{2,5,10,20,40,80\}$. (Note that the original author’s code uses the variable name sigma for these values but careful examination shows that they actually represent $\sigma^{2}$ in the usual parameterization of the Gaussian kernel.) As discussed at the end of SI Appendix, C, although the sum of Gaussians is a valid kernel, its square root is not. Thus, to preserve the geometry in the embedding space, we created a kernel that closely matched the original mixture of Gaussians over the data space. We used a generalized inverse-multiquadric kernel with two parameters: a length scale $\lambda$ and exponent $\alpha$, yielding:

[TABLE]

With the best-fit parameters, this kernel yielded an excellent match to the “fat-tailed” mixture of Gaussians and reproduced the original results nicely. Now we could use the (valid) square root of this kernel as the basis for the QPM, ensuring that the fundamental embedding geometry is preserved.

The second challenge was related to the MNIST dataset, which has many pairs of images that are very similar to one another, with Euclidean distance $\|\mathbf{x}-\mathbf{y}\|$ near zero, or, equivalently, a kernel inner product $k(\mathbf{x},\mathbf{y})$ very close to 1 (our kernels are normalized to $k(\mathbf{x},\mathbf{x})=1$). As a result, even if the Gram matrix for 1000 samples from the dataset is not truly singular, rounding error can cause Cholesky factorization to fail. We tried two approaches that yielded identical results: (1) numerical regularization by adding a small multiple of the identity ($\epsilon=10^{-4}$) to the Gram matrix, and (2) using a generalized eigensolver directly with the non-Hermitian matrix $CG$ ($C=\mathrm{diag}(c_{i})$; see the end of SI Appendix, C, for notation). The first approach is faster and is probably superior for most applications: the stochastic nature of the sampling and generation process is likely to outweigh any small errors in calculating the “true” eigenvalues (and their gradients). Nonetheless, for a robust implementation of QPM, the generalized eigensolver offers a backup in the event that Cholesky fails.

To generate the images in Figure 1, we trained for 200 epochs on the MNIST training set (50,000 images) using batches of 1000 and the PyTorch AdamW optimizer with a learning rate of 0.005 and default parameters. (See https://doi.org/10.5281/zenodo.17024967 for results and code.) A single training loop trained two models together with the only difference being the objective function (QPM versus MMD). The MMD results are qualitatively similar to those from [10], with many visible artifacts that do not improve beyond 200 epochs. In contrast, the images generated using QPM appear much more similar to the MNIST samples.

After training, we used a kernel two-sample test [2] to assess whether the generated images were similar to images from the dataset. We compared a batch of 1000 real images (from a held-out test set) with a batch of 1000 generated images. The baseline distance between these batches was measured using the same metrics as for training (MMD or QPM). Then, to determine if this distance was statistically significant, we randomly permuted all 2000 images into two new batches of 1000 each. We repeated this shuffling process 1000 times, calculating the distance for each random split. If the real vs. generated images were from the same distribution (null hypothesis), then the baseline distance should be typical of these random distances. This would imply that, based on that metric, the trained network met its goal of matching generated images to the true distribution.

Given $N=1000$ permutations, if $r$ of them exceed the baseline distance, the continuity-corrected $p$-value is $(r+1)/(N+1)$. Repeating each two-sample test for 30 distinct batches and taking the average, we find $p\approx 0.006$ for the MMD images and $p\leq 0.001$ for the QPM images. In both cases, these metrics can detect that the generated images are not drawn from the same distribution as the MNIST dataset—confirming our visual inspection.

For a more challenging test, we trained a GMMN on the CelebA dataset [12] (https://mmlab.ie.cuhk.edu.hk/projects/CelebA.html), pre-processed to yield $64\times 64$ RGB images—a data space of 12,288 dimensions. This dataset contains over 200,000 images with faces in various orientations, lighting conditions, and backgrounds. In place of the MLP used for MNIST, we adopted the convolutional generator DCGAN_G from https://github.com/martinarjovsky/WassersteinGAN. (Note there is nothing “adversarial” about our approach; there is no discriminator network.) The DCGAN_G network starts with a multivariate standard normal distribution on a 100-dimensional latent space and uses 5 convolutional layers with BatchNorm and ReLU after each hidden layer and tanh activation for the final layer (with pixel values scaled to $[-1,1]$). Two models were trained together (training set of 180,000 images, 50 epochs, batches of 1000), comparing QPM with MMD as objective functions. We used a simple Gaussian kernel to avoid the “kernel square root” problem (the square root of a Gaussian kernel is also a valid kernel), and found that the diversity of CelebA images ensured that Cholesky factorization typically succeeded without regularization.

Examining Figure 1, the images generated using these two different metrics appear quite different. Many of the MMD images have bizarre artifacts, while all the QPM images show clear facial structure. More strikingly, the two-sample test for the MMD-generated images, using the same metric, yielded $p\approx 0.23$ (average of 30 runs); this metric could not reject the null hypothesis that the generated images were drawn from the same distribution as CelebA. This result shows the limitations of MMD in high-dimensional spaces: the training has optimized the MMD-distance between the data and the generated images, and statistically the two sets of images are similar based on this metric, but the images are not even remotely similar. In contrast, QPM readily distinguishes its own generated images or the MMD-generated images from the CelebA data, with $p\leq 0.001$ in both cases. We note that the MMD metric became similarly saturated (unlike QPM) for Gaussian kernels with $\sigma$ ranging from 15 up to 80 (the median distance between pairs of image vectors drawn from the dataset); the results in Figure 1 took $\sigma$ as half this median distance.

While training the CelebA generator using QPM alone doesn’t generate photorealistic images, at least this metric can detect that the generated images are different from the dataset (using the two-sample test). In contrast, optimizing for MMD yielded very poor images that were nonetheless viewed by that metric as being statistically indistinguishable from the training data. This simple experiment shows that MMD can be unreliable in high dimensions where QPM is robust. The improved performance of QPM is most likely due to its rich class of dual functions (Theorem 5, SI Appendix, A.5). In contrast, the dual functions for MMD must vanish at infinity, so its witness functions (unit ball in the RKHS) effectively have limited support, even on compact domains. It is likely that these witness functions can only detect differences along a limited number of dimensions, leaving large parts of the high-dimensional space effectively unmonitored.