Interpretable Quantum Photonic Generative Models

Introduction

Many quantitative problems reduce to learning a joint distribution $P(\mathbf{x})$ of $N$ continuous variables from data and drawing new samples from it [1, 2]. Generative models solve this problem by capturing the dependence between variables and producing new samples cheaply once trained, as done by variational autoencoders, generative adversarial networks, normalizing flows, and diffusion models [3, 4, 5, 6, 7, 8]. However, the sampling cost grows with the data dimension and the model size, as inference requires a pass through the network or a sequence of such passes [6, 8]. This motivates the search for physical systems whose natural measurement statistics realize a trainable multivariate distribution, including coherent Ising machines [9] [10], probabilistic bits [11], thermodynamic computers [12], optical neural networks [13, 14], and polariton condensates [15].

Quantum generative models implement physical sampling through projective measurements of a parametrized quantum state [16, 17, 18]. The Born rule $p(m) = \operatorname{Tr}[|m\rangle\langle m|\hat{\rho}_{\boldsymbol{\theta}}]$ turns each measurement outcome $m$ into a sample from a distribution set by the state $\hat{\rho}_{\boldsymbol{\theta}}$ and its circuit parameters $\boldsymbol{\theta}$ [19, 20, 21, 22]. Models differ in the architecture, the circuit parametrization, the measured observable and its support, and the training strategy [23, 24]. Architectures include qubit-based circuits [16, 17, 21, 25], photonic samplers based on boson sampling and Gaussian boson sampling [26, 27, 28, 29], and fermionic linear-optical samplers [30, 31]. Parametrizations include quantum circuit Born machines [16, 32], quantum Boltzmann machines [33, 34], spectral generative models [23, 35, 36, 37], and models with physically parametrized dynamics [38]. Restricted architectures aid classical training while keeping sampling hard, as pursued in “train-on-classical, deploy-on-quantum” approaches [35, 39, 40, 41, 42]. Training strategies rely on tractable losses, including the maximum mean discrepancy (MMD) [16, 24, 43] and sample-based estimators [44].

Photonic systems with continuous-variable (CV) encoding offer a distinct architecture for quantum information processing [45, 46, 47]. Optical modes carrying coherent and squeezed light are manipulated with programmable interferometers [48, 49, 50], and the resulting Gaussian states are described by covariance matrices that are closed-form functions of the circuit parameters, so the correlation structure of the output can be designed. The description leaves the Gaussian family in two cases: when photon-number-resolving detection is used [51, 52], and when a non-Gaussian gate based on a Kerr or other nonlinearity is inserted [53, 54, 55]. The former underlies boson sampling [26, 27, 56, 57, 58], and the latter underpins non-Gaussian photonic computing [59, 60]. Experimentally, CV systems produce large multimode entangled states [61, 62], support multimode homodyne detection at 60 modes [63], and can be used to sample from distributions with over more than 400 modes [64]. From the learning perspective, photonic circuits were proposed for quantum neural networks [65, 66], and were recently used to show advantage in detecting signal from limited quantum data [67] as well as learning relational information [68]. Gaussian boson samplers were applied to optimization and graph problems [69, 70, 71]. For generative modeling, Born machines were built on boson sampling [28, 72, 73], and photonic circuits served as generators and latent samplers in adversarial networks [74, 75]. Closest to our work, the CV Born machine of Ref. [76] reads out homodyne samples of a variational CV circuit and trains it on these samples with a kernel-based loss, with demonstrations on single-mode targets and quantum states of up to three modes. Open questions are whether a structured CV circuit can be trained without sampling, on the covariance of the state alone, such that each optical layer has a stated role in the multivariate correlations of an empirical dataset, and whether this remains trainable at tens of modes.

In this work, we propose a CV photonic generative architecture built from pairwise squeezing, local phase rotations, and a programmable interferometer, read out by homodyne detection of the position quadratures. In our model squeezing creates correlations, phases set their sign and magnitude, and the interferometer redistributes them across mode pairs, such the effect of each layer on the output correlation matrix becomes more explainable. We train the circuit classically on an empirical correlation matrix using only the covariance of the simulated state, and obtain continuous samples by simulated homodyne detection of the trained circuit. To reach beyond the Gaussian family we include a Kerr nonlinearity [53, 54] and represent its weak-strength action as a coherent sum of Gaussian branches, which keeps training within Gaussian phase-space methods. We test the architecture on daily financial returns [77, 78] with up to twenty-mode circuits, reproducing an empirical correlation matrix up to small residuals and recovering the sign of correlations for most of the pairs.

Model

