(g + P)-sim: Extending g-sim by non-DLA observables and gates

Introduction

Lie-algebraic simulation techniques such as \(\mathfrak{g}\)-sim can be handy in the niche cases where the dynamical Lie algebra (DLA) scales polynomially with the number of qubits. Because those cases essentially boil down to the transverse field Ising model (TFIM) and variants thereof in 1D 5, we will do a case study on its DLA specifically.

We are interested in the case where we want to extend the DLA \(\mathfrak{g}\) by a few additional gates that are outside the DLA. For \(n\) qubits we get a DLA dimension of \(\text{dim}(\mathfrak{g}) = 2n(2n-1)/2\) for the TFIM (see here). Suppose we want to expand the DLA by a single operator \(p\) in order to use it as a gate, and let us assume that \(p\) is the product of two DLA operators that, itself, is not part of the DLA. Adding product operators to the TFIM DLA and computing their new Lie closure can lead to an exponential increase with a new dimension up to \(2(2^{2n-2}-1).\) In that worst case, we get the so-called associative algebra of \(\mathfrak{g};\) that is, the algebra from the closure over multiplication, i.e. which looks at all possible products of operators. This is also a Lie algebra.

Here, we show how to extend the DLA by such a \(p\) gate without going to the exponentially large associative algebra, but instead make use of the fact that \(p\) is a product of DLA elements. We do so by looking at moments of \(\mathfrak{g}\) instead. The \(m\)-th order moments are products of \((m+1)\) DLA elements. E.g. \(p = h_{\alpha_1} h_{\alpha_2} \notin \mathfrak{g}\) is a first order moment. Depending on their order, every non-DLA moment gate increases the highest moment order considered in the computation, \(m_\text{comp}\). The overall cost scales with the maximum order \(\text{dim}(\mathfrak{g})^{m_\text{comp}}.\)

In the worst case, each moment expands the space of operators by a factor \(\text{dim}(\mathfrak{g}),\) such that for \(m\) moments, we are dealing with a \(\text{dim}(\mathfrak{g})^{m+2}\) dimensional space. In that sense, this is similar to Clifford+T simulators where expensive \(T\) gates come with an exponential cost. A key difference is that for a finite dimensional DLA, there is a maximum moment \(m_\text{max}.\) This corresponds to simply constructing the full associative algebra again. In the case that the required \(m_\text{comp} = m_\text{max},\) we can just perform regular \(\mathfrak{g}\)-sim with the associative algebra. Here we will consider the case \(m_\text{comp} < m_\text{max}.\)

In 4, the authors already hint at the possibility of extending \(\mathfrak{g}\)-sim by expectation values of products of DLA elements. In this demo, we extend this notion to gates generated by moments of the DLA.

-sim\(\mathfrak{g}\)-sim

Let us briefly recap the core principles of \(\mathfrak{g}\)-sim. We consider a Lie algebra \(\mathfrak{g} = \{h_1, .., h_d\},\) which is closed under commutation (see lie_closure()). We know that gates \(e^{-i \theta h_\alpha}\) transform Lie algebra elements into Lie algebra elements,

This is the adjoint identity with the adjoint representation of the Lie algebra given by the structure_constants(), \(f^\mu_{\alpha \beta} = -i \left(\text{ad}_{h_\mu}\right)_{\alpha \beta}.\)

This lets us evolve any expectation value of DLA elements using the adjoint representation of the DLA. For that, we define the expectation value vector \((\boldsymbol{e})_\alpha = \text{tr}[h_\alpha \rho].\)

Also, let us write \(U = e^{-i \theta \text{ad}_{h_\mu}}\) corresponding to a unitary \(\mathcal{U} = e^{-i \theta h_\mu}.\) Using the adjoint identity above and the cyclic property of the trace, we can write an evolved expectation value vector as

Hence, the expectation value vector is simply transformed by matrix multiplication with \(U\) and we have

for some input expectation value vector \(\boldsymbol{e}^\text{in}.\)

A circuit comprised of multiple unitaries \(\mathcal{U}\) then simply corresponds to evolving the expectation value vector with \(U.\)

We are going to concretely use the DLA of the transverse field Ising model,

This is one of the few systems that yield a polynomially sized DLA. Let us construct its DLA via lie_closure().

import pennylane as qml
import numpy as np

from pennylane import X, Y, Z, I
from pennylane.pauli import PauliSentence, PauliWord
from scipy.linalg import expm

import copy

# TFIM generators
def TFIM(n):
    generators = [X(i) @ X(i+1) for i in range(n-1)]
    generators += [Z(i) for i in range(n)]
    generators = [op.pauli_rep for op in generators]

    dla = qml.lie_closure(generators, pauli=True)
    dim_g = len(dla)
    return generators, dla, dim_g

