Making Trotters Sprint: A Variational Imaginary Time Ansatz for Quantum Many-body Systems

TL;DR

This paper presents a variational imaginary time ansatz for quantum many-body ground states, demonstrating improved scaling and potential mitigation of the sign problem in certain models.

Contribution

The authors introduce a novel variational wavefunction using alternating imaginary time evolution, showing superior scaling properties and entanglement growth management.

Findings

Imaginary time scales logarithmically with system size.

Alternating evolution dynamics enable exponential entanglement growth.

Potential to mitigate the sign problem in quantum simulations.

Abstract

We introduce a variational wavefunction for many-body ground states that involves imaginary time evolution with two different Hamiltonians in an alternating fashion with variable time intervals. We successfully apply the ansatz on the one- and two-dimensional transverse-field Ising model and systematically study its scaling for the one-dimensional model at criticality. We find the total imaginary time required scales logarithmically with system size, in contrast to the linear scaling in conventional Quantum Monte Carlo. We suggest this is due to unique dynamics permitted by alternating imaginary time evolution, including the exponential growth of bipartite entanglement. For generic models, the superior scaling of our ansatz potentially mitigates the negative sign problem at the expense of having to optimize variational parameters.

Equations

$$|{\psi }_P(\mathbf{\alpha },\mathbf{\beta })⟩$$

$$H$$

$$\text{expectationvalue}\mathcal{O}{\psi }_P(\mathbf{\alpha },\mathbf{\beta })=\sum _{\{s\}}\tilde{\mathcal{O}}(s)p_{\mathbf{\alpha },\mathbf{\beta }}(s)$$

Full text

Introduction– Imaginary time plays a prominent role in multiple branches of physics, including cosmology, statistical mechanics and quantum field theory. The seemingly simple replacement of real time, $t$, with its imaginary counterpart, $\tau=-it$, leads to fundamental connections between quantum theory and statistical mechanics Wick (1954). Such connections enable the efficient simulation of many quantum systems using quantum Monte Carlo techniques Handscomb (1962); Blankenbecler and Sugar (1983); Sandvik and Kurkijärvi (1991). However, for many physically interesting models, these methods suffer from the prohibitive ‘negative sign problem’ Loh et al. (1990); Troyer and Wiese (2005), which requires an exponential amount of computational resources to obtain reasonable accuracy for quantum many-body systems. Many outstanding problems in condensed matter, such as those involving high temperature superconductors or topologically ordered phases, require an understanding of complex interacting models which are unsolved with present techniques.

One class of Monte Carlo methods that can avoid the sign problem are so-called variational Monte Carlo (VMC) methods Gutzwiller (1963); *gutzwillerCorrelationElectronsNarrow1965; Ceperley et al. (1977); Otsuka (1992); Baeriswyl (1987); Becca and Sorella (2017). In VMC, one assumes a sufficiently general trial state that depends on adjustable parameters. These parameters are then chosen to minimize the energy with respect to the given Hamiltonian. Finding an effective trial state such as Jastrow Becca and Sorella (2017), matrix product states Verstraete et al. (2004); Sandvik and Vidal (2007); Vidal (2007), or neural network states Carleo and Troyer (2017); Freitas et al. (2018); Inack et al. (2018); Pilati et al. (2019), can result in efficient simulation of interacting quantum systems. The key to the success of these techniques is a well-chosen ansatz that reflects the properties of the target phase and the existence of a viable optimization scheme Umrigar et al. (1988); Sorella (1998); *sorellaGeneralizedLanczosAlgorithm2001; *sorellaWaveFunctionOptimization2005.

