ValueError: sum(pvals[:-1]) > 1.0 in pulser_emulator

[Yoric]

Just encountered the following stack when running qutipemulator on a simulation with 20 qubits.

File ~/Documents/Code/quantum-evolution-kernel/qek/data/extractors.py:429, in QutipExtractor.run(self, max_qubits)
    427 logger.debug("Executing compiled graph # %s", id)
    428 simul = QutipEmulator.from_sequence(sequence=compiled.sequence)
--> 429 counter = cast(dict[str, Any], simul.run().sample_final_state())
    430 logger.debug("Execution of compiled graph # %s complete", id)
    431 raw_data.append(compiled.graph)

File ~/.local/share/hatch/env/virtual/quantum-evolution-kernel/CUtOsxF1/quantum-evolution-kernel/lib/python3.12/site-packages/pulser_simulation/simresults.py:157, in SimulationResults.sample_final_state(self, N_samples)
    147 def sample_final_state(self, N_samples: int = 1000) -> Counter:
    148     """Returns the result of multiple measurements of the final state.
    149 
    150     Args:
   (...)
    155         measured quantum states at the end of the simulation.
    156     """
--> 157     return self.sample_state(self._sim_times[-1], N_samples)

File ~/.local/share/hatch/env/virtual/quantum-evolution-kernel/CUtOsxF1/quantum-evolution-kernel/lib/python3.12/site-packages/pulser_simulation/simresults.py:532, in CoherentResults.sample_state(self, t, n_samples, t_tol)
    517 def sample_state(
    518     self, t: float, n_samples: int = 1000, t_tol: float = 1.0e-3
    519 ) -> Counter:
    520     """Returns the result of multiple measurements at time t.
    521 
    522     Args:
   (...)
    530         quantum states at time t.
    531     """
--> 532     sampled_state = super().sample_state(t, n_samples, t_tol)
    533     if self._meas_errors is None or (
    534         self._meas_errors["epsilon"] == 0.0
    535         and self._meas_errors["epsilon_prime"] == 0
    536     ):
    537         return sampled_state

File ~/.local/share/hatch/env/virtual/quantum-evolution-kernel/CUtOsxF1/quantum-evolution-kernel/lib/python3.12/site-packages/pulser_simulation/simresults.py:145, in SimulationResults.sample_state(self, t, n_samples, t_tol)
    132 """Returns the result of multiple measurements at time t.
    133 
    134 Args:
   (...)
    142     measured quantum states at time t.
    143 """
    144 t_index = self._get_index_from_time(t, t_tol)
--> 145 return self[t_index].get_samples(n_samples)

File ~/.local/share/hatch/env/virtual/quantum-evolution-kernel/CUtOsxF1/quantum-evolution-kernel/lib/python3.12/site-packages/pulser/result.py:82, in Result.get_samples(self, n_samples)
     73 def get_samples(self, n_samples: int) -> Counter[str]:
     74     """Takes multiple samples from the sampling distribution.
     75 
     76     Args:
   (...)
     80         Samples of bitstrings corresponding to measured quantum states.
     81     """
---> 82     dist = np.random.multinomial(n_samples, self._weights())
     83     return Counter(
     84         {
     85             np.binary_repr(i, self._size): dist[i]
     86             for i in np.nonzero(dist)[0]
     87         }
     88     )

File numpy/random/mtrand.pyx:4370, in numpy.random.mtrand.RandomState.multinomial()

ValueError: sum(pvals[:-1]) > 1.0

I got this result from the following notebook, https://github.com/pasqal-io/quantum-evolution-kernel/blob/afa74d99ac5b05a1e91b0ad90491115544c6b09c/examples/tutorial%201%20-%20Using%20a%20Quantum%20Device%20to%20Extract%20Machine-Learning%20Features.ipynb (you'll need to pip install quantum-evolution-kernel first), by replacing MAX_QUBITS with 50 for the sake of an experiment.

[HGSilveri]

@pasqal-io/emulation-team This might be a good first issue :)

[Yoric]

Versions used:

• numpy 1.26.4
• pulser 1.1.1
• pulser-core 1.1.1
• pulser-pasqal 0.20.1
• pulser-simulation 1.1.1

[pablolh]

I could reproduce the issue. For me it's a floating point error accumulation in the normalization of the weights done in qutip_result.py in pulsersim. e.g.:

>>> t = np.random.rand(30000000); sum(t / sum(t))
1.0000000000001992

The check in np.random.multinomial from result.py is sum(pvals[:-1]) < 1 and could throw when unlucky.

I think np.random.multinomial is maybe okay for small state vectors but for large distributions it's not good: the results are given as a 2^N sized array, which is very suboptimal when taking just a few samples (array is full of unused zeros). I really think this function is designed for emulating 6-faced dices, not 2^N distributions from quantum states.

I offer to write custom multinomial sampling code for get_samples in Results that directly creates the counter.

The algo:

• create empty bitstring counter
• generate a set of n_samples floats between 0 and 1 uniformly distributed
• iterate over the weights with index, and accumulate them in a variable
• when the accumulator goes above a uniform sampled float, add that index to the counter and remove that sampled float

Any opinion @HGSilveri ?

[HGSilveri]

No objections @pablolh , all I ask is a small benchmark once you're done, to see how it stacks up :)