generators, dla, dim_g = TFIM(n=4)

In regular \(\mathfrak{g}\)-sim, the unitary evolution \(\mathcal{U}\) of the expectation vector is simply generated by the adjoint representation \(U.\)

adjoint_repr = qml.structure_constants(dla)
gate = adjoint_repr[-1]
theta = 0.5

U = expm(theta * gate)

-sim\((\mathfrak{g}+P)\)-sim

We now want to extend \(\mathfrak{g}\)-sim by operators that are not in the DLA, but a product of DLA operators. Note that while the DLA is closed under commutation, it is not closed under multiplication, such that products of DLA elements are in general not in \(\mathfrak{g}.\)

Let us look at the adjoint action of a gate generated by \(p = h_{\mu_1} h_{\mu_2} .. \notin \mathfrak{g},\)

Here, \(\boldsymbol{P}^m\) correspond to the contributions of the \(m\)-th moments in \(\mathfrak{g}.\) Let us look at the case where we use a first order product, \(\mathcal{P} = e^{-i \theta h_{\mu_1} h_{\mu_2}},\) and only DLA operators and first order moments contribute to the adjoint action.

For that, let us construct a concrete example. First we pick two elements from \(\mathfrak{g}\) such that their product is not in \(\mathfrak{g}.\)

p = dla[-5] @ dla[-1]
p = next(iter(p)) # strip any scalar coefficients
dla_vspace = qml.pauli.PauliVSpace(dla, dtype=complex)
dla_vspace.is_independent(p.pauli_rep)

True

Now, we compute \(e^{i \theta h_{\mu_2} h_{\mu_1}} h_\alpha e^{-i \theta h_{\mu_1} h_{\mu_2}}\) and decompose it in the DLA and first moments.

Note that since the product basis is overcomplete, we only keep track of the linearly independent elements and ignore the rest.

def exppw(theta, ps):
    # assert that it is indeed a pure Pauli word, not a sentence
    assert (len(ps) == 1 and isinstance(ps, PauliSentence)) or isinstance(ps, PauliWord)
    return np.cos(theta) * PauliWord({}) + 1j * np.sin(theta) * ps

theta = 0.5

P = exppw(theta, p)
P_dagger = exppw(-theta, p) # complex conjugate with p just being a hermitian pauli word

P0 = np.zeros((dim_g, dim_g), dtype=float)

for i, h1 in enumerate(dla):
    res = P @ h1 @ P_dagger
    for j, h2 in enumerate(dla):
        # decompose the result in terms of DLA elements
        # res = ∑ (res · h_j / ||h_j||^2) * h_j
        value = (res @ h2).trace().real
        value = value / (h2 @ h2).trace()
        P0[i, j] = value

P1 = np.zeros((dim_g, dim_g, dim_g), dtype=float)

for i, h1 in enumerate(dla):
    res = P @ h1 @ P_dagger
    dla_and_M1_vspace = copy.deepcopy(dla_vspace)
    for j, h2 in enumerate(dla):
        for l, h3 in enumerate(dla):
            prod = h2 @ h3

            if not dla_and_M1_vspace.is_independent(prod):
                continue

            # decompose the result in terms of products of DLA elements
            # res = ∑ (res · p_j / ||p_j||^2) * p_j
            value = (res @ prod).trace().real
            value = value / (prod @ prod).trace().real
            P1[i, j, l] = value
            dla_and_M1_vspace.add(prod)

We want to confirm that the adjoint action of \(\mathcal{P}\) is indeed fully described by the zeroth and first moments.

For that, we reconstruct the transformed DLA elements and compare them with the decomposition.

for i, h1 in enumerate(dla):
    res = P @ h1 @ P_dagger
    res.simplify()

    reconstruct = sum([P0[i, j] * dla[j] for j in range(dim_g)])
    reconstruct += sum([P1[i, j, l] * dla[j] @ dla[l] for j in range(dim_g) for l in range(dim_g)])
    reconstruct.simplify()

    assert res == reconstruct

Now that we have successfully constructed a \(\mathcal{P}\) gate, let us look how entering it in a circuit transforms DLA elements (and therefore expectation value vector elements).

Here we have defined the expectation tensor \((\boldsymbol{E}^m)_{\beta_1 , .. , \beta_{m+1}} = \text{tr}\left[ h_{\beta_1} .. h_{\beta_{m+1}} \rho \right]\) for the \(m\)-th moment. Note that \(\boldsymbol{e} = \boldsymbol{E}^0\) is the expectation value vector for regular \(\mathfrak{g}\)-sim.

Such a computation corresponds to the branching off from the original diagram, with an extra contribution coming from the higher moments.