The recent advent of quantum computers and simulators has motivated the development of new variational approaches Jones et al. (2019). Such variational quantum algorithms involve applying a sequence of unitary operators, parameterized by several variables onto a easy-to-prepare initial state. The variables are chosen to optimize a given cost function involving the resulting wavefunction. For example, in the quantum approximate optimization algorithm (QAOA) Farhi et al. (2014); Zhou et al. (2018); Wang et al. (2018); Hadfield et al. (2019); Verdon et al. (2019) the cost function is a classical Hamiltonian encoding a combinatorial optimization problem, and the variational wavefunction is prepared by alternating between evolving with the Hamiltonian and a transverse field. The evolution times constitute variational parameters that are optimized to minimize the Hamiltonian cost function. This variational approach has been generalized for preparing both strongly correlated and highly-entangled states on near-term quantum devices Wecker et al. (2015); Ho and Hsieh (2019); Ho et al. (2019).

Motivated by the success of such variational approaches, in this Letter we propose a variational ansatz for ground states of quantum many-body systems which involves sequentially evolving with different Hamiltonians in imaginary time. In contrast to real time evolution with local Hamiltonians, which is limited by Lieb-Robinson bounds on the growth of correlation functions, imaginary time evolution does not have this constraint and can exhibit remarkable efficiency in traversing Hilbert space.

As proof of concept, we demonstrate the efficiency of our ansatz in representing the ground state of the transverse field Ising model at criticality. Whereas standard projector methods require imaginary time scaling with system size $L$ to reach the critical ground state, we show numerically that our ansatz requires time scaling logarithmically with $L$. Furthermore, we analyze how entanglement grows after each imaginary time operation in our ansatz, and we find an exponential growth that is a unique feature of imaginary time dynamics. We conclude by demonstrating that the ansatz continues to perform well in the presence of integrability-breaking perturbations, and we mention generalizations of our ansatz to other models, including those with sign problems. We envision the main purpose of this ansatz to be an efficient trial wavefunction for quantum many-body physics on (classical) computers, however, it is possible that one can also implement such imaginary time evolution natively on a quantum computer McArdle et al. (2018); Motta et al. (2019).

Ansatz– Many Hamiltonians are naturally a linear combination of two components $H=H_{A}+gH_{B}$, where $H_{A,B}$ are individually tractable to analyze. Examples include transverse field Ising, Hubbard, and the $J_{1}$–$J_{2}$ model. Motivated by the QAOA procedure, we consider the following variational imaginary time ansatz (VITA) for the ground state of $H$:

[TABLE]

where $P$ is the number of pairs of variational parameters $\boldsymbol{\alpha}=(\alpha_{1},...,\alpha_{P}),\boldsymbol{\beta}=(\beta_{1},...\beta_{P})$, $|\psi_{0}\rangle$ is an initial state, and $\mathcal{N}$ is a normalization factor. We further define the total imaginary time $\tau=\frac{1}{2}\sum_{p=1}^{P}(\alpha_{p}+\beta_{p})$. A circuit representation is shown in Fig. 1. The bang-bang, or square-pulse, style of the ansatz is optimal for quantum control in real-time QAOA as per Pontryagin’s principle Yang et al. (2017).

While VITA is applicable to any Hamiltonian, in specific cases there is explicit physical motivation for considering such an ansatz. For example, for the fermionic Hubbard model, the $P=1$ ansatz reduces to Otsuka’s generalization of the Gutzwiller variational wavefunction Otsuka (1992); Gutzwiller (1963); *gutzwillerCorrelationElectronsNarrow1965; Baeriswyl (1987), which seeks to balance single occupancy per site with itinerancy. The $P\leq 3$ case has been considered in Ref. Yanagisawa et al. (1998); *yanagisawaCrossoverWeaklyStrongly2016; *yanagisawaAntiferromagnetismSuperconductivityPhase2019 for the two-dimensional Hubbard model, but a systematic analysis of how its performance scales with system size and $P$ was not carried out. A related variational approach for the Hubbard model has also been considered in Ref. Vaezi and Vaezi (2018).

