Embedding classical dynamics in a quantum computer

We develop a framework for simulating measure-preserving, ergodic dynamical systems on a quantum computer. Our approach provides an operator-theoretic representation of classical dynamics by combining ergodic theory with quantum information science. The resulting quantum embedding of classical dynamics (QECD) enables efficient simulation of spaces of classical observables with exponentially large dimension using a quadratic number of quantum gates. The QECD framework is based on a quantum feature map that we introduce for representing classical states by density operators on a reproducing kernel Hilbert space, $\mathcal{H}$. Furthermore, an embedding of classical observables into self-adjoint operators on $\mathcal{H}$ is established, such that quantum mechanical expectation values are consistent with pointwise function evaluation. In this scheme, quantum states and observables evolve unitarily under the lifted action of Koopman evolution operators of the classical system. Moreover, by virtue of the reproducing property of $\mathcal{H}$, the quantum system is pointwise-consistent with the underlying classical dynamics. To achieve a quantum computational advantage, we project the state of the quantum system onto a finite-rank density operator on a $2^n$-dimensional tensor product Hilbert space associated with $n$ qubits. By employing discrete Fourier-Walsh transforms of spectral functions, the evolution operator of the finite-dimensional quantum system is factorized into tensor product form, enabling implementation through an $n$-channel quantum circuit of size $O\left(n\right)$ and no interchannel communication. Furthermore, the circuit features a state preparation stage, also of size $O\left(n\right)$, and a quantum Fourier transform stage of size $O\left(n^2\right)$, which makes predictions of observables possible by measurement in the standard computational basis. We prove theoretical convergence results for these predictions in the large-qubit limit, $n\to \infty$. In light of these properties, QECD provides a consistent simulator of the evolution of classical observables, realized through projective quantum measurement, which is able to simulate spaces of classical observables of dimension $2^n$ using circuits of size $O\left(n^2\right)$. We demonstrate the consistency of the scheme in prototypical dynamical systems involving periodic and quasiperiodic oscillators on tori. These examples include simulated quantum circuit experiments in Qiskit Aer, as well as actual experiments on the IBM Quantum System One.

Figure 1

Schematic of the relationship between the flow map ${\mathrm{\Phi }}^t$ that advances the (nonlinear) dynamical system in a state space $X$ in time and the linear Koopman operator $U^t$ that advances observables $f$ on $X$ in an infinite-dimensional Hilbert space.

Figure 2

Schematic representation of the QECD framework applied to states and observables of a classical dynamical system in five successive levels, leading to an $n$-qubit quantum computational system. These are the classical, classical statistical, quantum mechanical, matrix mechanical, and quantum computational levels. The horizontal arrows from top to bottom in the left- and right-hand columns represent the time evolution maps of states and observables, respectively. These are the flow map ${\mathrm{\Phi }}^t$ on the classical state space $X$, the transfer operator ${\mathrm{\Phi }}_{*}^t$ on the space of probability measures $\mathcal{P}\left(X\right)$, and the Koopman operator $U^t$ on the algebra of classical observables $\mathfrak{A}\subseteq C\left(X\right)$. They are followed by the unitary evolution map ${\mathrm{\Psi }}^t$ and the Heisenberg operator ${\mathcal{U}}^t$ on the space of density operators $Q\left(\mathcal{H}\right)$ and bounded linear operators $B\left(\mathcal{H}\right)$, respectively, on the reproducing kernel Hilbert space $\mathcal{H}$. The maps at the matrix mechanical level, ${\mathrm{\Psi }}_n^t$ and ${\mathcal{U}}_n^t$, are finite-rank projections of ${\mathrm{\Psi }}^t$ and ${\mathcal{U}}^t$, respectively, acting on operators on $2^n$-dimensional subspaces ${\mathcal{H}}_n$ of $\mathcal{H}$. The corresponding maps ${\widehat{\mathrm{\Psi }}}_n^t$ and ${\widehat{\mathcal{U}}}_n^t$, respectively, at the quantum computational level act on operators on the $2^n$-dimensional tensor product Hilbert space ${\mathbb{B}}_n$, which forms the basis of an $n$-qubit quantum computer. The vertical arrows correspond to maps that translate states (left-hand column) and observables (right-hand column) to the next representation level. Under the combined action of these maps, a classical state $x\in X$ is mapped to an $n$-qubit density matrix ${\widehat{\rho }}_{x,n}\in Q\left({\mathbb{B}}_n\right)$, and a classical observable $f\in \mathfrak{A}$ is mapped to a self-adjoint operator ${\widehat{S}}_n\in B\left({\mathbb{B}}_n\right)$. A loop of arrows represents a commutative diagram.

