Variational imaginary-time evolution

What it fills

• Estimate a Hamiltonian's ground-state energy

  Given a Hamiltonian you can query and some way of preparing trial states, return an estimate of its lowest eigenvalue to a stated additive error. The state achieving that energy may or may not come back with the number; almost every application wants the number.

When it applies

Yuan et al. review the variational principles this rests on and name them — the Rayleigh-Ritz method for static problems, and "the Dirac and Frenkel variational principle, the McLachlan's variational principle, and the time-dependent variational principle" for dynamics — then "discuss the connections of the three". Which principle is used is a real choice with real consequences, not a formality, and the paper is a theory of the family rather than a benchmark of one instance.

Requires

These do not move the route along. The method needs each of them alongside its own work, and the cost of getting them is part of what the method costs.

• Estimate an observable

  Given the ability to prepare $|ψ⟩$ and a description of an observable $O$, return a classical scalar within $\varepsilon$ of $⟨O⟩$ at confidence $1−δ$. The state is never returned; only the number is.

  What this route reads out is the coefficients of an $N\times N$ linear system rather than the energy alone. Imaginary-time McLachlan gives $\sum_j A^R_{i,j}\dot{\theta}_j=-C^R_i$ (Eq. 42), with $A_{i,j}=\frac{\partial\langle\phi|}{\partial\theta_i}\frac{\partial|\phi\rangle}{\partial\theta_j}$ and $C_i=\frac{\partial\langle\phi|}{\partial\theta_i}H|\phi\rangle$ (Eq. 10); the energy $E_\tau$ carried by $C'_i$ drops out, because normalisation forces $\Re\left(\frac{\partial\langle\phi|}{\partial\theta_i}E_\tau|\phi\rangle\right)=0$ and hence $C^R_i={C'}^R_i$. Expanding $H=\sum_j h_j\sigma_j$ and the gate derivatives makes each entry a sum of terms $a\Re\left(e^{i\theta}\langle\bar{0}|U|\bar{0}\rangle\right)$ with $U=R^\dagger_{k,i}R_{q,j}$ or $R^\dagger_{k,i}\sigma_jR$ (Eqs. 21-22); Fig. 1 evaluates one term, not one entry — an ancilla in $(|0\rangle+e^{i\theta}|1\rangle)/\sqrt{2}$, where $a$ and this $\theta$ are set by the coefficients and are not ansatz parameters, two controlled unitaries, a $|\pm\rangle$ measurement, and $2P_+-1$. approximation: each term is a sampled estimate; Yuan et al. say only that the precision of the estimated derivatives depends on the size of the system or the number of measurements and can affect the accuracy of the simulation, deferring the scaling to Li and Benjamin The energy itself is this same slot again, $E(\vec{\theta})=\sum_j h_j\langle\phi(\vec{\theta})|\sigma_j|\phi(\vec{\theta})\rangle$, measured term by term.

  approximation

Example

given  H = sum_j h_j sigma_j, assumed to be "a linear combination of tensor products
       of local operators sigma_j with coefficients hj"                     (Sec. 2.1)
           # the paper does NOT narrow sigma_j to Pauli strings and does not state
           #   h_j real here; "Pauli" appears only in the single-qubit numerical
           #   illustration of Sec. 3
       ansatz  |phi(theta)> = R_N(theta_N) ... R_1(theta_1) |0bar>, theta real, N params
           # "The ansatz is automatically normalised": <phi|phi> = 1 by construction
       each gate differentiable into unitaries:
           dR_k/dtheta_k = sum_i f_k,i R_k sigma_k,i                         (Eq. 18)
           # the general condition is only that the sigma_k,i be unitary and the
           #   f_k,i complex. Hermitian sigma is the paper's worked EXAMPLE
           #   (R_k = e^{-i theta_k sigma/2}, sigma Hermitian), not the requirement
           # real theta is not a restriction: the paper states real parameters
           #   suffice for quantum simulation, w.l.o.g. because the trial state can
           #   be built from gates e^{i theta_k sigma} with sigma Hermitian (Sec. 2.1)