The standard projector method for attaining the ground state of $H$ involves evaluating $e^{-\tau H}|\psi_{0}\rangle$ for $\tau\gtrsim 1/\Delta$ where $\Delta$ is the many-body spectral gap. This can be decomposed via Trotterization into a sequence of the VITA form, with parameters $\alpha_{p}=\beta_{p}=\tau/2P$ for large $P$. This guarantees that VITA can exactly represent the ground state in the $P\rightarrow\infty$ limit. However, the projector method is especially expensive for critical systems where $\Delta\sim 1/L$, and hence $\tau$ scales polynomially with $L$. One can consider Eq. (1) as a non-uniform Trotterization with large (and variable) times steps. We will show that remarkably high fidelities can be attained even with $\tau$ that is exponentially smaller compared to the aforementioned estimate from the standard projector method.

We first present some general considerations of why such an ansatz may be efficient. It is useful to first compare with the real-time analogue, which are QAOA-type circuits involving alternating real-time evolution between two Hamiltonians. The Lieb-Robinson bound dictates that real-time evolution with local Hamiltonians can generate correlations only within a light cone, and thus there are lower bounds on the time it takes to prepare highly correlated states starting from unentangled product states. For example, in one dimension, the total time to prepare the GHZ (“cat”) state $\frac{1}{\sqrt{2}}\left(|11\ldots 1\rangle+|00\ldots 0\rangle\right)$ scales at least linearly with system size $L$ Bravyi et al. (2006); Ho and Hsieh (2019). In contrast, by evolving with the GHZ parent Hamiltonian $H_{ZZ}=-\sum_{i=1}^{L}Z_{i}Z_{i+1}$ in imaginary time ($Z$ is the Pauli-$Z$ matrix), the GHZ state can be prepared with imaginary time scaling as $\log L$ 111See Supplemental Material.

Application– We test VITA on the transverse-field Ising model (TFIM)

[TABLE]

with periodic boundary conditions on a system with $N$ spins. $Z,X$ are the Pauli matrices, and $h$ is the transverse field strength. Our ansatz in this case starts from the paramagnetic ground state of $H_{X}$, $|+\rangle$ and alternates between $H_{A}=H_{ZZ}$ and $H_{B}=H_{X}$.

The Ising chain can be mapped to free fermions via the Jordan-Wigner transformation Lieb et al. (1961) which allows for the efficient evaluation of Eq. (1). We can thus optimize our ansatz for very large system sizes and several pulses. This allows us to properly characterize how efficient the VITA ansatz is without introducing sampling error. For a fixed $P$, optimization involves finding the minima of the energy cost function $E_{P}(\boldsymbol{\alpha},\boldsymbol{\beta})=\expectationvalue{H}{\psi_{P}(\boldsymbol{\alpha},\boldsymbol{\beta})}$.

We first focus on approximating the critical ground state of $H$ at $h=1$. Figure 2a shows the relative error in energy $\epsilon_{\rm rel}=|(E_{P}(\boldsymbol{\alpha},\boldsymbol{\beta})-E_{\rm exact})/E_{\rm exact}|$ where $E_{\rm exact}$ is the exact ground state energy at the critical point, for various $P$. Evidently, increasing $P$ dramatically improves the accuracy in the energy, even for large system sizes. Optimized parameters for various $P$ are provided in supplementary material.

Since the exact ground state for the TFIM is known, we also compare the fidelity, $f$, of the optimized trial state with the target state. The error in fidelity, $1-f\equiv 1-|\langle\psi_{\rm exact}|\psi_{P}(\boldsymbol{\alpha},\boldsymbol{\beta})\rangle\rangle|^{2}$, is shown in Fig. 2b for various $P,L$. The efficiency is quite remarkable; for example, for $L=64$, $P=2$ is already sufficient to approximate the critical state to within around $10^{-4}$ in relative energy and $10^{-2}$ in fidelity. Recall that the error in fidelity provides an upper-bound for the error in any observable 222See Supplemental Material.

