Sampling solution vectors from a circuit

Hi,
It is possible to get the empirical distribution of computational basis states using probs() function. Alternatively, one of Pennylane examples (can’t find it anymore) used a custom observable for sampling

basis_obs = qml.Hermitian(np.diag(range(2 ** num_qubits)), wires=range(num_qubits))
qml.sample(basis_obs)

However, both methods are not very scalable because they require initialization of an exponentially large data structure. For example, the custom Hermitian requires a construction of the numpy matrix of size 2**num_qubits by 2**num_qubits.

Is there a better (preferably quick) way to sample basis states and associated energies without creating exponentially large data structures?
Ideally, I would like to get a set of samples that roughly looks like this

[
    [binary_string_1 or its decimal, associated_energy_1], 
    ...
    [binary_string_n or its decimal, associated_energy_n], 
]

where n is a number of shots on a device.
Any suggestions?

Also, I would like to thank everyone for quick and detailed answers to my previous questions :heart:.

Hi @Einar_Gabbassov :slight_smile:

You can use any observable with the qml.sample function.

For example, you could create a tensor product of Pauli’s:

obs = qml.PauliZ(wires=0) @ qml.PauliX(wires=1) @ qml.PauliY(wires=2)

And then sample that:

qml.sample(obs)

Hi Christina,
Thank you for your reply. I’m aware of this possibility and sorry, but I’m not sure if your reply answers my question.

Perhaps, I should elaborate more.
Here is the Pennylane example of preparing a quantum solution
In the example, we sample the solution using the custom Hermitian which takes in 2^num_qubits by 2^num_qubits matrix. Each eigenvalue of the observable represents a solution vector encoded as decimal.
I’m looking for a scalable way to sample solutions without constructing huge matrices or doing any operations that scale exponentially with the number of qubits.

Creating a tensor observable would only scale linearly in the number of wires. It wouldn’t store the full matrix, just the necessary operation for each wire.

You could even do so iteratively to create the observable for an arbitrary number of wires:

obs = qml.PauliZ(0)
for wire in range(1, n_wires):
    obs = obs @ qml.PauliZ(wire)

You also bypass the return value entirely and query the samples from the device directly.

If the circuit ended in your computational basis state, you can look up the stored dev._samples or query dev.generate_samples

Does any of that help?

For the release coming out on Monday night (you can get now via master branch), PennyLane will have the Projector observable. You can also sample from this.

Hi Christina,
Unless I miss something obvious, the approach you propose can’t recover a solution vector. Let’s suppose that we have a two qubit system, then possible solutions to my problem are 00, 01, 10, 11. I don’t know apriori what the correct solution is (otherwise the projector observable could be pretty handy). I proceed as you suggested,

obs = qml.PauliZ(0) @ qml.PauliZ(1)
return qml.sample(obs)

This will yield a set consisting of two possible values 1 or -1.
For example, with 3 measurements I could get {1,1,-1}.
Both 1 and -1 are eigenvalues of the Z tensor Z observable. However, I have 4 possible solutions, but I can only observe two distinct eigenvalues.
We don’t have one to one relation with sampled eigenvalues and possible solution set 00, 01, 10, 11. Therefore, there is no way to recover the solution using the proposed method.

The Pennylane code in preparation of a quantum solution circumvents this issue by explicitly building a matrix with eigenvalues that range from 0 to 2^n. Each eigenvalue maps to a unique quantum solution. This means when I draw a value 10, then I know that I observed the solution 1010, and if I draw value 1 then this is the same as observing the solution 0001.

The problem with this approach that it is not practical. The matrix explodes pretty quickly into unrealistically big dimensions even with few qubits.

Given that you twice pointed me to the same solution (tensoring Pauli matrices) I feel like I’m missing something obvious here :sweat_smile:

I was missing something. Sorry. Now I understand what you want.

I see two solutions to get what you want (If I properly understand what you are asking for this time):

Option 1

Return a tuple of individual samples:

dev = qml.device('default.qubit', wires=2, shots=10)

@qml.qnode(dev)
def circ():
    qml.Hadamard(wires=0)
    qml.CNOT(wires=(0,1))
    return qml.sample(qml.PauliZ(0)), qml.sample(qml.PauliZ(1))

print(circ())
array([[ 1, -1,  1,  1, -1,  1, -1, -1,  1,  1],
       [ 1, -1,  1,  1, -1,  1, -1, -1,  1,  1]])

Option 2

Prepare the device in the state you want, return anything in your preferred computational basis, and then query device samples. For example, after running the above, you can examine:

>>> dev._samples
array([[0, 0],
       [1, 1],
       [0, 0],
       [0, 0],
       [1, 1],
       [0, 0],
       [1, 1],
       [1, 1],
       [0, 0],
       [0, 0]])

If you are on a QubitDevice, like "default.qubit", you can generate a new set of samples from the device via dev.generate_samples():

>>>  new_samples = dev.generate_samples()
>>>  new_samples
array([[0, 0],
       [1, 1],
       [1, 1],
       [0, 0],
       [1, 1],
       [1, 1],
       [0, 0],
       [1, 1],
       [1, 1],
       [1, 1]])

You need to execute the circuit before querying the samples.

Using a simple qnode return statement

If you query the samples on the device, you would just need to return some unimportant observable on at least one wire, but could then get the sample information on all the wires.

dev2 = qml.device('default.qubit', wires=5, shots = 2)

@qml.qnode(dev2)
def circ2():
    qml.Hadamard(wires=0)
    qml.CNOT(wires=(0,1))
    qml.CNOT(wires=(1,2))
    return qml.sample(qml.PauliZ(0))

circ2()

print(dev2.generate_samples()
array([[1, 1, 1, 0, 0],
       [1, 1, 1, 0, 0]])

There may be a better way of doing this, and we may add a better way of doing this in the future.

1 Like

Great!
Option 2 is exactly what I wanted!
Although, I have to say, it is completely none-obvious way of doing this.
I expected it would be a return statement from the circuit and not the device.

For example, something like this

def circuit():
    ....
    return qml.samples()

Thank you!

Hi @Einar_Gabbassov! As the current approach for returning samples always requires an observable which is to be sampled, it is slightly awkward to simply return computational basis state samples (as you note).

Adding qml.samples(), which simply returns raw device samples, is a great idea, and something we will look into adding to PennyLane.