When inserting an arbitrary DLA gate \(U\) before and \(V\) after the \(\mathcal{P}\) gate, we obtain the following diagram. Note that \(U\) and \(V\) can be compositions of multiple DLA gates again.

Note that in one vertical column the \(U\) correspond to the same matrix.

Example

Let us compute an example. For that we start by computing the initial expectation vector and tensor.

E0_in = np.zeros((dim_g), dtype=float)
E1_in = np.zeros((dim_g, dim_g), dtype=float)
E_in = [E0_in, E1_in]

for i, hi in enumerate(dla):
    rho_in = qml.prod(*(I(i) + Z(i) for i in hi.wires))
    rho_in = rho_in.pauli_rep

    E_in[0][i] = (hi @ rho_in).trace()

for i, hi in enumerate(dla):
    for j, hj in enumerate(dla):
        prod = hi @ hj
        if prod.wires != qml.wires.Wires([]):
            rho_in = qml.prod(*(I(i) + Z(i) for i in prod.wires))
        else:
            rho_in = PauliWord({}) # identity
        rho_in = rho_in.pauli_rep

        E_in[1][i, j] = (prod @ rho_in).trace().real

Now we need to compute the two branches from the diagram above.

# some other DLA gate V
V = expm(0.5 * adjoint_repr[-2])

# contract first branch
# V - P - U - E^0_in
res0 = U @ E_in[0]
res0 = P0 @ res0
res0 = V @ res0

# contract second branch

# --U-==-+--------+   -+------+
#        | E^1_in | =  | res  |
# --U-==-+--------+   -+------+
res = np.einsum("ij,jl->il", U, E_in[1])
res = np.einsum("kl,il->ik", U, res)

#        +-----+-==-+------+
#  -V-==-| P^1 |    | res  |
#        +-----+-==-+------+
res = np.einsum("ijl,jl->i", P1, res)
res = V @ res

res = res + res0

As a sanity check, let us compare this to the same circuit but using our default state vector simulator in PennyLane.

@qml.qnode(qml.device("default.qubit"))
def true():
    qml.exp(-1j * theta * dla[-1].operation())
    qml.exp(-1j * 0.5 * p.operation())
    qml.exp(-1j * 0.5 * dla[-2].operation())
    return [qml.expval(op.operation()) for op in dla]

true_res = np.array(true())

np.allclose(true_res, res)

True

We find that indeed both results coincide and expectation value vectors are correctly transformed in \((\mathfrak{g}+P)\)-sim.

Higher moments

We can extend the above approach by a second \(P\) gate in the circuit. This leads to contributions from up to the third order. The diagram for a circuit \(P V P U\) has the following five branches.

First, the zeroth- and first-order contribution. This can be seen as the branching off from the first previous diagram.

We also obtain the third-order diagram containing both \(\boldsymbol{P}^1\) tensors.

To get the correct results, we also obtain intermediate second-order diagrams.

Moment vector space

The above diagrams are handy to understand the maximum moment order required for adding \(P\) gates of a certain order. There is, however, a lot of redundancy due to the overcompleteness of naively looking at moments as all possible products between DLA elements.

Instead, we can also work in the vector space of unique moments

that is iteratively built from \(\mathcal{M}^0 = \mathfrak{g}.\)

Even though these spaces in general do not form Lie algebras, we can still compute their (pseudo) adjoint representations and use them for \(\mathfrak{g}\)-sim as long as we work in a large enough space with the correct maximum moment order.

We now set up the moment vector spaces starting from the DLA and keep adding linearly independent product operators.

def Moment_step(ops, dla):
    MomentX = qml.pauli.PauliVSpace(ops.copy())
    for i, op1 in enumerate(dla):
        for op2 in ops[i+1:]:
            prod = op1 @ op2
            # ignore scalar coefficient
            pw = next(iter(prod.keys()))

            MomentX.add(pw)

    return MomentX.basis

Moment0 = dla.copy()
Moment = [Moment0]
dim = [len(Moment0)]
for i in range(1, 5):
    Moment.append(Moment_step(Moment[-1], dla))
    dim.append(len(Moment[-1]))

dim

[28, 98, 126, 127, 127]

We see that for the considered example of \(n = 4\) we reach the maximum moment already for the second order (the additional operator in the third moment space is just the identity).

We can repeat our original computation for the first moments using the \(98\)-dimensional first-order moment vector space \(\mathcal{M}^1.\)

The recipe follows the exact same steps as \(\mathfrak{g}\)-sim but using \(\mathcal{M}^1\) instead. First, we compute the input expectation value vector.

comp_moment = 1

e_in = np.zeros((dim[comp_moment]), dtype=float)

for i, hi in enumerate(Moment[comp_moment]):
    rho_in = qml.prod(*(I(i) + Z(i) for i in hi.wires))
    rho_in = rho_in.pauli_rep

    e_in[i] = (hi @ rho_in).trace()