requires  McLachlan's variational principle specifically -- for real parameters it is
    the only one of the three that is stated to guarantee a real theta_dot.
    Dirac-Frenkel leaves theta_dot possibly complex, since A and C' are complex
    (Eq. 39, Sec. 5.1.1); the time-dependent variational principle gives
    i sum_j A^I_{i,j} theta_dot_j = C^R_i, whose "solution is imaginary and hence
    not physical"                                                (Eq. 44, Sec. 5.1.3)

# --- tangent-vector circuits, defined once from the ansatz ------------------
for each gate R_k, for each term sigma_k,i in its derivative (Eq. 18):
    R_k,i = R_N R_{N-1} ... R_{k+1} R_k sigma_k,i ... R_2 R_1                (Eq. 20)
        # sigma_k,i sits between R_{k-1} and R_k, i.e. in circuit order it is
        #   applied immediately BEFORE R_k -- which is where Fig. 1 draws it
        # for most one- and two-qubit gates R_k the sum over i has only one or two
        #   terms, and sigma_k,i is itself a one- or two-qubit gate       (Sec. 2.3)
    d|phi>/dtheta_k = sum_i f_k,i R_k,i |0bar>                               (Eq. 19)

theta <- theta_0                                            # initial parameters
loop over imaginary-time steps, until E(theta) settles:  # tau = 0, dtau, 2dtau, ...
    # --- the entries wanted, written out term by term ------------------------
    A_k,q = sum_{i,j} conj(f_k,i) f_q,j <0bar| R^dagger_k,i R_q,j |0bar>     (Eq. 21)
    C_k   = sum_{i,j} conj(f_k,i) h_j <0bar| R^dagger_k,i sigma_j R |0bar>   (Eq. 22)
        # k, q index PARAMETERS; i, j index the derivative terms inside gate k or q.
        #   Eq. 42 below writes the same objects with parameter indices i, j

    for each single term of Eq. 21 or Eq. 22 -- its real and imaginary parts are
        each of the form  a * Re( e^{i*angle} <0bar| U |0bar> ), with
        U = R^dagger_k,i R_q,j  or  R^dagger_k,i sigma_j R, and with amplitude a
        and phase angle "determined by the coefficients", not ansatz parameters:
                                                                         (Sec. 2.3)
        # one Hadamard-test circuit per TERM, not per matrix entry
        prepare ancilla in (|0> + e^{i*angle}|1>)/sqrt(2); register in |Psi_0>
        run R_1, R_2, ..., R_N on the register UNCONTROLLED, and insert exactly two
            ancilla-CONTROLLED gates: U_k = one sigma_k,i, placed just before R_k,
            and U_q = one sigma_q,j (or sigma_j, by taking q = N+1 so sigma_j sits
            to the left of R_N), placed just before R_q, assuming k < q     (Fig. 1)
            # only these two local insertions are controlled; the ansatz gates run
            #   plainly. "Usually, Uk and Uq are unitary operators on only one or
            #   two register qubits rather than the entire register" -- controlling
            #   a whole R_k,i circuit is NOT what Fig. 1 does
            # the circuit is N register gates, two X gates on the ancilla, and two
            #   controlled unitaries; register gates after the second controlled
            #   unitary can be omitted                              (Fig. 1 caption)
        measure the ancilla in the |+>, |-> basis; let P_+ = prob(ancilla = |+>)
        value = 2*P_+ - 1 = Re( e^{i*angle} <0bar| U |0bar> )               (Sec. 2.3)
        # Ref. [58] is cited for reducing this to a direct measurement on the
        #   register with no ancilla; that variant is not given here

    A^R_k,q = Re(A_k,q),   C^R_k = Re(C_k)                          (Eq. 10, 21, 22)
        # this C is the plain C of Eq. 10, carrying no -E_tau term. Normalisation
        #   forces Re( d<phi|/dtheta_i * E_tau * |phi> ) = 0, so C^R equals the
        #   real part of the imaginary-time C' of Eq. 40: E_tau is never measured
        #   or subtracted here                                          (Sec. 5.1.2)

    # --- solve the projected imaginary-time equation of motion ---------------
    solve  sum_q A^R_k,q theta_dot_q = -C^R_k   for real theta_dot           (Eq. 42)
        # McLachlan's principle, delta || (d/dtau + H - E_tau)|phi(theta(tau))> || = 0
        #   (Eq. 41), applied to the Wick-rotated flow d|psi>/dtau =
        #   -(H - E_tau)|psi> (Eq. 36) projected onto the ansatz tangent space (Eq. 37)
        # the paper gives no regularisation scheme for this system when A^R is
        #   singular or near-singular; its only singularity discussion concerns
        #   A^I in the real-time TDVP branch                                (Sec. 3)

    theta <- theta + theta_dot * dtau
        # RENDERING CHOICE, not the paper's: Sec. 5.1 states Eqs. 36, 37, 41, 42 as
        #   continuous relations and never discretises tau. Both numerical examples
        #   in the paper are real-time