Entanglement dynamics– While there is no Lieb-Robinson bound limiting the rate for generating long-range correlations in our ansatz, entanglement considerations provide lower bounds on the iterations $P$ required to prepare the critical state. Ignoring normalization, imaginary time evolution with a local Hamiltonian can be represented by a (non-unitary) quantum circuit; each iteration of our ansatz corresponds to three layers shown in Fig. 1. After $P$ pulses, the bipartite entanglement entropy between the left and right halves the system, hereafter abbreviated EE, can be attained by bisecting the circuit through $P$ bonds. Hence, EE after $P$ iterations is at most $P\log{D}$, where $D$ is the Schmidt rank (number of singular values) upon decomposing a single two-qubit imaginary time operator.

In order to generate the EE of the critical state, which scales as $S\propto\log{L}$, we need the number of pulses scaling at least as $P\propto\log{L}$. We observe that Fig. 2b is consistent with this scaling form. For example, for a target fidelity error of $10^{-10}$, each additional $P$ in the ansatz can represent a system approximately twice as large.

The entanglement dynamics in imaginary time evolution can be considerably different from its real time counterpart. For real time evolution, EE across a bipartition can increase only by acting with an operator supported on both sides of the partition. If the circuit in Fig. 1 were unitary, any increase in EE from one layer to the next would be bounded by a constant depending on the two-qubit unitary but not on the state being acted on (see the “small incremental entangling theorem” Mariën et al. (2014)). In contrast, even imaginary time operators acting on one side of the bipartition can generate entanglement across the cut. As a very simple example involving two spins, the action of $e^{-\beta Z_{1}}$ on $\eta_{+}|11\rangle+\eta_{-}|00\rangle$ can increase EE as long as $|\eta_{+}|>|\eta_{-}|>0$. This illustrates that the more entangled the initial state, the more imaginary time evolution can change the entanglement. The change is not simply bounded by a state-independent constant. This allows in principle an exponential growth of EE, as long as the total EE after $P$ steps lies below $P\log{D}$.

Our ansatz exhibits such dynamics. We take the $P=5$ ansatz and analyze the EE of the states at intermediate steps, $p$, of our protocol using the technique of Peschel (2003); Vidal et al. (2003); *latorreGroundStateEntanglement2004. We find that the EE increases exponentially with imaginary time (Fig. 3a). Moreover, for every intermediate state, we plot the mutual information ($S_{A}+S_{B}-S_{AB}$) between two spins $A,B$ as a function of their separation (Fig. 3b). The power law decay for every step is in stark contrast to any local real time evolution and illustrates the ability of imaginary time evolution to generate long-range correlations 333See Supplemental Material at [URL will be inserted by publisher] for considerations of a local quench in imaginary time in a conformal field theory . We find that under imaginary-time evolution, entanglement starts to grow immediately following a local quench, in contrast to real-time evolution Calabrese and Cardy (2005); *calabreseEntanglementCorrelationFunctions2007 where it takes a time proportional to $\ell$ before growing, $\ell$ being the distance between the location of the local quench and the entanglement cut.

Scaling– To compare directly with the projector method, we also investigate the total imaginary time $\tau\equiv\frac{1}{2}\sum_{p=1}^{P}(\alpha_{p}+\beta_{p})$ required to achieve a target fidelity. Motivated by the $P\propto\log{L}$ scaling for achieving a target fidelity, we propose a scaling form of $1-f=G(\tau(\log L)^{-\nu})$ for some exponent $\nu$. We perform a scaling collapse for $L\in[4,6,\ldots,1024]$, and $P\in[1,\ldots,7]$ and find the optimal exponent $\nu=2.3\pm 0.1$. This logarithmic scaling is an exponential speedup compared to the linear scaling of the projector method, $1-f=F(\tau L^{-1})$ Liu et al. (2018).

