Scalable quantum simulation of continuous-time generative models via tensor networks

Nathan X. Kodama; L. Andrew Wray; Sam Cochran; Chad Rigetti; Shravan Veerapaneni; Michael J. Keiser

quant-phAug 22, 2026

Scalable quantum simulation of continuous-time generative models via tensor networks

Authors: Nathan X. Kodama, L. Andrew Wray, Sam Cochran, Chad Rigetti, Shravan Veerapaneni, Michael J. Keiser

Organizations: Sygaldry Technologies, Inc., Ann Arbor, MI, USA. · University of Michigan, Ann Arbor, MI, USA.

Abstract

Continuous-time flow and diffusion models are widely used across many application domains, from large-scale deployment in computer vision and protein folding to emerging adoption for modeling language, time series, and quantum states. After training, inferring statistical properties from continuous-time models is costly. Wavefunction flows target this cost by recasting learned transport as unitary evolution, whose final Born distribution approximates the target distribution. This prepares a coherent amplitude encoding (a qsample) that can be post-processed by quantum algorithms offering a quadratic advantage over Monte Carlo sampling. We present the first numerical study of these flows, in which we represent time-dependent potentials and states as tensor networks. At spatial dimension $d=8$ , storage falls by $\sim 10^7\times$ relative to the dense grid of $N^d$ points, and evolution wall-clock time falls by $\gtrsim 10^3\times$ against a baseline extrapolated from the measured $d\le 5$ scaling. We validate our pipeline by reproducing the $O(1/\sqrt{p_{\rm rare}})$ scaling of rare-event sampling.

Figures & tables

Figure 1: Wavefunction flows and their tensor network realization. (a) A Gaussian source wavefunction

\Psi_{0}

is evolved through an intermediate state

\Psi_{t}

into a target

\Psi_{T}

via the Trotterized unitary

U(T)=\mathcal{T}\exp(-i\int_{0}^{T}\hat{\mathcal{H}}^{\mathrm{c}}_{t}\,dt)

\hat{\mathcal{H}}^{\mathrm{c}}_{t}=i[\hat{K},{V_{t}}]

(time-ordered because

{V_{t}}

, hence

\hat{\mathcal{H}}^{\mathrm{c}}_{t}

, is time-dependent), along the generation-time axis

t

. Each Trotter step alternates the kinetic step

\hat{K}

(separable, applied via one-dimensional Fourier transforms) with the position-diagonal potential step

e^{i\beta{V_{t}}}

. The final Born distribution supplies the coherent amplitude representation required by amplitude amplification algorithms when a suitable controlled preparation circuit is available. (b) A single velocity potential

V_{t}

, learned once by joint action matching, drives two realizations of the same transport. The probability flow realization integrates the ordinary differential equation (ODE)

\dot{x}_{t}=v_{t}(x_{t})

v_{t}=\nabla V_{t}

, to transport samples; the wavefunction flow realization evolves the wavefunction under

\hat{\mathcal{H}}^{\mathrm{c}}_{t}=i[\hat{K},{V_{t}}]

with tensor network compression and samples

|\Psi_{T}|^{2}

. This work develops and benchmarks a classical tensor network simulation of the wavefunction realization. (c) State encoding for the flow wavefunction

{\Psi_{t}}

: a matrix product state (MPS) with one rank-

\leq D

core per spatial axis (physical dimension

N

), so the network has one site per spatial dimension

d

and memory

\mathcal{O}(dND^{2})

in place of the

N^{d}

dense grid. (d) The position-diagonal potential step

e^{i\beta{V_{t}}}

is applied from the real generator

{V_{t}}

rather than as a generic matrix product operator (MPO). Tensor cross interpolation (TCI) interpolates

{V_{t}}

, and a time-dependent variational principle (TDVP) step evolves

{\Psi_{t}}

under it (the two-site variant 2TDVP, and the one-site variant with subspace expansion 1TDVP+SE, grow the bond adaptively; one-site 1TDVP alone is fixed-rank); a trained MPS

{V_{t}}

supplies the cores directly and removes the runtime TCI construction.

