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 .