E(theta) = sum_j h_j <phi(theta)| sigma_j |phi(theta)>       # measured term by term
return  E(theta) read off wherever the trajectory has settled
# the paper states no bound on how much imaginary time a ground state needs; the
#   only tau it fixes is tau = 1/(2T), taking the maximally mixed state to the
#   Gibbs state at temperature T                                        (Sec. 5.3)
# mixed-state imaginary time (two registers plus a controlled-SWAP, Fig. 3) and
#   variational Gibbs-state preparation are a separate extension (Sec. 5.2-5.3),
#   not part of this pure-state hop
# shot and accuracy scaling for each Hadamard test is not derived here; it is
#   deferred to Ref. [41], Li and Benjamin                               (Sec. 6.2)

Cost, as the source states it

Yuan, Endo, Zhao, Li and Benjamin cost one integration step, never a run. For an $N$-parameter ansatz $R_N(\theta_N)\cdots R_1(\theta_1)|\bar{0}\rangle$ and $H=\sum_j h_j\sigma_j$, a linear combination of tensor products of local operators, imaginary-time McLachlan is an $N\times N$ linear system: $\sum_j A^R_{i,j}\dot{\theta}_j=-C^R_i$ pure (Eq. 42), $\sum_j M_{i,j}\dot{\theta}_j=Y_i$ mixed (Eq. 49). They write it as an ODE and never discretise imaginary time, so the per-step unit is ours. Sec. 6.2 requires every entry: $O(N^2)$ matrix entries, $N$ vector entries, each of the latter summing the $M$ Hamiltonian terms (Eqs. 21-22, pure state). Each overlap is one ancilla circuit of at most $N$ register gates (Fig. 1); mixed states instead need Fig. 3, two registers holding $\rho_0$ and a controlled swap. No asymptotic complexity appears; shot and accuracy scaling is deferred to Li and Benjamin. Nothing bounds $\tau$ for a ground state, though $\tau=1/(2T)$ takes the maximally mixed state to the Gibbs state at temperature $T$ on $2n$ qubits (Sec. 5.3). Both numerical runs are real time; there are no imaginary-time numerics.

Implementations