Figure 2: Dense-grid wavefunction evolution at $N{=}64$ across generation time for the two-dimensional Swiss roll target. Top row: the wavefunction $\Psi_{t}$ visualized by hue $\!=\!\arg\Psi_{t}$ and brightness $\!=\!|\Psi_{t}|/\max|\Psi_{t}|$ . Bottom row: the probability mass $|\Psi_{t}|^{2}$ . Columns show $t\in\{0,\,0.40,\,0.70,\,0.85,\,0.95,\,1.00\}$ , spaced more finely toward $t{=}1$ where the transport concentrates. The top-row hue ( $\arg\Psi_{t}$ ) stays near $0$ throughout, consistent with the trivial-phase evolution of the wavefunction flow (Supplementary Section S1.1 ). The dynamics under $\hat{\mathcal{H}}^{\mathrm{c}}_{t}=i[\hat{K},V_{t}]$ smoothly transport the Gaussian source into the prescribed target distribution. SW and MMD reported in the remaining figures of the paper are computed against samples drawn from $|\Psi(t{=}1)|^{2}$ .

Figure 3: Probability flow view of the two-dimensional orthogonal Gaussian mixture target. Source samples $x_{0}\sim\mathcal{N}(0,\sigma_{0}^{2}I)$ are integrated under $\dot{x}_{t}=v_{t}(x_{t})$ , with $v_{t}=\nabla V_{t}$ taken from the analytic potential of Supplementary Eq. ( S13 ). (Left) Sample positions at four time points $t\in\{0,0.33,0.67,1.0\}$ , colored by time and overlaid on the faded target density. (Right) Trajectories of the same samples from $t=0$ to $t=1$ as gray lines, with $t=1$ endpoints highlighted; the radial structure of the transport is set by the four orthogonal modes of the target. This is the classical limit of the wavefunction flow: $|\Psi_{T}|^{2}=p_{1}$ as an exact identity in the continuum (see Supplementary Section S1.1 ).

Figure 4: Error scaling of the dense-grid wavefunction flow matches the predicted decomposition. Two grid-free error sources, isolated on exact

d{=}2

evolutions (no tensor network truncation) and measured with a phase-aligned state metric. (a) Grid/representation error

\varepsilon_{\mathrm{grid}}=\|(1{-}P_{N})\Psi_{T}\|/\|\Psi_{T}\|

, the relative

L^{2}

Fourier-truncation error of the ideal target amplitude

\Psi_{T}\propto\sqrt{p_{T}}

, versus grid resolution

N

. For a Laplace target the amplitude

\sqrt{p}=e^{-|x|}

has a kink, so

\Psi_{T}\in H^{s}

for every

s<3/2

and the derivative discontinuity produces Fourier coefficients consistent with

\varepsilon_{\mathrm{grid}}\sim N^{-3/2}

(fit

N^{-1.47}

). (b) Trotter error

\varepsilon_{\mathrm{time}}=\|\Psi_{K}-\Psi_{\infty}\|

, the state distance between the

K

-step product-formula flow and its

K\to\infty

limit, versus Trotter steps

K

(at

N{=}16

, where

K\gg N^{2}

is reachable). The global order is

O(K^{-1})

(fit

K^{-0.94}

)—the accumulation of the local

O(\Delta t^{2})

step error over

K=T/\Delta t

steps. (c) The same

\varepsilon_{\mathrm{time}}

versus dimensionality

d

(fixed

K,N

): the realized growth (

\sim d^{1.5}

) stays within the

d^{2}

coefficient upper bound. Dense-grid propagation isolates these two terms; the tensor network ansatz adds a third, compression

\sum_{k}\eta_{k}

(Section 2.4 ).

Figure 5: Low-dimensional visualization (t-SNE) overlays at

d=3,5,7

for orthogonal Gaussian mixtures. Each panel fits a panel-local t-SNE and overlays, on a magma scale from darkest to lightest, the true target samples (darkest), an independent target draw as a finite-sample reference (labeled Exact ), and the generated TCI+1TDVP samples (lightest); the sliced Wasserstein of the plotted run is annotated. That is a per-run value and need not coincide with Table 1 , which minimizes over cells after averaging across seeds. The TCI+1TDVP samples resolve the same cluster structure as the target at every

d