Next, we compute the (pseudo-)adjoint representation of \(\mathcal{M}^1.\)

adjoint_repr = qml.structure_constants(Moment[comp_moment])

We can now choose arbitrary DLA gates and a maximum of one \(P\) gate to evolve the expectation value vector.

P_index = Moment[comp_moment].index(1.*p) # pick the gate in the Moment space that p corresponds to

e_t = e_in
e_t = expm(0.5 * adjoint_repr[dim_g-1]) @ e_t # the same U gate
e_t = expm(0.5 * adjoint_repr[P_index]) @ e_t      # the same P gate
e_t = expm(0.5 * adjoint_repr[dim_g-2]) @ e_t # the same V gate

The final result matches the state vector result again

np.allclose(e_t[:dim_g], true_res)

True

Limitations

We saw how we can make use of moment vector spaces to extend \(\frak{g}\)-sim by non-DLA elements under certain conditions. The upside is that while the Lie closure or construction of the associative algebra leads to an exponential DLA in \(n,\) we get away with a polynomial cost in \(n,\) as we have \(O(\text{dim}(\mathfrak{g})^{m_{\text{comp}}})\)-sized objects for some fixed maximum moment order \(m_{\text{comp}}\) in the computation, with some additional reductions due to the redundancies in the moment spaces.

However, we argue that while this is interesting in theory, there is little practical utility. To show that, we plot the dimensions of the first and second moments against the associative algebra dimension.

dims_dla = []
dims_moment = []
dims_tempdla = []

ns = np.arange(2, 6)

for n in ns:
    _, dla, dim_g = TFIM(n)

    Moment0 = dla.copy()
    Moment = [Moment0]
    dim = [len(Moment0)]
    for i in range(1, 5):
        Moment.append(Moment_step(Moment[-1], dla))
        dim.append(len(Moment[-1]))

    tempdla = qml.lie_closure(dla + [Moment[1][-1]], pauli=True)

    dims_dla.append(dim_g)
    dims_moment.append(dim)
    dims_tempdla.append(len(tempdla))

import matplotlib.pyplot as plt

plt.title("Dimensions of $\\langle g + P \\rangle_{{Lie}}$ vs $\mathcal{M}^m$")

plt.plot(ns, dims_tempdla, "o--", label="${{dim}}(\\langle g + P \\rangle_{{Lie}})$", color="tab:blue")
plt.plot(ns, 2 * (2**(2*ns - 2) - 1), "-", label="$2(2^{{2n-2}}-1)$", color="tab:blue")

color = ["tab:orange", "tab:green", "tab:cyan"]
dims_moment = np.array(dims_moment)
for i in range(3):
    plt.plot(ns, dims_moment[:, i], "x--", label=f"$dim(\mathcal{{M}}^{i})$", color=color[i])
    fitcoeff = np.polyfit(ns, dims_moment[:, i], i+2) # polynomial fit of order m+1
    plt.plot(ns, np.poly1d(fitcoeff)(ns), "-", label=f"$O(n^{{{i}+2}})$", color=color[i])

plt.xticks(ns)
plt.yscale("log")

plt.legend()
plt.xlabel("n")
plt.show()

/home/runner/work/qml/qml/_build/demonstrations/tutorial_liesim_extension.py:474: RankWarning: Polyfit may be poorly conditioned
  fitcoeff = np.polyfit(ns, dims_moment[:, i], i+2) # polynomial fit of order m+1

First, we note that the maximum moment is quickly reached for small system sizes. In that case we have \(m_\text{comp} = m_\text{max}\) so we might as well look at the associative algebra and perform regular \(\mathfrak{g}\)-sim. Secondly, we also note that the dimensions quickly explode and become hard to handle in reasonable times.

For example, in all cases here there are no non-trivial third-order moments. We would have to go to \(n = 6,\) for which there are \(1980\) elements in \(\mathcal{M}^3,\) which corresponds to iterating over \(1980^3 \approx 2^{32}\) commutators to compute the (pseudo-)adjoint representation. This is already a stretch to accomplish with the available tools.

Hence, this method is effectively restricted to very few non-DLA gates of Ising-type DLAs, rendering the method overall niche for practical applications. On the other hand, we gained some new theoretical insights into the relationship between the simulability of a quantum circuit and its DLA. In particular, we constructed a polynomial algorithm in the number of qubits \(n,\) for simulating DLA circuits with up to \(m_\text{comp}\) moments in \(\mathcal{O}(\text{dim}\left(\mathfrak{g}(n)\right)^{m_\text{comp} + 2}) = \mathcal{O}(\text{poly}(n)).\)

References

1

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

2

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.

3

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

4(1,2)

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.

5

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.