We build a photonic circuit whose homodyne samples follow a multivariate distribution with tunable correlations. We use Gaussian gates with known effects on the position-quadrature correlations, add local Kerr gates to introduce non-Gaussianity, and train the circuit classically to match a target correlation matrix.

Photonic variables and structured ansatz

Each optical mode $j=1,\ldots,N$ carries position and momentum quadratures $\hat{x}_j$ and $\hat{p}_j$ with commutation relations $[\hat{x}_j,\hat{p}_k]=2i\delta_{jk}$. We collect the quadratures into a vector $\hat{\mathbf{R}}=(\hat{x}_1,\hat{p}_1,\ldots,\hat{x}_N,\hat{p}_N)^{\mathsf{T}}$. A Gaussian state is specified by the first moments $\mathbf{d}=\langle\hat{\mathbf{R}}\rangle$ and the $2N\times2N$ covariance matrix

with $\{\cdot,\cdot\}$ the anticommutator, and Gaussian gates act on $\Sigma$ as symplectic maps [46] (App. A). The vacuum then has $\Sigma/2=1$, and Pearson correlations do not depend on this choice of normalization.

The model outputs are obtained by homodyne detection of $\hat{x}_j$ on every mode. The outcome $\mathbf{x}\in\mathbb{R}^{N}$ is distributed according to the position marginal of the Wigner function, which for a Gaussian state is a multivariate normal distribution with covariance given by the $N\times N$ position block $\Sigma_{ij}^{(x)}=\Sigma_{2i-1,2j-1}$. The output correlations are then determined by

Note that our choice of measurement is distinct from boson sampling models, where photon-number-resolving

detection of the same state returns integer counts with probabilities governed by hafnians [27].

We prepare a CV photonic state by applying an ansatz $U(\boldsymbol{\theta})$ to the vacuum $|\mathbf{0}\rangle=|0\rangle^{\otimes N}$, where

