Feature request: more efficient state sampling

[kwkbtr]

We currently use QuantumState::sampling for state sampling.
The current implementation involves random number generation for sampling_count times:

qulacs/src/cppsim/state.hpp

Lines 668 to 674 in 00040b9

	for (UINT count = 0; count < sampling_count; ++count) {
	double r = random.uniform();
	auto ite =
	std::lower_bound(stacked_prob.begin(), stacked_prob.end(), r);
	auto index = std::distance(stacked_prob.begin(), ite) - 1;
	result.push_back(index);
	}

For our use cases, it is often the case that the sampling_count is much larger than the dimension of a state vector (e.g. sampling_count = 10^12 while qubit_count is around 10). In such a case, the random number generation is quite time-consuming.

We only need statistics of sampling for our use cases (i.e. 00...0 occurred m_1 times, 00...1 occurred m_2 times, ...).
For such use cases, it might be more efficient if we sample the statistics by using a multinomial distribution with parameters given by probabilities of a single measurement on the quantum state.

[m-ymzk]

Hi, @kwkbtr san,
Since uniform random numbers are used for sampling,
If the qubit_count is as small as 10, you might get the same result in the following steps.

1. Get normalized state_vector
2. Calculates the square of the norm of each element ( Or element * conj(element) )
3. multiply each probability by 10^12

Or am I misunderstanding something?

[kwkbtr]

@m-ymzk san,
That procedure will give expectation values of counts for each measurement result, which are not random.
What we expect from sampling is a single realization of counts for repeated measurements in a probabilistic simulation.
The latter is random, hence a procedure generating it necessarily uses a random number generator, while your procedure does not.

[m-ymzk]

So there is a case that a value including a stochastic variation is required instead of an ideal value.
I understand. Thank you.

[snuffkin]

How about using scipy multinomial?
I made a sample.

from scipy.stats import multinomial
import numpy as np
import time

# arrange
n = 1e12 # number of sampling
uniform = np.random.rand(2**10)
p = uniform / np.sum(uniform) # probability distribution

# act
start_time = time.time()
multinomial.rvs(n, p)
end_time = time.time()

print(f"{end_time - start_time:.6f} [s]")

When I ran this it completed within 1 millisecond, sampling 10^12 times at 10 qubits.

In practice, we may need to call scipy's C++ interface from qulacs.

[kwkbtr]

Yes, we are currently using multinomial in numpy.
I just thought that this feature would be useful for many users and having this in qulacs may be good for better performance.