• QuEST simulation of imaginary-time ground states: H2 and LiH

  • Variational ansatz-based quantum simulation of imaginary time evolution

  The paper that introduced this algorithm and first ran it — "We propose a variational algorithm for simulating imaginary time evolution on a hybrid quantum computer" — posted to arXiv in April 2018, eight months before the method's own citation, which in turn opens its imaginary-time section on "the pure state case studied in Ref. [55]", that reference being this paper. It is written by four of Yuan et al.'s five authors together with Sam McArdle and Tyson Jones, and appeared as npj Quantum Information 5, 75 (2019). It starts from the obstacle rather than the principle: "because the imaginary time operator is non-unitary, it is not straightforward to decompose it into a sequence of unitary gates using Trotterization, and thus directly realise it with a quantum circuit." The answer is to move the non-unitarity onto the parameters and leave the circuit unitary. What the paper then asks is not whether that converges in principle but whether it beats a classical optimiser on molecules, so the reported comparison is head-to-head: "We compare the LiH results to those obtained using the VQE, with gradient descent as the classical optimisation routine."

  McLachlan's variational principle $\delta\|(\partial/\partial\tau + H - E_\tau)|\psi(\tau)\rangle\| = 0$ (Eq. 3) gives the linear system $\sum_j A_{ij}\dot{\theta}_j = C_i$ (Eq. 4), integrated by the Euler rule $\vec{\theta}(\tau + \delta\tau) \simeq \vec{\theta}(\tau) + A^{-1}(\tau)\cdot\vec{C}(\tau)\delta\tau$ (Eq. 7), repeated $N_T = \tau_{total}/\delta\tau$ times. $A$ is generally not invertible, and the Supplementary Materials say what was done about it: the system is "generally underdetermined and leaves us unable to invert $A$, we instead update the parameters under Tikhonov regularisation", minimising $\|C - A\dot{\vec{\theta}}\|^2 + \lambda\|\dot{\vec{\theta}}\|^2$ (Eq. 44), with $\lambda$ estimated at each step "by selecting the corner of a 3-point L-curve" and then forced into $\lambda \in [10^{-4}, 10^{-2}]$, each bound for its own reason: above, because "too large a $\lambda$ over-restricts the change in the parameters and was seen to lead to eventual convergence to non-ground states"; below, because "no regularisation ($\lambda = 0$) saw residuals in $A^{-1}$ disrupt the monotonic decrease in energy". That choice was made against alternatives: "Using Tikhonov regularisation affords us a larger time-step than other tested methods, which included LU decomposition, least squares minimisation, singular value decomposition (SVD) and truncated SVD." The two molecules are not run the same way, and the difference is not cosmetic: "While our H2 tests simulate the full experimental routine, our LiH tests are optimised by individually computing each ansatz derivative $\partial|\psi(\vec{\theta})\rangle/\partial\theta_j$", which populates $A$ with $N_p$ evaluations "in contrast to the $N_p^2$ evaluations involved in a full experimental routine" and is available only off a device: "While obtaining this wavefunction from an experiment would require a number of measurements that grows exponentially with the number of qubits, we can access it directly in our numerical simulations." Shot noise is injected analytically, $A_{ij} \sim N\left(a_{ij}, \frac{1}{N_A}(1/16 - a_{ij}^2)\right)$ (Eq. 46); gate error is modelled as a mixing toward the fully mixed state, "We assume an experimental error rate of $10^{-4}$ per gate", and "For our LiH simulations $D \approx 100$, so the effect of gate error is to skew our expected values by $\epsilon \approx 0.99$."

  Both molecules in the minimal STO-3G basis, and both Hamiltonians written out. H2 has four spin-orbitals; the Bravyi-Kitaev four-qubit Hamiltonian of Eq. (40) reduces, because it "only acts off diagonally on qubits 0 and 2", to the two-qubit $H = g_0 I + g_1 Z_0 + g_2 Z_1 + g_3 Z_0 Z_1 + g_4 Y_0 Y_1 + g_5 X_0 X_1$ (Eq. 41), and "In our work, we consider an internuclear distance of $R = 0.75$ Å and hence $g_0 = 0.2252$, $g_1 = 0.3435$, $g_2 = -0.4347$ $g_3 = 0.5716$, $g_4 = 0.0910$, $g_5 = 0.0910$", with a universal ansatz of eight parameters (Fig. 6). LiH contributes twelve spin-orbitals, cut to an active space of eight: the one-electron reduced density matrix comes from "a classically tractable CISD (configuration interaction, single and double excitations) calculation" at 1.45 Å, is rotated to natural molecular orbitals with diagonal $1.99992, 1.96201, 0.03459, 0.00005, 0.00172, 0.00172$ (Eq. 43), the first orbital frozen as doubly occupied and the fourth dropped as empty, and the remainder mapped with Jordan-Wigner. "All of these steps were carried out using OpenFermion [54], an electronic structure package to transform computational chemistry problems into a form that is suitable for investigation using a quantum computer." The result is 181 Hamiltonian terms against a 137-parameter ansatz — "in our LiH simulations, $N_p = 137$ and $N_H = 181$" — the ansatz being "inspired by the low depth circuit ansatz (LDCA)" at depth $M = 3$, "In total there are $(3 \times 5 \times 7) + (4 \times 8) = 137$ parameters" (Fig. 8).

  Everything reported is classical simulation on a named simulator: "We simulate the aforementioned quantum circuits using the Quantum Exact Simulation Toolkit (QuEST) [55]. True ground states are found by diagonalising the considered Hamiltonians with the GNU Scientific Library (GSL) [56]". No quantum processor appears anywhere in the paper; hardware is named only as a prospect, "it can be realised with current and near-term quantum processors". **H2** (Fig. 2), random initial parameters, $\delta\tau = 0.01$, $R = 0.75$ Å: the variational trajectory is plotted against exact imaginary-time evolution with the fidelity between them inset, and "Our simulations were able to converge to the ground state in all trials." **LiH** (Fig. 3), noise-free, $R = 1.45$ Å: two starts, "a small random perturbation (of at most, $\Delta\theta_j = \pi/50$) from the Hartree-Fock state" and "uniformly random parameters", scored by "the fraction of 1280 simulations which, by the given iteration, have converged to within 1 mHartree of the true ground state". Each method ran at its own ceiling — "We use the maximum stable stepsize $\delta\tau$ for each method such that energy monotonically decreases in the first 200 iterations. The stable timestep for imaginary time was 0.225, and for gradient descent it was 0.886." — and "Fig. 3 shows the imaginary time method outperforming gradient descent. It is able to locate the ground state more quickly, and accurately. This advantage is most noticeable for the case of random start states". **Under noise** (Fig. 4), "the fraction of 100 trials which, after 2000 iterations from uniformly random initial parameter states, finished within 1 mHartree of the true ground state of LiH", at a $10^{-4}$ error rate per gate: imaginary time tolerates far fewer shots on $A$ than on the gradient, but the tolerance has a floor — "if the gradient is not known to sufficient accuracy ($N_C < 2 \times 10^4$), the reliability of imaginary time evolution cannot be improved by increasing $N_A$, and can even perform less effectively than gradient descent." The two toy Hamiltonians of Eq. (10) fix the limit the molecules cannot: "While imaginary time avoids all local minima in system A, both methods can become trapped in local minima for the adversarial system B."