; the embedding is qualitative (t-SNE does not preserve global distances or density), and the quantitative agreement is the annotated sliced Wasserstein.

Method	$d=2$	$d=3$	$d=4$	$d=5$	$d=6$	$d=7$	$d=8$
Target–target (finite-sample)	0.069	0.068	0.066	0.065	0.063	0.062	0.060
JAM	0.131	0.143	0.133	0.102	0.106	0.101	0.094
Dense	0.107	0.092	0.086	0.104	–	–	–
TCI+1TDVP	0.115	0.081	0.095	0.092	0.093	0.092	0.100
TCI+2TDVP	0.118	0.084	0.098	0.099	0.111	0.091	0.100

Table 1: Best-cell sliced Wasserstein for the dense-grid baseline and TDVP V-step variants, by spatial dimension. Values are the minimum over the hyperparameter search box (mean over

\leq 3

seeds when available); – means the method is infeasible or unsampled at that

d

. The target–target reference is the SW between two independent finite samples of the target (

n=500

samples, 128 random projections; see Section 4.15 ) — a finite-sample scale, not an irreducible floor; the method rows use the same

n=500

. Cells cluster tightly near the minimum, and the gap between the two tensor network integrators is comparable to the spread of an individual cell (median

0.005

across cells, combining seed-to-seed variation with the Monte Carlo spread of the SW estimator itself), so over the searched box TCI+1TDVP and TCI+2TDVP are not distinguishable on accuracy, consistent with Section 3 , which separates them on cost. A separate replicated

n{=}10

audit under a trained velocity potential (Supplementary Section S2.2 , JAM and TCI+1TDVP only) validates, with a second metric (MMD), that TCI+1TDVP tracks the JAM baseline. That audit includes neither Dense-grid nor TCI+2TDVP, so it does not re-rank them.

Figure 6: Tensor network costs scale polynomially in

d

at fixed accuracy. Best-cell trajectory of the TDVP V-step variants (minimum mean sliced Wasserstein over the hyperparameter search box, analytic

V_{t}

) relative to the Dense-grid baseline across

d\in\{2,\ldots,8\}

; the underlying data are the same runs that produce Table 1 . (a) Compression ratio

\mathcal{R}=N^{d}/(2N\chi+(d{-}2)N\chi^{2})

, the Dense-grid amplitude count over the MPS parameter count at the achieved bond dimension

\chi

;

\mathcal{R}

rises by roughly an order of magnitude per added

d

, reaching

{\sim}10^{7}

d{=}8

. The shaded band in all three panels marks the regime where the Dense-grid baseline is feasible and was run (

d\leq 5

); every tensor network marker is a measured run. (b) Evolution acceleration factor

\mathcal{S}=t_{\mathrm{Dense\text{-}grid}}/t_{\mathrm{MPS}}

, the Dense-grid evolution wall-clock time over the tensor network’s, so that values above unity mean the tensor network is faster. The wall-clock time is summed per Trotter step and excludes the per-step metric callback; that callback is under

3\%

of the enclosing loop at

d\geq 4

but over

90\%

d{=}2

, enough to invert the apparent ordering of

d{=}2

and

d{=}3

if included. TCI+1TDVP crosses above unity near

d{=}4

and TCI+2TDVP near

d{=}5

, after which the tensor network methods are cheaper than Dense-grid at matching accuracy, the one exception being TCI+2TDVP at

d{=}6

(

0.9\times

). The Dense-grid numerator is measured at

d\in\{3,4,5\}

and is a log-linear fit through those three points elsewhere, the

N^{d}

grid being infeasible above

d{=}5

; the

d{=}6

exception falls in the projected regime. The tensor network denominators are measured from runs at each shown

d

. (c) Best-cell accuracy

\mathrm{SW}(\Psi_{T})

(sliced Wasserstein of the endpoint state) versus

d

: TCI+1TDVP and TCI+2TDVP stay within a factor of

1.7

of the target–target reference across the full range, without upward drift in

d

, and track the Dense-grid baseline wherever it is feasible. The cost–accuracy Pareto frontiers are given in the supplement (Supplementary Section S2 , Figure S1 ); a replicated

n{=}10

audit with a second metric appears in Supplementary Section S2.2 .

