How to speed up expval measurement?

I need to measure from 4k to 10k expvals in a 12 to 15 qubits circuit. It looks like

def circuit():
  # quantum gates
  measurements: List[qml.Hamiltonian]
  return [qml.expval(measurement) for measurement in measurements]

Unfortunately I cannot use the lightning.qubit because some of the measurements is a linear combination with a SparseHamiltionian.

Is my network too large? Is there a way to speed it up?

I tried setting the shots parameter when declaring the device. That returns the error Measurement expval(0.5 * (X(0) @ X(1) @ X(2) @ X(3) @ ... + SparseHamiltonian ...) not accepted with finite shots on default.qubit

Hi @mchau ,

Can you please help me understand your use case?

lightning.qubit has no problem working with SparseHamiltonian, but SparseHamiltonian does have a problem working with shots as you can see in the docs.

I made a minimal working example using the first example in SparseHamiltonian and it runs with no issues.

import pennylane as qml

# Example from the docs of qml.SparseHamiltonian
# ----------------
wires = range(10)
coeffs = [1 for _ in wires]
observables = [qml.Z(i) for i in wires]
H = qml.Hamiltonian(coeffs, observables)
Hmat = H.sparse_matrix()
H_sparse = qml.SparseHamiltonian(Hmat, wires)
# ----------------

# Now create a list of 4k measurements 
measurements = [i*qml.X(0)+ H_sparse for i in range(4000)]

# Use lightning.qubit
dev = qml.device("lightning.qubit",wires=10)
@qml.qnode(dev)
def circuit():
  # Add the measurements
  return [qml.expval(measurement) for measurement in measurements]

# Run the circuit
circuit()
1 Like

Hi @CatalinaAlbornoz , thank you for the effort. Here is something closer to what I am trying to do. It would take way longer than your example

from scipy.sparse import csr_matrix
import numpy as np
import pennylane as qml
from pennylane import X, Y, Z, I

n = 2**12
rows = np.arange(n)
cols = n - 1 - rows
data = np.ones(n)

sparse_matrix = csr_matrix((data, (rows, cols)), shape=(n, n))

wires = range(12)
coeffs = [1 for _ in wires]

test_ops = X(0) @ X(1) @ X(2) @ X(3) @ X(4) @ X(5) @ X(6) @ I(7) @ Z(8) @ Z(9) @ Y(10) @ Y(11)
H_sparse = qml.SparseHamiltonian(sparse_matrix, wires)
observables = [test_ops + H_sparse @ test_ops @ qml.adjoint(H_sparse)]

dev = qml.device("lightning.qubit", wires=12)
@qml.qnode(dev)
def circuit():
  return [qml.expval(i*observables[0]) for i in range(4000)]

circuit()

As far as I can see the problem is I am measuring test_ops + H_sparse @ test_ops @ qml.adjoint(H_sparse) so it takes long

Hi @mchau,

Remember that every additional qubit that you add is exponentially harder and slower to simulate. So using 12 qubits is going to be much slower than using 10. I ran your code for 12 qubits and 4 measurements and it took a few minutes. I then ran it for 10 qubits and 4 measurements and it took a couple of seconds. The difference is massive. That’s the motivation for using actual quantum computers in the future.

The easiest change you can make here is reduce the number of qubits. If you cannot, then maybe you can try to use a Python profiler (e.g. SnakeViz) to see where everything is taking so long.

You could potentially try using Circuit Cutting but I don’t know whether it will work or not. Make sure to check out the documentation here and here to learn how to use this feature.

I hope this helps.

1 Like

I see, thank you. I two follow up questions

In the code, we have

return [qml.expval(i*observables[0]) for i in range(4000)]
  1. I understand that after the 1st measurement, the state collapses so should the subsequent measurement should be fast enough? Hence the expval should be faster too?
  2. If 1. is wrong, does it violate any QM laws to parallel compute (python’s ProcessPoolExecutor for example) to speed this up. I noticed we are only using 1 CPU core at this expvals

Hi!

  1. No, it doesn’t work that way. Take this simpler example where I try to find three different expectation values.
import pennylane as qml

dev = qml.device("default.qubit", wires=2)

@qml.qnode(dev)
def circuit(x):
    qml.RX(x, wires=0)
    qml.Hadamard(wires=0)
    qml.CNOT(wires=[0, 1])
    return qml.expval(qml.Y(1)),qml.expval(qml.Y(0)),qml.expval(qml.Y(0)@qml.Y(1))

print(circuit(0.5))

The output of qml.expval(qml.Y(0)) is independent of the fact that I had measured qml.expval(qml.Y(1)).

  1. I would like to know more about your project. It would help me in providing better answers. What is the intention of this line return [qml.expval(i*observables[0]) for i in range(4000)]? It looks as though you were calculating the same expectation 4000 times but multiplying it by a scalar each time. Is that so? if so, why is it done this way?

The way I see it is that you are trying to find a nontrivial expectation value for a problem with high dimensionality and the only way to make it faster will be by reducing the amount of qubits. However, if you explain a little more, I might be able to provide more help.