• VarQITE Gibbs-state preparation on ibmq_johannesburg

  • Variational Quantum Boltzmann Machines

  This method run on a quantum processor, reached through the branch the method's own paper opens at Sec. 5.3: "Instead, we can input a maximally entangled state $|\Phi\rangle_d = 1/\sqrt{d}\sum_i|ii\rangle_{AE}$ of system $AE$ and evolve the whole system with Hamiltonian $H \otimes I_d$ under imaginary time $\tau$. Then, the state of system $A$ at time $\tau$ will be the thermal state with temperature $T = 1/2\tau$." Zoufal, Lucchi and Woerner needed exactly those thermal states — a quantum Boltzmann machine is trained on $\rho^{Gibbs} = e^{-H_\theta/(k_BT)}/Z$ — and took the algorithm from this method's own literature: "As introduced in [36, 37], an approximate ITE can be implemented on a gate-based quantum computer by using McLachlan's variational principle", references 36 and 37 being McArdle et al. and Yuan et al. The machine-learning framing does not change the machinery — the abstract calls VarQITE "a technique that is typically used for ground state computation" — and what this artefact adds to the record is that the machinery was executed on a device rather than emulated.

  The maximally mixed input is reached by purification rather than by mixed-state evolution. The trial state is $|\psi_0\rangle = V(\omega(0))|0\rangle^{\otimes 2n} = |\phi^+\rangle^{\otimes n}$ over two $n$-qubit subsystems $a$ and $b$, one qubit of each Bell pair in each, so that tracing out $b$ "results in an $n$-dimensional maximally mixed state"; the effective Hamiltonian is $H_{eff} = H^a_\theta + I^b$, VarQITE propagates to $\tau = 1/2(k_BT)$, and $\rho^{Gibbs}_\omega = \mathrm{Tr}_b[|\psi(\omega(\tau))\rangle\langle\psi(\omega(\tau))|]$. Algorithm 1 pairs each integration step with $p$ further solves of the same matrix, $A(\partial_{\theta_i}\dot{\omega}(t)) = \partial_{\theta_i}C - (\partial_{\theta_i}A)\dot{\omega}(t)$ (Eq. 10), so the Boltzmann-machine gradient falls out of the linear system that the imaginary-time step already builds. Regularisation was chosen empirically: Tikhonov and Lasso and "an $\epsilon$-perturbation of the diagonal" were tested with the parameter set by L-curve fitting, "It turns out that all regularization methods perform similarly well", and "The results discussed in this section employ Tikhonov regularization." Every run carries one extra qubit: "All experiments employ an additional qubit $|0\rangle_{add}$ and parameter $\omega_{add}$ to circumvent a potential phase mismatch between the target $|\psi_\tau\rangle$ and the trained state $|\psi(\omega(\tau))\rangle$", which "increases the dimension of $A$ and $C$ by one". On the device, "readout error-mitigation [78-80] is used to obtain the final results run on real quantum hardware" and "the fidelity for the quantum backend evaluations employ state tomography".

  Two Hamiltonians whose Gibbs states are known exactly, and no dataset: $H_1 = 1.0Z$ and $H_2 = 1.0ZZ - 0.2ZI - 0.2IZ + 0.3XI + 0.3IX$, with $\rho_1^{Gibbs}$ and $\rho_2^{Gibbs}$ printed as explicit $2 \times 2$ and $4 \times 4$ density matrices to compare against. The temperature is fixed rather than swept: "The effective temperature, which in principle acts as a scaling factor on the Hamiltonian parameters, is set to $(k_BT) = 1$ in all experiments", so the integration runs to $\tau = 1/2$. The parameterised circuits that carry each of the two states are given in Fig. 2, and the initial parameters are what makes them Bell pairs.

  Fig. 3 is the run: "Fidelity between trained and target Gibbs state with VarQITE for (a) $\rho^{Gibbs}_1$ (b) $\rho^{Gibbs}_2$ trained with an ideal simulator and real quantum hardware, i.e., the ibmq_johannesburg 20-qubit backend. Each simulation used 10 time steps." With readout error mitigation applied and the fidelity read off state tomography, "the method approximates the states, we are interested in, reasonably well and that also the real quantum hardware achieves fidelity values over 0.99 and 0.96, respectively" — the one-qubit Gibbs state first, the two-qubit one second. The device appears at this step only. The Boltzmann-machine training that consumes these states is simulation: the abstract separates the two, "The variational Gibbs state approximation is demonstrated with numerical simulations and experiments run on real quantum hardware provided by IBM Quantum. Furthermore, we illustrate the application of this variational QBM approach to generative and discriminative learning tasks using numerical simulation", and the generative task "runs on an ideal simulation of a quantum computer using AMSGrad [81] with initial learning rate 0.1, maximum number of iterations 200, first momentum 0.7, and second momentum 0.99 as optimization routine".