Figure 7: Deviation from the probability flow grows as the flow reaches the target and the V-step integrator sets the cost–accuracy balance. Gaussian mixture at

d{=}8

\sigma{=}0.5

; the trained MPS potential cores are fed directly to the V-step, with no runtime TCI. Two encodings run through panels (a–c): color denotes the Trotter count

K\in\{40,80,120,160\}

(dark to light), and line style denotes the integrator—pure 2TDVP (two-site update every Trotter step, thick lines) versus 1TDVP with demand-driven subspace expansion (1TDVP+SE, thin lines); panel (d) instead marks the integrator by symbol (circles: 1TDVP+SE; squares: 2TDVP). (a) Instantaneous deviation from the probability flow,

\delta\mathrm{SW}(t)=\mathrm{SW}^{\mathrm{wf}}({\Psi_{t}})-\mathrm{SW}^{\mathrm{pf}}(p_{t})

, where

\mathrm{SW}(q)

is the distance from samples of the source

q

to samples of the target—

{\Psi_{t}}

sampled by Born measurement, the marginal density

p_{t}

by transport—the two track closely until late in the flow; (b) its cumulative integral

\int_{0}^{t}|\delta\mathrm{SW}|\,dt^{\prime}

. (c) Bond dimension

\chi_{\max}(\Psi_{t})

—pure 2TDVP self-limits

\chi

by truncation (

\chi^{\star}\!\approx\!16

) while 1TDVP+SE grows

\chi

at the larger

K

needed to suppress Trotter error. (d) Total deviation

\Delta\mathrm{SW}=\int_{0}^{1}|\delta\mathrm{SW}|\,dt