Monte Carlo approach– While the TFIM model admits a dual representation as free fermions, for a general model sampling methods are crucial for estimating the energy cost function. As a proof of concept, we also use Monte Carlo sampling for stochastically optimizing VITA.

The quantum-classical correspondence maps quantum observables to dual classical observables of a classical anisotropic Ising model on an $L\times(2P+1)$ lattice. We denote the classical spin configurations by $\{s\}$ and the spatial and imaginary time by $(i,p)$, respectively 444 This is a slight abuse of notation since $1\leq p\leq P$ in Eq. (1). However, this may be permitted since there are only $p$ unique time couplings $J_{t}(p)$ due to the lattice symmetry..

The expectation value of a quantum observable $\mathcal{O}$ is

[TABLE]

where $\tilde{\mathcal{O}}$ are dual classical observables, and $p_{\boldsymbol{\alpha},\boldsymbol{\beta}}(s)$ is the Boltzmann weight corresponding to the Ising model with couplings $J_{x}(p)=\alpha_{p},J_{\tau}(p)=\frac{1}{2}\ln\coth\beta_{p}$ between nearest neighbors in space and imaginary time, respectively. The couplings vary with imaginary time, $p$, but are uniform in space. In this way, the energy of the trial wavefunction can be sampled efficiently with Monte Carlo.

We use this scheme with $P=1,2$ to target the ground states for various values of the transverse field $h$ in both the (integrable) one-dimensional and the (non-integrable) two-dimensional TFIM. We use stochastic natural gradient descent (stochastic reconfiguration) to optimize the parameters Sorella (1998). We find rapid convergence for $P=1$, while higher-$P$ becomes more difficult, especially for with a noisy objective function.

The relative error in energy achieved is shown in Fig. 5, with the free fermion results for comparison. For $P=1$ the VMC achieves the same accuracy as the free fermion method. However, for $P=2$ away from the critical point $h=1$, the VMC performance is limited by sampling error 555 Of course the statistical error in Monte Carlo can be made arbitrarily small given long enough run-time since the sampling error goes as $\mathcal{O}(N_{\rm MC}^{-1/2})$ for $N_{\rm MC}$ Monte Carlo sweeps. This guarantees that VMC will converge to the free fermion solution in Fig. 5a if the globally optimal parameters can be found. The computational time of the MC sampling is $\mathcal{O}(\tau_{\rm corr}N_{\rm MC})$ where the autocorrelation time $\tau_{\rm corr}\sim(PN)^{\gamma}$ is a polynomial function of $PN$ for $N$ spins and $P$ pulses of VITA. .

Discussion– We have introduced a variational technique that is motivated by both projector methods and recently developed quantum algorithms. It provides substantial shortcuts to the usual Trotterization of imaginary time evolution, at the expense of making the procedure variational. Using TFIM as a first testbed, we have demonstrated that this ansatz is viable for sampling methods and highly efficient. In particular, the number of variational parameters required to represent the TFIM critical state scales as $\sim\log{L}$, in contrast to other variational methods such as density matrix renormalization group (DMRG), which in this critical case requires bond dimension scaling with $L$ and thus number of parameters scaling with $L^{2}$. One reason for this efficiency is the fact that imaginary time evolution is not subject to many bounds for real time evolution; despite being generated by local Hamiltonians, our ansatz exhibits an exponential growth of entanglement entropy and rapid generation of long-range correlations, features unique to imaginary time evolution.

Our variational approach is potentially useful in the many situations where imaginary time Trotterization involves a prohibitively large number of steps. For example, in models with a sign problem, the computational cost scales exponentially with space and imaginary time $\mathcal{O}(\tau L^{d})$. Our ansatz provides a variational shortcut that significantly reduces $\tau$ (from $\tau\sim L$ to $\tau\sim(\log L)^{2.3}$ in the critical 1d TFIM) which could enable the study of larger systems even with a sign problem. For few variational parameters ($P=1,2$), the optimization may be feasible, and we leave these investigations to future work.