• VarQITE in Qiskit Algorithms

  The released implementation, and the one a reader can actually run. Qiskit's variational time-evolution module carries this method as a class, `VarQITE`, alongside its real-time sibling `VarQRTE` under a shared `VarQTE` base whose docstring names one reference for both: "[1] Benjamin, Simon C. et al. (2019). Theory of variational quantum simulation.", linked to https://doi.org/10.22331/q-2019-10-07-191 — the Quantum article of this method's own citation, though the docstring's attribution is off, that paper's first author being Xiao Yuan. The library's tutorial states the same lineage in its own words: "In VarQTE, the time evolution of the _state_ $|\Psi(t)\rangle$ is replaced by the evolution of _parameters_ $\theta(t)$ in a variational ansatz", citing "Yuan et al. Quantum 3, 191".

  Construction takes an ansatz circuit, initial parameter values and a variational principle, which "Defaults to ``ImaginaryMcLachlanPrinciple``" - the class that "aims to minimize the distance between both sides of the Wick-rotated Schrödinger equation with a quantum state given as a parametrized trial state. The principle leads to a system of linear equations handled by a linear solver." Its two arguments are the two objects that populate that system, and both have defaults: for the quantum geometric tensor, "If ``None`` provided, ``LinCombQGT`` is used", and for the state gradient, "If ``None`` provided, ``LinCombEstimatorGradient`` is used". The integration and the linear solve are replaceable too: `ode_solver` is `ForwardEulerSolver` unless a SciPy `OdeSolver` or a SciPy method name is passed; for `lse_solver`, "If ``None``, the default ``np.linalg.lstsq`` solver is used"; and for `num_timesteps`, "If ``None``, it is automatically selected to achieve a timestep of approximately 0.01. Only relevant in case of the ``ForwardEulerSolver``." Two tolerances name the numerical facts of the imaginary-time branch directly: `imag_part_tol`, "Allowed value of an imaginary part that can be neglected if no imaginary part is expected", and `num_instability_tol`, "The amount of negative value that is allowed to be rounded up to 0 for quantities that are expected to be non-negative", both `1e-7` by default.

  `VarQITE` lives in `qiskit_algorithms/time_evolvers/variational/var_qite.py` in https://github.com/qiskit-community/qiskit-algorithms, is Python, and is licensed Apache 2.0 ("Copyright 2017 IBM and its contributors"). It is installed as a package rather than vendored — the README says `pip install qiskit-algorithms` — and PyPI carries `qiskit-algorithms` 0.4.0, uploaded 2025-08-29, licence Apache-2.0, requiring Python >= 3.9. The README also carries a warning that belongs next to any recommendation of it: "**Qiskit Algorithms is no longer officially supported by IBM**. Like any other Apache 2 licensed code, you are free to use it or/and extend it, but please be aware that it is under your own risk." The imaginary-time principle is a separate file one directory down, `variational_principles/imaginary_mc_lachlan_principle.py`, and the worked runs are in `docs/tutorials/11_VarQTE.ipynb`.

  The repository's own tutorial runs it and keeps the output, which is what makes this an implementation rather than a listing. The Hamiltonian is a two-spin Ising chain, `SparsePauliOp(["ZZ", "IX", "XI"], coeffs=[-0.2, -1, -1])`, from $H = -J\sum_{i=0}^{L-2}Z_iZ_{i+1} - h\sum_{i=0}^{L-1}X_i$ "Taking $L=2$, $J=0.2$, and $h=1$"; the ansatz is `efficient_su2(reps=1)` with every initial parameter set to $\pi/2$; the imaginary time is `time = 5.0`; and expectation values are taken with Qiskit's `StatevectorEstimator` primitive. The stored cell output reads `Ground state energy -2.0097479079521197`. Re-running the same problem with the classically efficient gradients, `ImaginaryMcLachlanPrinciple(qgt=ReverseQGT(), gradient=ReverseEstimatorGradient())`, gives `Ground state energy -2.0097479079521183`, and "The execution time is reduced by about 33% (this may vary for each execution, but generally results in a speedup)". The reference curve it is checked against is not another quantum run but `SciPyImaginaryEvolver` at 501 time steps, and the notebook's own verdict is "As the above figure indicates, we have obtained the converged ground state energy." This is a documented example with a stored output, not a benchmark sweep: no molecule, no device, no error bars.