versus measured wall-clock time: a cost–accuracy frontier—pure 2TDVP attains lower deviation at a modest (

\sim\!2

–

7\times

) wall-clock cost over 1TDVP+SE (up to

\sim\!10

–

40\times

at small

K

$d$	params	SW ${}^{\mathrm{wf}}_{T}$ (mean $\pm$ 95% CI)	SW ${}^{\mathrm{pf}}_{T}$	$\chi^{\star}$
8	15k	0.0363 $\pm$ 0.0027	0.0299	5–20
12	23k	0.0341 $\pm$ 0.0026	0.0293	6–31
16	31k	0.0380 $\pm$ 0.0020	0.0291	6–19
32	64k	0.0306 $\pm$ 0.0009	0.0250	7–11

Table 2: Trained MPS potential tensor network scaling on the orthogonal Gaussian mixture. The learned real-valued MPS velocity potential (bond

D{=}8

) fed directly to the 2TDVP V-step, eliminating the runtime TCI construction; fixed

K{=}160

Trotter steps across

d

. Each row averages

n{=}10

independent initializations of the trained MPS potential. Columns: spatial dimension

d

; oracle parameter count; endpoint sliced Wasserstein SW

{}^{\mathrm{wf}}_{T}

(mean

\pm

95% CI); the probability flow action matching floor SW

{}^{\mathrm{pf}}_{T}

(

\dot{x}=\nabla V

); and the achieved bond

\chi^{\star}

Figure 8: Rare-event sampling cost and yield on the

d{=}8

orthogonal Gaussian mixture (

\sigma{=}0.5

>\!4\sigma

a\approx 4.24\%

). Both pipelines are driven by the same learned velocity potential—the flow ODE pipeline integrates

\dot{x}=\nabla V

of the oracle whose cores the wavefunction flow consumes—so state preparation error is common to both pipelines. The comparison isolates the transport and the sampling method. Both are classical simulations. Cost is quoted in state preparations, i.e. complete transports. (a) Sampling cost. State preparations per accepted rare sample against the rare-event probability

{p_{\mathrm{rare}}}

, each pipeline plotted at its own measured tail mass. Markers are measurements and lines are the cost models. Prefactors are fitted with the exponents held at their structural values: rejection over all points, amplification on the tail (

{p_{\mathrm{rare}}}\leq 0.2

, i.e.

k=3.5

–

6.5

) where its

1/\sqrt{{p_{\mathrm{rare}}}}

law is asymptotic; the fitting protocol and the round-count rule are given in Methods. (b) Rejection sampling. Drawn from the flow ODE state: at a budget of

Q{=}500

state preparations and a

4\sigma

threshold, it yields

Qp=42

rare samples. (c) Amplitude amplification. Drawn from the wavefunction flow state:

Q/\mathrm{cost}=104

rare samples at the same budget, a factor

2.48

more, which reproduces the lift of panel (a) up to the rounding of integer counts. Each method is thus shown on the state its own pipeline prepares. The two panels come from a single t-SNE fit applied to both clouds at once, so they share a frame, and the ring densities compare directly; color encodes the nearest mode. Untransported residual mass—closer to the origin than to any mode—is excluded throughout, matching the oracle’s indicator; it is

3.0\%

of the wavefunction flow’s mass, against

0.2\%

for a finite sample of the target itself. Tail fractions are measured on

4\times 10^{4}

Born samples per pipeline; the Grover outcome statistics are analytic (Methods), with no compilation, controlled state preparation, or tensor network contraction costs included.

$d$	$N_{\max}$	$K_{\max}$	$D_{\max}$
2	256	64	128
3	64	64	64
4	32	64	64
5	32	64	32
6	32	64	32
7	32	64	32

Table 3: Hyperparameter search-grid maxima per

d

. The

N

cap is set by container memory (dense-grid path); the

D

cap is the largest bond dimension we Sobol-sampled.

Figure S1: Cost–accuracy Pareto frontiers across spatial dimensions. Top row (panels a–c): SW vs MPS/Dense-grid memory ratio at

d=3

4

5

. Bottom row (panels d–f): SW vs evolution wall-clock time (s) at

d=3

4

5

— summed per Trotter step, the same quantity as Figure 6 (b), so the two are directly comparable. Each point is a hyperparameter cell; staircases are per-method Pareto fronts; dashed lines mark the Dense-grid reference SW and dotted lines the JAM reference. JAM and the TDVP variants occupy the lower-left frontier in every panel. The MPS-to-Dense-grid advantage on both axes widens as

d

increases.

method	$d{=}2$	$d{=}3$	$d{=}4$	$d{=}5$
Target–target (finite-sample)	$0.023{\pm}0.001$	$0.025{\pm}0.002$	$0.023{\pm}0.001$	$0.024{\pm}0.001$
JAM	$0.111{\pm}0.006$	$0.084{\pm}0.005$	$0.073{\pm}0.004$	$0.077{\pm}0.004$
TCI+1TDVP	$0.206{\pm}0.007$ (0.051)	$0.142{\pm}0.006$ (0.066)	$0.146{\pm}0.007$ (0.095)	$0.156{\pm}0.005$ (0.065)

Table S1: Replicated

n{=}10

best-cell audit (trained velocity potential), sliced-Wasserstein as mean

{\pm}

std over the

10

seeds, with the seed-mean MMD (a point value, no spread) in parentheses; lower is better. The best cell

(N,K^{*},D^{*})

per entry is located by coordinate descent; the located corners are interior (

K^{*}{\sim}40

D^{*}{\sim}8

), not at the box maxima.

Figure S2: Sharpness route: instantaneous deviation from the probability flow versus mode width. One panel per mode width

\sigma\in\{0.35,0.5,0.75,1.0\}

at fixed

N{=}32

(orthogonal Gaussian mixture, analytic

V_{t}

d{=}8

); color encodes the Trotter count

K

, with 1TDVP+SE drawn light and 2TDVP dark. The deviation

\mathrm{SW}^{\mathrm{wf}}(\Psi_{t})-\mathrm{SW}^{\mathrm{pf}}(t)

grows as the modes sharpen (smaller

\sigma

, smaller

\sigma/\mathrm{d}x

Figure S3: Resolution route: instantaneous deviation versus grid resolution. As in Figure S2 but at fixed

\sigma{=}0.5

with the grid

N\in\{22,32,48,64\}

varied, reaching the same

\sigma/\mathrm{d}x

ratios by refining

\mathrm{d}x

. The deviation grows with

N

, opposite to the sharpness route at matched

\sigma/\mathrm{d}x

: the two routes do not collapse because a finer grid also raises the Trotter constant.

Figure S4: Cost–accuracy across Trotter count and mode sharpness. Total deviation versus measured wall-clock time (orthogonal Gaussian mixture, analytic

V_{t}

d{=}8

N{=}32

). Color encodes the Trotter count

K

, marker shape the mode width

\sigma

, and fill the

V

-step (open: 2TDVP; filled: 1TDVP+SE). Lower-left is better; sharper modes and larger

K

both move points up and to the right.

Figure S5: Analytic versus trained

V_{t}

\sigma{=}0.5

. Orthogonal Gaussian mixture,

d{=}8

N{=}32

. Columns: analytic Gaussian mixture velocity (left) and pre-trained MPS potential (right). Top row: instantaneous deviation from the probability flow,

\delta\mathrm{SW}(t)=\mathrm{SW}^{\mathrm{wf}}({\Psi_{t}})-\mathrm{SW}^{\mathrm{pf}}(p_{t})

, where

\mathrm{SW}(q)

is the distance from samples of the source

q

to samples of the target (

{\Psi_{t}}

by Born measurement, the marginal density

p_{t}

by transport), with

K

by color. Bottom row: the deviation accumulated over the flow,

\Delta\mathrm{SW}=\int_{0}^{1}|\delta\mathrm{SW}|\,dt

, versus wall-clock time (1TDVP+SE and 2TDVP). The trained oracle reproduces the analytic-

V

flow, isolating velocity learning error as small.

Quantum Generative Modeling

Quantum generative modeling uses quantum circuits and quantum-inspired methods to learn and sample data distributions and quantum states.

See all 54 papers

An Efficient Quantum Circuit for Flow Model Execution Using Quantum Neural Networks
Oct 8, 2026Rui Che, Ludvig af KlintebergQuantum Generative Modeling Quantum Circuit Simulation
Department of Business and Mathematics, Mälardalen University, Universitetsplan 1, Väster˚as, 721 23, Sweden.
Quantum Flow Matching
Aug 17, 2025Zidong Cui, Pan Zhang, Ying TangFlow Matching Quantum Generative Modeling
Institute of Fundamental and Frontier Sciences, University of Electronic Science and Technology of China, Chengdu 611731, China · Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China · School of Fundamental Physics and Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China +3
Diffusion-warm sampling of the XY model enables fast thermalization at scale
Jun 29, 2026Sehmimul Hoque, Roger Melko, Pooya RonaghMarkov Chain Monte Carlo Diffusion Sampling
Institute for Quantum Computing, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada · Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada · Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada +1

Oct 8, 2026J/K move · Enter open · S save

An Efficient Quantum Circuit for Flow Model Execution Using Quantum Neural Networks

Rui Che, Ludvig af Klinteberg

Department of Business and Mathematics, Mälardalen University, Universitetsplan 1, Väster˚as, 721 23, Sweden.

Quantum Generative Modeling Quantum Circuit Simulation Flow Matching

Open paper AskPDF arXiv

Flow models generate trajectories from an initial distribution to a target distribution by solving an ordinary differential equation defined by a velocity field. Flow matching learns this velocity field by modeling the transport dynamics between the two distributions. Wavefunction flow establishes a formal connection between flow models and quantum dynamics by introducing a continuity Hamiltonian, which drives the Schrödinger evolution of quantum states. In this paper, we investigate accurate and efficient quantum simulation of the wavefunction flow, thereby realizing the efficient implementation of flow models on quantum computers. We first leverage a quantum read-only memory (QROM)-based phase kickback framework for the wavefunction flow simulation, generating probability densities that closely match those produced by the corresponding conventional flow model. To address the high circuit-resource cost, we further incorporate a trained quantum neural network (QNN) into the phase kickback framework, replacing QROM for data encoding. Numerical experiments demonstrate that our proposed method implements flow models on quantum computers more efficiently, since it maintains the accuracy of wavefunction flow simulation compared with the QROM-based framework, and significantly reduces the circuit resources.

Abstract

Quantum Generative Modeling

Explore similar work

An Efficient Quantum Circuit for Flow Model Execution Using Quantum Neural Networks

Quantum Flow Matching

Diffusion-warm sampling of the XY model enables fast thermalization at scale