Figure 3

Ergodic (a) and nonergodic (b) linear flows on the two-dimensional torus ${\mathbb{T}}^2$. In (a) the ratio of the frequency parameters ${\alpha }_1/{\alpha }_2$ is irrational, and the trajectory starts to fill the torus surface. In (b) the ratio of the frequencies is rational, and the trajectory is closed. The corresponding frequency parameters ${\alpha }_1$ and ${\alpha }_2$ are given to the right of each figure.

Figure 4

Structure constant matrices $c_{jl}$ for reproducing kernel Hilbert algebras on (a) the circle with $d=1$ and (b) the 2-torus with $d=2$. In both cases we use the parameter values $p=1/4$ and $\tau =1/4$ as given in (12). In (a), we consider indices in the range $-2^{n-1}\le j,l\le 2^{n-1}$ with $n=3$. In (b), the multi-indices $j=(j_1,j_2)$ and $l=(l_1,l_2)$ satisfy $-2^{n/2-1}\le j_i,l_i\le 2^{n/2-1}$ with $n=8$. In both (a) and (b), we map $j$ and $l$ into standard matrix indices $1,2,...,{(2^{n/d}+1)}^d$ [which results in ${(2^4+1)}^2=289$ for (b)] by lexicographical ordering. The matrix in (b) is thus equal to the Kronecker product of the matrix in (a) with itself.

Figure 5

Matrix elements ${\left(A_{{\psi }_1}\right)}_{ij}$ (a) and ${\left(S_{{\psi }_1}\right)}_{ij}$ (b) of the multiplication operator $A_{{\psi }_1}$ [Eq. (26)] and the self-adjoint operator $S_{{\psi }_1}$ [Eq. (27)] representing the basis function ${\psi }_1$ for dimension $d=1$. As in Fig. 4, we consider the reproducing kernel Hilbert algebra $\mathfrak{A}$ on the circle from with $p=1/4$ and $\tau =1/4$, and map the indices $i$ and $j$ into standard matrix indices $1,2,...,2^n$ with $n=3$. The matrix in (a) has nonzero elements only in the first lower diagonal, $i-j=1$. The matrix in (b) is a symmetric bidiagonal matrix with elements in the first upper and lower diagonals, $i-j=\pm 1$.

Figure 6

Quantum circuit implementation of the three- and seven-qubit approximation of a circle rotation with frequency ${\alpha }_1=2\pi$ in the ideal Qiskit Aer environment. (a) Circuit diagram with $n=3$ qubits, comprising (from left to right) of state vector load, Koopman evolution over time $t$ using ${\mathsf{R}}_z$ gates, quantum Fourier transform (QFT), and measurement. (b) Empirical distribution of an ensemble of $K={10}^6$ projective measurements (shots) of the projection-valued measure (PVM) associated with the computational basis vectors $|\mathit{b}\rangle \equiv |b\rangle$ for $n=3$ and $t=0.94$. (c) Temporal evolution of the empirical probability distributions for $n=3$ and 7. (d) Reconstruction of the classical observable $f^{\left(t\right)}\left(x\right)=sin(x\left(t\right))=sin({\theta }^1+{\alpha }_{1}t)$ from the ensemble means, ${\widehat{f}}_n^{\left(t\right)}\left(x\right)$. The analytical result $f^{\left(t\right)}\left(x\right)$ is plotted as a solid line. In panels (b)–(d), the initial condition is $x={\theta }^1=2.5$ and the reproducing kernel Hilbert algebra (RKHA) parameters are $p=\tau =1/4$. Measurements are performed at a fixed time step $\mathrm{\Delta }t=0.02$. In panels (b) and (c), the computational basis vectors $|b\rangle$ are indexed by an integer $b$ in the range $0,\cdots ,2^{n-1}$.

Figure 7