What it needs

1. Choose a parameterised trial state 13 methods

   Fix the gate structure of a circuit family and leave its rotation angles open. What comes back is not a circuit but the set of states the later optimisation is allowed to search — which is why this is a slot of its own and not a paragraph in one method's write-up.

Other ways to fill the same slot

Different approaches

• Variational quantum eigensolver

  Prepare a parameterised trial state on the quantum computer, measure the Hamiltonian's expectation value in it, and let a classical optimiser move the parameters. The quantum computer never runs a long coherent evolution; it runs a short one many times, and the loop closes through a classical number.

• Quantum imaginary-time evolution

  Approximate evolution in imaginary time, which damps every excited state faster than the ground state, by a sequence of unitaries determined from measurements on a local domain. There is no variational ansatz and no classical optimiser: the method converges because imaginary time converges.

• Ground-state energy by phase estimation

  Prepare a state that already overlaps the ground state, evolve it under the molecular Hamiltonian, and read the energy off the accumulated phase. Nothing is optimized and nothing is varied — the answer is a measured eigenvalue, and the whole difficulty moves into the starting state.

In the Atlas

• Variational imaginary-time evolution

  McLachlan-style projected imaginary-time dynamics update parameters toward low-energy states.

Sources

• Theory of variational quantum simulation↗Xiao Yuan, Suguru Endo, Qi Zhao, Ying Li, Simon Benjamin · 2018theory yes · simulation yes · hardware no