$(j,j')$ runs over the pairs $(1,2)$, $(3,4),\ldots$ and, for odd $N$, the last mode receives a single-mode squeezer $S(r_s)$ (Figure 1). The trainable parameters $\boldsymbol{\theta}$ are the squeezing amplitudes, the local phases $\phi_i$, and the angles $\boldsymbol{\vartheta}$ of the $N(N-1)/2$ beam splitters of the Clements mesh [49]. The internal phases of the beam splitters are fixed at zero, and the mesh represents a restricted real-orthogonal mixing layer. The homodyne measurement follows the Born rule and reads $p(\mathbf{x})=|\langle\mathbf{x}|U(\boldsymbol{\theta})|\mathbf{0}\rangle|^2$. In each run we then get samples $\mathbf{x}=(x_1,x_2,\ldots,x_N)\sim p(\mathbf{x})$.

The three Gaussian layers shape the correlations in distinct ways. Two-mode squeezing correlates the position quadratures of a vacuum pair, with the magnitude set by the squeezing amplitude $r$. A local phase rotation changes which quadrature is measured, which tunes the sign and magnitude of the correlation at fixed $r$. The Clements mesh then redistributes these correlations across all modes (see III A). While the ansatz is simple, it spans a variety of correlated distributions and returns one $N$-dimensional sample per shot, with a potentially short time per sample.

Once the model is defined, we test it in silico with established simulators of photonic systems, Strawberry Fields [79] and Piquasso [80]. The simulation workflow splits into a Gaussian backend based on symplectic propagation, which is fast, and a Fock backend that works in a truncated Fock basis, which is numerically heavy and depends on the cutoff. We use Piquasso for performance, with the Gaussian backend for training and the Fock backend for testing the effects of the Kerr nonlinearity.

Image not extracted

Kerr nonlinearity and Gaussian-branch representation

Gaussian circuits with Gaussian measurements are efficiently simulable classically [81, 82], and the simulation cost increases once non-Gaussian gates, such as a Kerr gate, enter the preparation circuit. These gates lead to output distributions outside the normal family and increase the expressivity of the generative model. We use the single-mode Kerr gate $K(\kappa) = \exp(i\kappa\hat{n}^{2})$ with $\hat{n} = \hat{a}^{\dagger}\hat{a}$ [83, 84], whose quartic generator imprints an intensity-dependent phase and distorts the Wigner function.

To keep the simulation within Gaussian methods we represent the Kerr gate as a linear combination of Gaussian operations. This follows from the Hubbard–Stratonovich identity [85, 86], which writes $\exp(i\kappa\hat{n}^{2})$ as an integral over $\exp(ix\hat{n})$ with a Gaussian kernel,

for $\kappa> 0$ with the principal branch of the square root. A quadrature rule for the oscillatory integrand [87] gives a finite sum (App. B), and the action on a Gaussian state $|G\rangle$ becomes

where the complex weights $c_m$ and real nodes $x_m$ come from the discretization of Eq. (4). Each branch $|G_m\rangle$ is Gaussian and propagates through the remaining layers within the Gaussian formalism. The branches add as amplitudes, so expectation values contain interference terms between every pair of branches,

where the denominator normalizes the truncated sum. Overlaps and matrix elements of quadrature monomials between Gaussian states are available in closed form, so the correlators needed for Eq. (2) follow by choosing $\hat{O}=\hat{x}_i\hat{x}_j$ and $\hat{O}=\hat{x}_i$ and subtracting the means. The accuracy of the representation is controlled by the approximation parameters, in particular the number of nodes $M$, and the number of branches multiplies across modes. More details on this decomposition are provided in App. B.

In the following, we work in the regime of weak Kerr nonlinearity, with $\kappa= 10^{-5}$ in the five-mode circuit and vanishing nonlinearity in the twenty-mode circuit, where the representation of Eq. (6) requires a few thousand nodes per gate (App. B).

Training, sampling, and evaluation

We train the model on an empirical dataset $\mathcal{D}_{\mathrm{target}}$ by fixing its output correlations to the target computed from data. Inference then returns real-valued multivariate samples with these correlations. From $T$ joint observations $\{\mathbf{x}^{(t)}\}_{t=1}^{T}$ of $N$ variables we compute the sample Pearson correlation matrix $C_{\mathrm{target}}$ and assign each variable to one mode. For a parameter vector $\boldsymbol{\theta}$ we evaluate the output state, extract $\Sigma^{(x)}(\boldsymbol{\theta})$, and form $C_{\mathrm{gen}}(\boldsymbol{\theta})$ by Eq. (2). The loss is a distance between the model and empirical matrices, measured as a maximum mean discrepancy (MMD) [43]. With $r_i$ the $i$-th row of a correlation matrix and the kernel $k(\mathbf{u},\mathbf{v})=\exp(-\gamma\|\mathbf{u}-\mathbf{v}\|^{2})$ with $\gamma=1$, the MMD loss reads

We minimize Eq. (7) with the bound-constrained quasi-Newton method L-BFGS-B [88] under the box constraints of App. C. Training uses only the covariance of the simulated state and requires no sampling, so it corresponds to moment matching. Minimizing the distance between correlation matrices does not guarantee a match of the full model and empirical distributions, e.g. in total variation distance, but it reproduces the correlations and serves as a pre-training route. The marginals are adjusted in a post-processing stage applied to the generated samples.

Post-processing with affine maps. After training we draw position-quadrature homodyne samples $\mathbf{x}:=\mathbf{x}_{\mathrm{raw}}$ from the optimized circuit. The correlation loss is invariant under a rescaling of each mode, so the raw samples carry the learned correlations but arbitrary scales. We apply a per-mode affine map

with $a_i>0$ and $b_i$ chosen to match the mean and variance of each generated marginal to the target. For positive $a_i$ the Pearson correlations are unchanged. The map adjusts location and scale and cannot alter the shape of a marginal.

We evaluate the models with quantities distinct from the training loss. For correlation agreement we use the mean squared error $\mathrm{MSE}_{\mathrm{full}}=N^{-2}\sum_{i,j}[(C_{\mathrm{target}})_{ij}-$ $(C_{\mathrm{gen}})_{ij}^{2}$ over all entries, the entrywise residuals $\Delta C_{ij} = (C_{\mathrm{target}})_{ij} - (C_{\mathrm{gen}})_{ij}$, and the sign-agreement rate, the fraction of the $N(N - 1)/2$ pairs $i < j$ for which $C_{\mathrm{gen}}$ and $C_{\mathrm{target}}$ have the same sign. For the marginals we use the total variation distance $\mathrm{TVD} = \frac{1}{2}\sum_b |p_b - q_b|$ between histograms of generated and target samples and the two-sample Kolmogorov–Smirnov (KS) statistic $D$ with its $p$-value. Estimator settings are given in App. C.

Results

We now test the approach, which requires continuous multivariate datasets with a nontrivial correlation structure. We choose financial datasets, which exhibit positive and negative correlations, multivariate trends, and heavy tails, while retaining a significant normal component [1, 77]. Specifically, we collect daily returns of several assets over a fixed observation window. Daily adjusted closing prices $P_t$ are retrieved with yfinance [89], converted to log returns $r_t = \log P_t - \log P_{t-1}$, and standardized to zero mean and unit variance per asset. The target $C_{\mathrm{target}}$ is the sample Pearson correlation matrix over the window, and we study joint return statistics at a fixed window without any forecasting. Two ensembles are used: a five-asset set consisting of the equities AAPL and MSFT, the gold fund GLD, the long-maturity Treasury fund TLT, and the bitcoin-dollar pair BTC-USD over $T = 251$ daily returns, and a set of twenty S&P 500 equities over $T = 48$ daily returns. We present more details on acquisition and observation windows in App. C.

Two-mode mechanism

Before training in the general $N$-mode case, we examine the properties of the model on a two-mode example. In Figure 2 we show the position-quadrature covariance and the Pearson correlation of a two-mode squeezed vacuum for $r \in[0, 2]$. The covariance grows exponentially with $r$, while the correlation rises steeply at small $r$, exceeds 0.9 for $r \gtrsim0.8$, and saturates toward unity without reaching it at finite squeezing. Squeezing therefore creates correlation between a mode pair and sets its magnitude. With the squeezing fixed, a phase rotation applied to one mode of the pair and swept through one period modulates the measured correlation as a cosine of the applied phase, from its maximal positive value through zero at $\pi/2$ to its maximal negative value at $\pi$ (App. § 3.8). A local phase thus sets the sign of a measured correlation and tunes its magnitude continuously. The mesh then mixes the modes and redistributes these correlations across all pairs.

Using this interpretation, the model becomes a design tool for distributions in which positive and negative correlations are defining properties. For instance, an inspection of an empirical correlation matrix shows which pairs are strongly correlated and with which sign, which suggests a suitable initialization. Based on this, we select and pair modes for two-mode squeezing, choose the initial squeezing amplitudes, and decide which local phases start near 0 or $\pi$. Such a physically informed initialization can precede numerical optimization or serve as a pre-training stage before the mesh parameters are refined.

Five-mode demonstration

We proceed to model a multivariate distribution from the financial domain, instantiating the ansatz with $N = 5$. Two-mode squeezers are applied on pairs $(1, 2)$ and $(3, 4)$, and a single-mode squeezer applied on mode 5. Local phase rotations and single-mode Kerr gates with $\kappa= 10^{-5}$ act on every mode, followed by a $C$lements mesh of ten beam splitters. The circuit has 18 trainable parameters, corresponding to five phases, three squeezing amplitudes, and ten beam-splitter angles. We minimize Eq. (7) with L-BFGS-B from one fixed initialization, which is the run shown in Figure 3 and Figure 4, as well as from ten random perturbations around that initialization (App. C).

We compare the target and generated matrices in Figure 3. The optimized circuit reproduces the sign and magnitude of all ten pairwise correlations, including the negative correlations of TLT with the two equities and the weak correlations involving GLD. The residuals are bounded within $4.1 \times10^{-5}$, with $\mathrm{MSE}_{\mathrm{full}} = 3.7 \times10^{-10}$. All ten perturbed initializations reach $\mathrm{MSE}_{\mathrm{full}}$ below $3.1 \times10^{-8}$ within 4600 function evaluations (App. § 3.8).

Figure 4 shows the post-affine homodyne samples beyond their second moments. The raw generated samples are unimodal and centered near zero, with widths that differ across assets because the loss does not constrain the scale of each mode (App. § 3.8, Table II). The affine map brings every generated marginal onto the target mean and variance and visibly improves the overlap in the central region. We note that the targets have sharper peaks and heavier tails than a normal distribution, while the marginals of the trained circuit, which is Gaussian up to the weak Kerr distortion, remain close to normal under any affine map. The mismatch is largest for the most volatile asset, BTC-USD. Quantitatively (App. § 3.8, Table II), the affine map reduces the asset-averaged TVD from 0.370 to 0.245 and the averaged KS statistic from 0.190 to 0.121, with BTC-USD retaining the largest TVD of 0.330. For GLD, whose raw marginal already lies close to the target, the affine map raises the TVD from 0.165 to 0.225 and the KS statistic from 0.066 to 0.087, and for TLT it leaves the TVD essentially unchanged (0.273 and 0.271). The correlation structure is thus learned to high accuracy. Learning the tails of the marginals and higher-order correlations requires further increasing the expressivity of the generative model with non-Gaussian gates, at the expense of more complex simulation.

Twenty-mode demonstration

We extend our tests to $N = 20$ multivariate distributions without structural changes. The model includes ten two-mode squeezers on consecutive pairs, a local phase on every mode, and a mesh of 190 beam splitters. With zero Kerr nonlinearity this is a Gaussian circuit with 220 trainable parameters (20 phases, 10 squeezing amplitudes, 190 beam-splitter angles), trained on the empirical correlation matrix of twenty equities from a fixed initialization, which is the run shown in Fig. 5, and from its five random perturbations.

In Fig. 5 we compare the matrices for one optimized circuit. The generated matrix reproduces the block of strongly correlated equities around modes 8 to 12 and the weaker, partly negative correlations elsewhere. Averaged over seeds, the sign of 93.8% ± 0.9% of the 190 pairs is recovered (standard deviation across seeds, individual runs from 92.6% to 94.7%), 86% of the pairs have residuals below 0.10, and no residual exceeds 0.5 (App. § 3.8, Table III). Since equity correlations are predominantly positive, a constant positive sign would also score well above 50%, so we read the sign agreement together with the residual distribution. All five runs reach a comparable plateau of $\mathrm{MSE}_{\mathrm{full}}$ (App. § 3.8, Figure 9). The target is estimated from 48 observations of 20 variables and is a noisy estimate of the dependence between the

caption: Correlation learning with a twenty-mode Gaussian circuit (220 trainable parameters) on twenty S&P 500 equities. (a) Empirical Pearson correlation matrix $C_{\mathrm{target}}$ of the standardized daily log returns over $T = 48$ trading days. (b) Correlation matrix $C_{\mathrm{gen}(\theta^{*})}$ of the optimized circuit for the run started from the fixed initialization of App. C. (c) Entrywise residuals $\Delta C_{ij} = (C_{\mathrm{target}})_{ij} - (C_{\mathrm{gen}})_{ij}$ on a color scale of $\pm0.36$. The block of strongly correlated equities around assets 8 to 12 and the weaker, partly negative correlations elsewhere are reproduced. Aggregate performance over five random initializations (sign agreement and residual distribution) is reported in the text and in App. § 3.8, Table III.

:::

assets, and the model is trained to reproduce this empirical matrix. The experiment shows that the ansatz and the correlation loss remain trainable with 220 parameters and 190 target entries at this size.

Discussion

We first place the model relative to boson sampling. Gaussian boson sampling [27, 56], its Fock-state and non-linear variants [26, 90], and the advantage experiments built on it [57, 58] use squeezers and an interferometer (as in Figure 1) but read the output with photon-number-resolving detectors. Their hardness is based on the readout. We read the same states with homodyne detection, which gives a Gaussian output whose covariance we compute and sample classically in polynomial time [81, 82]. This makes the training classically tractable, while hardness results for homodyne-sampled circuits require non-Gaussian elements of finite strength [55, 91, 92]. To gain a computational advantage one therefore needs models with larger $\kappa$, at the expense of simulation complexity.

Next, we discuss experimental benefits and limitations of CV optics. Continuous-variable systems offer a versatile optical platform [47, 93], but photon loss is the main limitation for CV-based processing [46, 94]. Consider a loss channel that admixes vacuum into each mode with transmission $\eta$. At $\eta= 1$ and for 8.7 dB ($r = 1$) of squeezing the correlation reaches 0.96. At $\eta= 0.5$ it drops to 0.76, and recovering 0.96 then requires about 17 dB of squeezing. Bulk optics reaches 15 dB [95], while integrated sources reach 1 dB measured (4 dB inferred on chip) in silicon nitride [96] and 4.9 dB measured (11 dB inferred) in thin-film lithium niobate [97], which corresponds to end-to-end transmissions of roughly 0.3 to 0.7. As Gaussian error correction cannot remove this constraint [94], the model favors a shallow circuit architecture, possibly with a modified interferometric mesh.

At the same time, optical systems offer short times per sample at large width. Homodyne receivers run at gigahertz bandwidth [98], multimode homodyne detection reaches 60 modes [63], and every shot returns a full $N$-dimensional sample. For a Gaussian circuit a classical sampler competes at $O(N^{2})$ per sample, and larger width does not provide a large separation on its own. With the addition of the Kerr layer this changes. The branch count (Sec. II B) multiplies across modes, and phase-space expansions for stronger non-Gaussianity carry a cost that grows with the non-Gaussian content [54, 99, 100], whereas the physical sampling cost does not depend on $\kappa$. We regard a wide, shallow circuit with Kerr gates strong enough to make the branch sum expensive, benchmarked by time per sample at fixed correlation accuracy, as the target. The Gaussian boson sampler of Ref. [64], with more than 400 frequency modes used as a reservoir in which inter-mode correlations carried the learned information, shows that Gaussian optics at this width is available and that its correlations are useful for learning.

Conclusion

We proposed a photonic generative model based on continuous-variable encoding and a readout of position quadratures by homodyne detection. The model is built from pairwise squeezing, local phase adjustment, and interferometric mixing, where each layer controls the correlations in an interpretable way. We trained it by minimizing the distance between empirical correlation matrices and those generated by the model, reproducing correlations between multiple stochastic variables corresponding to log returns of financial assets. We also suggested an effective simulation method for a wide Kerr gate based on a coherent sum of Gaussian branches, which keeps the whole pipeline within symplectic simulation methods. Given that fast optical measurements provide real-valued samples from a joint multivariate distribution, this generative model represents a native sampler for certain distribution families at a small time-per-sample cost. As the next step, we consider models with significant non-Gaussian content, via stronger Kerr gates or measurement-induced nonlinearities, and the transfer of trained circuits onto photonic hardware with realistic loss. The single non-Gaussian element. Its generator $\hat{n}^{2}$, with $\hat{n}=\frac{1}{4}(\hat{x}^{2}+\hat{p}^{2})-\frac{1}{2}$, is quartic in the quadratures and has no symplectic representation. In general this breaks the Gaussian structure of the model, and in Sec. II B we represent it as a sum of Gaussian operations.

Acknowledgments

We thank Zoltán Kolarovszki and Zoltán Zimborás for insightful discussions on simulation backends. We are grateful to Raj B. Patel and Shang Yu for the discussions of experimental realizations. We also thank Yuan Wang for helping with the submission, and Pieter Kok for useful suggestions. O. K. acknowledges the support from the UK EPSRC award for QCi3 Hub (EP/Z53318X/1) and Open Plus Fellowship award “Quantum Generative Modelling: from Foundations to Applications”.

Contributions and AI statement. The original idea was proposed by O. K. and developed by D. C. L. as a part of the MSc project starter over a year ago, before agents became widespread. When finalizing the project, we used AI tools for writing parts of the manuscript and polishing the code. The authors take responsibility for results and presentation.

Gaussian states and optical gates

In this appendix we collect the phase-space conventions behind Eqs. (1) and (2), and summarize the action of the gates that enter the ansatz of Eq. (3). This appendix supports Sec. II A.

The state of $N$ optical modes lives in $\mathcal{H}=\bigotimes_{j=1}^{N}\mathcal{H}_{j}$. A state is Gaussian when its Wigner function is a Gaussian function of the phase-space variables, in which case the first moments $\mathbf{d}$ and the covariance matrix of Eq. (1) determine it completely [45, 46]. With the commutator $[\hat{x},\hat{p}]=2i$ (Sec. II A) the vacuum covariance matrix is $\Sigma/2=\mathbb{1}$, and a two-mode squeezed vacuum has $\Sigma_{12}^{(x)}=2\sinh2r$, the curve in Figure 2(a). Measuring a subset of modes yields outcomes distributed according to the corresponding Wigner marginal, and the position marginal of any state is the density of position-quadrature homodyne outcomes. Gaussian unitaries are generated by Hamiltonians at most quadratic in the quadratures and act as $\mathbf{d}\mapsto S\mathbf{d}$ and $\Sigma\mapsto S\Sigma S^{T}$ with $S$ symplectic.

In Table I we list the gates of Eq. (3) together with the rule that each one imposes on the measured correlations. The first four are Gaussian and account for the correlation mechanism of Sec. III A. The Kerr gate is

Operation	Symbol	Rule for the measured correlations
Two-mode squeezing	$S_{2}(r)$	Creates a positive correlation between the position quadratures of a pair; the magnitude grows monotonically with $r$ and saturates toward unity [Fig. 2(b)].
Single-mode squeezing	$S(r_{s})$	Rescales the quadratures of an unpaired mode; its correlations with other modes arise only after mixing in the mesh.
Phase rotation	$R(\phi)$	Multiplies the correlation of a squeezed pair by $\cos\phi$: $\phi=0$ keeps it, $\phi=\pi/2$ removes it, $\phi=\pi$ flips its sign (Fig. 7).
Beam splitter	$B(\vartheta)$	Mixes two modes; the elements mesh redistributes the pair correlations across all modes, with $\vartheta=\pi/4$ a balanced split.
Kerr gate	$K(\kappa)$	Intensity-dependent phase $e^{i\kappa\hat{n}^{2}}$; adds non-Gaussian shape to the output, while at small $\kappa$ the second moments are nearly unchanged.

Kerr decomposition: construction, settings, and checks

In this appendix we describe the integral representation of a single-mode Kerr gate and the finite-sum construction behind the branch expansion of Eq. (5), and we report the errors of the decomposition against exact Kerr evolution on Fock states. The scope is limited to the weak-nonlinearity regime used in Sec. III B. A systematic study of the representation is the subject of forthcoming work.

Integral representation. The Hubbard-Stratonovich identity of Eq. (4) linearizes the quadratic exponent of the Kerr gate at the cost of an auxiliary integral over $x$. Each $e^{i\kappa\hat{n}^{2}}$ is a phase rotation $R(x)$, which is the origin of the Gaussian branches in Eq. (5). Choosing sufficiently large integration limits and a quadrature rule, we convert the integral into a sum of Gaussian operations.

Finite-sum construction. The integrand oscillates rapidly, and to mitigate this we discretize the integral with Levin’s collocation method [87], which converts the oscillatory integral into a small linear system and yields the finite sum $e^{i\kappa\hat{n}^{2}}\approx\sum_{m=1}^{M}c_{m}R(x_{m})$ with complex weights $c_{m}$ and real nodes $x_{m}$. Applied to a Gaussian state this sum produces the branches $[G_{m}]=R(x_{m})|G\rangle$ of Eq. (5), and expectation values follow from Eq. (6). We set the limits to $[-100,100]$ and introduce a Gaussian window $e^{-\gamma_{\mathrm{w}}x^{2}}$, with $\gamma_{\mathrm{w}}=5\times10^{-4}$, which regularizes the integrand.

Approximation quality. In Figure 6 we compare the decomposition with exact Kerr evolution on Fock states, where $K(\kappa)|n\rangle= e^{i\kappa n^2}|n\rangle$. As a reference we use the Kerr gate simulated in the Fock backend of Piquasso, and track the error of the sum $s_{\kappa}(n) = \sum_m c_m e^{i x_m n}$ against the exact phase. The plotted scalar error is

which for $n \ge1$ lies between 1.0 and 1.5 times the absolute error $|e^{i\kappa n^2} - s_{\kappa}(n)|$. The error grows with $n$ and with $\kappa$. At $\kappa= 0.1$ it remains below $1.3 \times10^{-3}$ over the displayed photon numbers, and at $\kappa= 0.6$ it reaches about $6 \times10^{-2}$ at $n = 8$. The test establishes that the finite sum reproduces the Kerr phases accurately for weak nonlinearity and low occupation, which is the regime of the five-mode circuit. The operating value $\kappa= 10^{-5}$ lies below the scanned range and was checked directly with the same Levin weights. The unnormalized sum carries a global factor of modulus about 2.2, so $s_{\kappa}(n)$ differs from $e^{i\kappa n^2}$ by order one at this $\kappa$, but dividing by its value at $n = 0$, which is the normalization of Eq. (6), gives an error below $3 \times10^{-7}$ for $n \le8$.

Datasets, training settings, and estimators

In this appendix we list the details needed to reproduce the training experiments of Sec. III: the construction of the two datasets and their correlation targets, the optimizer settings and parameter bounds used with the loss of Eq. (7), and the conventions of the evaluation metrics defined in Sec. II C.

Datasets. The five-asset dataset contains $T = 251$ daily log returns of AAPL, MSFT, GLD, TLT, and BTC-USD. The twenty-equity dataset contains $T = 48$ daily log returns of twenty S&P 500 constituents. Prices are daily adjusted closes from Yahoo! Finance via yfinance [89]. Log returns are formed per asset, standardized to zero mean and unit variance over the window, and the target $C_{\mathrm{target}}$ is their sample Pearson correlation matrix. The five-asset window runs from 2020-03-02 to 2021-03-01, with 252 closing prices and $T = 251$ returns. BTC-USD closes are sampled on the equity trading days, so its returns span weekends and holidays. The twenty-equity window runs from 2019-11-01 to 2020-01-13, with 49 closing prices and $T = 48$ returns, and the tickers in model order are AAPL, MSFT, GOOGL, NVDA, ED, NEM, VZ, XOM, CVX, JPM, BAC, P G, KO, JNJ, UNH, NEE, DUK, BA, WMT, and TSLA.

Training. All circuits are trained with L-BFGS-B on Eq. (7) with kernel bandwidth $\gamma= 1$ under box constraints that keep the simulator states in a numerically stable range: the local phases $\phi$ are unconstrained, the two-mode and single-mode squeezing amplitudes $r$ and $r_s$ lie in $[0,1.5]$, and the beam-splitter angles $\vartheta$ of the Clements mesh lie in $[0,\pi]$. The beam-splitter internal phases are fixed at zero, and the Kerr strength is fixed at $\kappa= 10^{-5}$ in the five-mode circuit and $\kappa= 0$ in the twenty-mode circuit. Each circuit is trained from a fixed initialization and from random perturbations of it. The fixed initialization sets all local phases to 0, all beam-splitter angles to $\pi/4$, and all squeezing amplitudes to 0.35 in the five-mode circuit and to 0.70 in the twenty-mode circuit. It is the run shown in Figure 3, Figure 4, and 5. The random initializations add to the phases and the beam-splitter angles uniform offsets in $[-0.3,0.3]$ and to the squeezing amplitudes uniform offsets in $[-0.1,0.1]$, with the amplitudes and angles clipped to their bounds. The five-mode circuit is trained from ten random initializations for up to about 4600 function evaluations, and the twenty-mode circuit from five random initializations for up to about 15000 function evaluations. The convergence plots count every function evaluation. The resulting convergence traces are shown in App. § 3.8.

Estimators. $MSE_{\mathrm{full}}$ is averaged over all $N^2$ entries, including the diagonal, which contributes zero. $TVD$ is computed between histograms of target and generated samples on a common binning of 50 equal-width bins spanning the combined range of both samples, and the KS test is the two-sample test on the same samples. The generated marginals use $10^4$ homodyne samples drawn from the normal distribution with covariance $\Sigma^{(x)}$, and the target marginals use the first 70% of the standardized returns (175 days, standardized with the mean and standard deviation of this training split).

Additional results

In this appendix we collect results that support the findings in the main text: the phase control of the two-mode correlation (Sec. III A), the marginal statistics of the five-mode model (Sec. III B), the residual distribution of the twenty-mode model (Sec. III C), and the convergence of the training runs at both sizes.

Phase control. Figure 7 shows the second control discussed in Sec. III A. For a two-mode squeezed vacuum at fixed squeezing, a phase rotation applied to one mode and swept over one period modulates the measured position-quadrature correlation as a cosine of the phase. Together with Figure 2, this shows that squeezing sets the magnitude of a pairwise correlation and a local phase sets its sign.

Five-mode marginal statistics. Table II quantifies the marginal comparison of Figure 4 with the TVD and KS statistics of Sec. II C, before and after the affine map of Eq. (8). The affine map lowers the asset-averaged TVD from 0.370 to 0.245 and the averaged KS statistic from 0.190 to 0.121. The largest post-affine TVD, 0.330, belongs to BTC-USD, the most volatile asset. For GLD the raw marginal has the smallest TVD (0.165) and a KS $p$-value of 0.43, and the affine map raises its TVD to 0.225 and lowers the $p$-value to 0.14. For TLT the affine map changes neither the TVD (0.273 and 0.271) nor, appreciably, the $p$-value (0.0085 and 0.0059). The residual TVD across all assets reflects the heavy tails of the returns, which second-moment control does not reproduce.

Twenty-mode residuals. Table III gives the distribution of residual magnitudes $|\Delta C_{ij}|$ over the 190 off-diagonal pairs of the twenty-mode model of Fig. 5, averaged over five random initializations. More than half of the pairs have residuals below 0.05, 86% lie below 0.10, and a single pair on average lies between 0.20 and 0.50. The seed-to-seed spread of each fraction is about two percentage points.

Convergence. Figure 8 and Figure 9 show $\mathrm{MSE}_{\mathrm{full}}$ between $C_{\mathrm{gen}}(\theta)$ and $C_{\mathrm{target}}$ during training. Nine of the ten five-mode runs fall below $10^{-2}$ within 860 function evaluations and below $10^{-6}$ within 1960, with transient excursions of individual seeds during line searches. The tenth run remains near $1.4 \times10^{-2}$ for about 3500 evaluations before it follows the others, and falls below $10^{-2}$ after 3554 and below $10^{-6}$ after 4276 evaluations. All ten runs end between $5.0 \times10^{-10}$ and $3.0 \times10^{-8}$, and the median is $4.3 \times10^{-9}$. The twenty-mode traces decrease from initial values of 0.12 to 0.15, cross 0.05 after roughly 4000 evaluations, and approach a common plateau below 0.01 after roughly 10000 evaluations. At both sizes the runs end at comparable values, which supports the multi-seed statements in Secs. III B and III C.

Asset	Before affine map			After affine map
	TVD	KS $D$	$p$-value	TVD	KS $D$	$p$-value
AAPL	0.310	0.138	0.0026	0.202	0.090	0.112
MSFT	0.406	0.230	$1.8 \times10^{-8}$	0.195	0.102	0.051
GLD	0.165	0.066	0.431	0.225	0.087	0.138
TLT	0.273	0.125	0.0085	0.271	0.129	0.0059
BTC-USD	0.696	0.390	$5.3 \times10^{-24}$	0.330	0.195	$3.4 \times10^{-6}$
Average	0.370	0.190	—	0.245	0.121	—

Residual $|\Delta C_{ij}|$	Pairs (mean)	Fraction	Cumulative
[0, 0.05)	111	$58.5\% \pm2.1\%$	58.5%
[0.05, 0.10)	53	$27.7\% \pm1.6\%$	86.2%
[0.10, 0.20)	25	$13.1\% \pm2.1\%$	99.3%
[0.20, 0.50)	1	$0.7\% \pm0.7\%$	100%