Eigenvalues $s_j$ (left-hand column) and representative eigenfunctions $u_j\text{:}{S}^1\to \mathbb{R}$ (right-hand column) of the self-adjoint operators $S_n=T_{n}f$ representing the classical observable $f\left(x\right)=sinx$ on the circle for the qubit numbers $n=3$ (top row) and 7 (bottom row). The RKHA parameters are $p=\tau =1/4$ as in Fig. 6. The index $j$ runs from 1 to $2^n$. Notice that as $n$ increases the spectra of $S_n$ provide an increasingly dense sampling of the range of values of $f$ (i.e., the interval $[-1,1]$), and the eigenfunctions $u_j\left(x\right)$ become increasingly localized around values of $x$ for which $f\left(x\right)\approx s_j$.

Figure 8

As in Fig. 6, but for an eight-qubit approximation of a quasiperiodic rotation on the 2-torus with frequency parameters ${\alpha }_1=3\sqrt{2}\pi$ and ${\alpha }_2=2\pi$. (a) Quantum circuit for the quasiperiodic system, composed as two parallel copies of the circuit in Fig. 6 for the one-dimensional case, with 4 qubits allocated to each dimension of the 2-torus. An empirical probability distribution obtained from $K={10}^6$ shots is shown to the right of the circuit diagram, where the integers $b=0,...,2^8-1=255$ index the computational basis vectors $|b\rangle$ of the 256-dimensional Hilbert space ${\mathbb{B}}_n$ with $n=8$. The RKHA parameters are again $p=\tau =1/4$. (b) Snapshots of the probability distribution at three representative evolution times, combined in a single surface plot. The horizontal axes labeled $|k\rangle$ and $|l\rangle$ correspond to the basis vector indices for each of the four-qubit spaces associated with each torus dimension through the factorization ${\mathbb{B}}_8={\mathbb{B}}_4\otimes {\mathbb{B}}_4$. Note that the indices $k$ and $l$ range from 0 to $2^4-1=15$. (c, d) Evolution of the marginal distributions obtained by measurement of the PVMs of each of the two four-qubit spaces, i.e., one of the two torus dimensions only. The initial condition is $x=({\theta }^1,{\theta }^2)=(1.0,2.5)$, and measurements are performed at a fixed time step $\mathrm{\Delta }t=0.02$. The slopes of the probability contours in panels (c) and (d) are proportional to the frequency parameters ${\alpha }_2$ and ${\alpha }_1$, respectively. Notice that the slopes in panel (c) are shallower than those in panel (d) since ${\alpha }_2<{\alpha }_1$, and are equal to the corresponding slopes in Fig. 6 since ${\alpha }_2$ is equal to the frequency parameter of the one-dimensional example. (e) Reconstruction of the classical observable $f^{\left(t\right)}\left(x\right)=f^{\left(t\right)}(x_1,x_2)=cos({\theta }^2+{\alpha }_{2}t)sin({\theta }^1+{\alpha }_{1}t)$ from the ensemble means ${\widehat{f}}_n^{\left(t\right)}\left(x\right)$ output from the quantum computer. The true evolution $f^{\left(t\right)}\left(x\right)$ is plotted as a solid line.

Figure 9

Comparison of three-qubit approximations of a circle rotation with frequency ${\alpha }_1=2\pi$ from simulated circuit experiments in the ideal Qiskit Aer environment (ideal) and actual quantum computing experiments on the IBM Quantum System One (ibmq). (a) Temporal evolution of the empirical probability similar to Fig. 6 in the ideal Qiskit circuit simulation, using amplitude encoding with QuantumCircuit.initialize for the state preparation and the RKHA parameters $p=\tau =1/4$. (b) Temporal evolution of the empirical probability distributions for $n=3$ on the quantum computer starting with a uniform superposition state $|\mathrm{\Omega }\rangle$ at $t=0$. (c) Reconstruction of the classical observable $f^{\left(t\right)}\left(x\right)=sin(x\left(t\right))=sin({\theta }^1+{\alpha }_{1}t)$ from the ensemble means, ${\widehat{f}}_n^{\left(t\right)}\left(x\right)$. The analytical result $f^{\left(t\right)}\left(x\right)$ is plotted as a solid line. In all panels, the initial condition is $x={\theta }^1=0$. Measurements are performed at a fixed time step $\mathrm{\Delta }t=0.02$. The number of shots is $K=2^{18}=\text{262}\text{144}$ in both cases.