Shadow Hamiltonian simulation

Shadow Hamiltonian simulation — Definition

In common quantum Hamiltonian simulation, we evolve a state vector \(|\psi(t)\rangle\) according to the Schrödinger equation,

by some Hamiltonian \(H,\) and then compute expectation values of the evolved state through measurement. In shadow Hamiltonian simulation, we encode a set of expectation values in the amplitudes of a quantum state, and evolve those according to some shadow Schrödinger equation.

For that, we first need to define the shadow state,

for a set of operators \(S = \{O_m\}\) and normalization constant \(A = \sum_m |\langle O_m \rangle|^2.\) This means that we can encode these \(M\) operator expectation values into \(n_S\) qubits, as long as \(2^{n_S} \geq M.\) Note that \(\langle O_m \rangle = \text{tr}[O_m \rho],\) so we can have mixed or even infinite-dimensional states \(\rho.\)

The shadow state evolves according to its shadow Schrödinger equation,

The Hamiltonian matrix \(H_S\) is given by the commutation relations between the system Hamiltonian \(H\) and the operators in \(S = \{O_m\},\)

Let us solve for the matrix elements \((H_S)_{m m'}.\) To do this, recall that a vector \(\boldsymbol{v}\) can always be decomposed in an orthogonal basis \(\boldsymbol{e}_j\) via \(\boldsymbol{v} = \sum_j \frac{\langle \boldsymbol{e}_j, \boldsymbol{v}\rangle}{||\boldsymbol{e}_j||^2} \boldsymbol{e}_j.\) Since the operators under consideration are elements of the vector space of Hermitian operators, we can use this to compute \(H_S.\)

In particular, with the trace inner product, this amounts to

from which we can read off the matrix elements of \(H_S,\) i.e.,

Now, we can see that the operators \(O_m\) need to be chosen such that all potentially new operators \(\mathcal{O} = [H, O_m]\), resulting from taking the commutator between \(H\) and \(O_m,\) are decomposable in terms of \(O_m\) again. In particular, the operators \(O_m\) need to form a basis for \(\{\mathcal{O} | \mathcal{O} = [H, O_m] \}.\) Another way to say this is that \(\{O_m\}\) need to contain all nested commutators \([[[H, O_m], O_m'], .. ],\) which is similar to lie_closure() but weaker because it revolves around just \(H.\) In the paper this is called the invariance property.

How this relates to g-sim

In g-sim 2 3 4 5, we have operators \(\{ g_i \}\) that are generators or observables for a parametrized quantum circuit, e.g. \(U(\theta) = \prod_\ell \exp(-i \theta_\ell g_\ell)\) and \(\langle g_i \rangle.\) For that, we are looking at the so-called dynamical Lie algebra (DLA) of the circuit, \(\mathfrak{g} = \langle \{ g_i \} \rangle_\text{Lie} = \{ g_1, .., g_{|\mathfrak{g}|} \},\) as well as the adjoint representation \((-i \text{ad}_{g_\gamma})_{\alpha \beta} = f^\gamma_{\alpha \beta},\) where \(f^\gamma_{\alpha \beta}\) are the structure_constants() of the DLA. They are computed via

The operators in \(\frak{g}\) can always be orthonormalized via the Gram–Schmidt process, in which case we can drop the denominator. Further, by means of the cyclic property of the trace, we can rewrite this expression to obtain

From this, we see how \(H_S\) corresponds to the adjoint representation \(i \text{ad}_H\) (but we don’t require the full Lie algebra here, see below). For further details on the concept of the adjoint representation, see our demo on g-sim that makes extensive use of it.

In g-sim, we also evolve expectation vectors \((\vec{g})_i = \langle g_i \rangle.\) In particular, the circuit of evolving a state according to \(U(\theta)\) and computing expectation values \(\langle g_i \rangle\) then corresponds to evolving \(\vec{g}\) by \(\prod_\ell \exp(-i \theta_\ell \text{ad}_{g_\ell}).\)

Shadow Hamiltonian simulation can thus be viewed as g-sim with a single, specific gate \(U(\theta) = e^{-i \theta H}\) and parameter \(\theta = t,\) and run on a quantum computer.

One striking difference is that, because we only have one specific “gate”, we do not need the full Lie closure of the operators whose expectation values we want to compute. Instead, here it is sufficient to choose \(O_m\) such that they build up the full support for all \([H, O_m].\) This could potentially be a significant difference, as the Lie closure in most cases leads to an exponentially large DLA 6 7, though the scaling of the span of all \([H, O_m]\) is unclear at this point.

A simple example

The abstract concepts of shadow Hamiltonian simulation are best illustrated with a simple and concrete example. We are interested in simulating the Hamiltonian evolution of

after a time \(t = 1\) and computing the expectation values of \(S = \{X, Y, Z, I \}.\) In the standard formulation, we simply evolve the initial quantum state \(|\psi(0)\rangle = |0\rangle\) by \(H\) in the following way.

import pennylane as qml
import numpy as np
from pennylane import X, Y, Z, I

dev = qml.device("default.qubit")

S = [X(0), Y(0), Z(0), I(0)]
H = X(0) + Y(0)

@qml.qnode(dev)
def evolve(H, t):
    qml.evolve(H, t)
    return [qml.expval(Om) for Om in S]

t = 1.
O_t_standard = np.array(evolve(H, t))
O_t_standard

array([ 0.21783962, -0.21783962, -0.95136313,  1.        ])

We evolved a \(2^n = 2\) dimensional quantum state and performed \(3\) independent (non-commuting) measurements.

In shadow Hamiltonian simulation, we encode \(4\) expectation values in a \(2^2 = 4\)-dimensional quantum state, i.e., \(n_S = 2.\)

For this specific example, the number of operators is larger than the number of qubits, leading to a shadow system that is larger than the original system. This may or may not be a clever choice, but the point here is just to illustrate the conceptual difference between both approaches. The authors in 1 show various examples where the resulting shadow system is significantly smaller than the original system. It should also be noted that having a smaller shadow system may not always be its sole purpose, as there are conceptually new avenues one can explore with shadow Hamiltonian simulation, such as sampling from the distribution \(p_m = |\langle O_m \rangle |^2.\)

Let us first construct the initial shadow state \(\boldsymbol{O}(t=0)\) by computing \(\langle O_m \rangle_{t=0} = \text{tr}\left(O_m |\psi(0)\rangle \langle \psi(0)| \right)\) with \(|\psi(0)\rangle = |0\rangle.\) The pauli_rep attribute of PennyLane operators returns a PauliSentence instance and lets us efficiently compute the trace, where we use the trick that \(|0 \rangle \langle 0| = (I + Z)/2.\)

S_pauli = [op.pauli_rep for op in S]

O_0 = np.zeros(len(S))

for m, Om in enumerate(S_pauli):
    psi0 = (I(0) + Z(0)).pauli_rep

    O_0[m] = (psi0 @ Om).trace()


O_0

array([0., 0., 1., 1.])

There are a variety of methods to encode this vector in a qubit basis, but we will just be using StatePrep later.

We now go on to construct the shadow Hamiltonian \(H_S\) by computing the elements \((H_S)_{m m'} = \frac{\text{tr}\left( O_{m'} [H, O_m] \right)}{|| O_{m'} ||^2},\) and we again make use of the trace() method.

H_pauli = H.pauli_rep

H_S = np.zeros((len(S), len(S)), dtype=complex)

for m, Om in enumerate(S_pauli):
    com = H_pauli.commutator(Om)
    for mt, Omt in enumerate(S_pauli):
        # v = ∑ (v · e_j / ||e_j||^2) * e_j

        value = (Omt @ com).trace()
        value = value / (Omt @ Omt).trace()
        H_S[m,mt] = value

H_S = -H_S # definition eq. (2) in [1]

In order for the shadow evolution to be unitary and implementable on a quantum computer, we need \(H_S\) to be Hermitian.

np.all(H_S == H_S.conj().T)

True

Knowing that, we can write the formal solution to the shadow Schrödinger equation as

from scipy.linalg import expm

O_t = expm(-1j * t * H_S) @ O_0
O_t

array([ 0.21783962+0.j, -0.21783962+0.j, -0.95136313+0.j,  1.        +0.j])

Up to this point, this is equivalent to g-sim if we were doing classical simulation. Now, the main novelty for shadow Hamiltonian simulation is to perform this on a quantum computer by encoding the expectation values of \(\langle O_m \rangle\) in the amplitude of a quantum state, and to translate \(H_S\) accordingly.

This can be done by decomposing the numerical matrix \(H_S\) into Pauli operators, which can, in turn, be implemented on a quantum computer.

H_S_qubit = qml.pauli_decompose(H_S)
H_S_qubit

(
    -1.0 * (X(0) @ Y(1))
  + -1.0 * (Y(0) @ I(1))
  + 1.0 * (Y(0) @ X(1))
  + -1.0 * (Y(0) @ Z(1))
)

Using all these ingredients, we now are able to formulate the shadow Hamiltonian simulation as a quantum algorithm. For the amplitude encoding, we need to make sure that the state is normalized. We use that normalization factor to then later retrieve the correct result.

A = np.linalg.norm(O_0)

@qml.qnode(dev)
def shadow_evolve(H_S_qubit, O_0, t):
    qml.StatePrep(O_0 / A, wires=range(2))
    qml.evolve(H_S_qubit, t)
    return qml.state()

O_t_shadow = shadow_evolve(H_S_qubit, O_0, t) * A

print(O_t_standard)
print(O_t_shadow)

[ 0.21783962 -0.21783962 -0.95136313  1.        ]
[ 0.21783962+0.j -0.21783962+0.j -0.95136313+0.j  1.        +0.j]

We see that the results of both approaches match.

The first result is coming from three independent measurements on a quantum computer after evolution with system Hamiltonian \(H.\) This is conceptually very different from the second result where \(\boldsymbol{O}\) is encoded in the state of the shadow system (note the qml.state() return), which we evolved according to \(H_S.\)

In the first case, the measurement is directly obtained, however, in the shadow Hamiltonian simulation, we need to access the amplitudes of the underlying state. This can be done naively with state tomography, but in instances where we know that \(\langle O_m \rangle \geq 0,\) we can just sample bitstrings according to \(p_m = |\langle O_m\rangle|^2.\) The ability to sample from such a distribution \(p_m = |\langle O_m\rangle|^2\) is a unique and new feature to shadow Hamiltonian simulation.

We should also note that we made use of the abstract quantum sub-routines evolve() and StatePrep, which each warrant their specific implementation. For example, StatePrep can be realized by MottonenStatePreparation and evolve() can be realized by TrotterProduct, though that is not be the focus of this demo.

Conclusion

We introduced the basic concepts of shadow Hamiltonian simulation and learned how it fundamentally differs from the common approach to Hamiltonian simulation.

We have seen how classical Hamiltonian simulation is tightly connected to g-sim, but run on a quantum computer. A significant difference comes from the fact that the authors in 1 specifically look at Hamiltonian simulation, \(\exp(-i t H),\) which allows us to just look at operators \(O_m\) that support all commutators \([H, O_m],\) instead of the full Lie closure. There may be some advantage to this feat, because Lie algebras in quantum computing typically scale exponentially 6 7. However, the scaling of such sets of operators is unclear at this point and needs further investigation.

Note that even in the case of an exponentially sized set of operators, we have — at least in principle — an exponentially large state vector to store the \(M \leq 2^{n_S}\) values. In the absolute worst case we have \(\mathfrak{su}(2^n)\) with a dimension of \(2^{2n}-1,\) so \(n_S = 2n\) and thus it is just doubling the number of qubits.

The biggest potential to this new persepctive on Hamiltonian simulation most likely lies in finding interesting applications like 8 or 9 that naturally encode the problem and allow for efficient retrieval of all the relevant information.

References

1(1,2,3,4)

Rolando D. Somma, Robbie King, Robin Kothari, Thomas O’Brien, Ryan Babbush “Shadow Hamiltonian Simulation” arXiv:2407.21775, 2024.

2

Rolando D. Somma “Quantum Computation, Complexity, and Many-Body Physics” arXiv:quant-ph/0512209, 2005.

3

Rolando Somma, Howard Barnum, Gerardo Ortiz, Emanuel Knill “Efficient solvability of Hamiltonians and limits on the power of some quantum computational models” arXiv:quant-ph/0601030, 2006.

4

Victor Galitski “Quantum-to-Classical Correspondence and Hubbard-Stratonovich Dynamical Systems, a Lie-Algebraic Approach” arXiv:1012.2873, 2010.

5

Matthew L. Goh, Martin Larocca, Lukasz Cincio, M. Cerezo, Frédéric Sauvage “Lie-algebraic classical simulations for variational quantum computing” arXiv:2308.01432, 2023.

6(1,2)

Roeland Wiersema, Efekan Kökcü, Alexander F. Kemper, Bojko N. Bakalov “Classification of dynamical Lie algebras for translation-invariant 2-local spin systems in one dimension” arXiv:2309.05690, 2023.

7(1,2)

Gerard Aguilar, Simon Cichy, Jens Eisert, Lennart Bittel “Full classification of Pauli Lie algebras” arXiv:2408.00081, 2024.

8

Ryan Babbush, Dominic W. Berry, Robin Kothari, Rolando D. Somma, Nathan Wiebe “Exponential quantum speedup in simulating coupled classical oscillators” arXiv:2303.13012, 2023.

9

Alice Barthe, M. Cerezo, Andrew T. Sornborger, Martin Larocca, Diego García-Martín “Gate-based quantum simulation of Gaussian bosonic circuits on exponentially many modes” arXiv:2407.06290, 2024.