cirq.sample_state_vector fails when the number of qubits > 32

Description of the issue

I use the cirq.sample_state_vector to extract counts from the statevector. But it appears to not scale beyond 32 qubits due to a NumPy limitation, which hardcodes the maximum number of dimensions to be 32: numpy/numpy#5744.

How to reproduce the issue

I'm running the following code in an NVIDIA cuQuantum Docker instance, because a CPU simulation of the code would be slow. However, this shouldn't affect that the error still happens.

import cirq
import qsimcirq

gpu_mode = 2
options = qsimcirq.QSimOptions(gpu_mode=gpu_mode)
simulator = qsimcirq.QSimSimulator(options)

n = 33
qc = cirq.Circuit()
qs = cirq.LineQubit.range(n)
for i in range(n):
    qc.append(cirq.H(qs[i]))
    qc.append(cirq.measure(qs[i]))

result = simulator.simulate(qc)
sv = result.final_state_vector
cirq.sample_state_vector(sv, range(n))

Error message:

Traceback (most recent call last):
  File "cirq_huge_circuit.py", line 21, in <module>
    cirq.sample_state_vector(sv, range(n))
  File "/opt/conda/lib/python3.8/site-packages/cirq/sim/state_vector.py", line 214, in sample_state_vector
    probs = _probs(state_vector, indices, shape)
  File "/opt/conda/lib/python3.8/site-packages/cirq/sim/state_vector.py", line 322, in _probs
    tensor = np.reshape(state, qid_shape)
  File "<__array_function__ internals>", line 180, in reshape
  File "/opt/conda/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 298, in reshape
    return _wrapfunc(a, 'reshape', newshape, order=order)
  File "/opt/conda/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 57, in _wrapfunc
    return bound(*args, **kwds)
ValueError: maximum supported dimension for an ndarray is 32, found 33

Cirq version

cirq-core 0.14

But should apply to 1.1.x as well:

Cirq/cirq-core/cirq/sim/state_vector.py

Lines 326 to 354 in 8b97aa5

    
           def _probs(state: np.ndarray, indices: Sequence[int], qid_shape: Tuple[int, ...]) -> np.ndarray: 
        
               """Returns the probabilities for a measurement on the given indices.""" 
        
               tensor = np.reshape(state, qid_shape) 
        
               # Calculate the probabilities for measuring the particular results. 
        
               if len(indices) == len(qid_shape): 
        
                   # We're measuring every qudit, so no need for fancy indexing 
        
                   probs = np.abs(tensor) ** 2 
        
                   probs = np.transpose(probs, indices) 
        
                   probs = probs.reshape(-1) 
        
               else: 
        
                   # Fancy indexing required 
        
                   meas_shape = tuple(qid_shape[i] for i in indices) 
        
                   probs = ( 
        
                       np.abs( 
        
                           [ 
        
                               tensor[ 
        
                                   linalg.slice_for_qubits_equal_to( 
        
                                       indices, big_endian_qureg_value=b, qid_shape=qid_shape 
        
                                   ) 
        
                               ] 
        
                               for b in range(np.prod(meas_shape, dtype=np.int64)) 
        
                           ] 
        
                       ) 
        
                       ** 2 
        
                   ) 
        
                   probs = np.sum(probs, axis=tuple(range(1, len(probs.shape)))) 
        
               # To deal with rounding issues, ensure that the probabilities sum to 1. 
        
               return probs / np.sum(probs)

. np.reshape is not the only operation that has the <=32 dimension constraint.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	def _probs(state: np.ndarray, indices: Sequence[int], qid_shape: Tuple[int, ...]) -> np.ndarray:
	"""Returns the probabilities for a measurement on the given indices."""
	tensor = np.reshape(state, qid_shape)
	# Calculate the probabilities for measuring the particular results.
	if len(indices) == len(qid_shape):
	# We're measuring every qudit, so no need for fancy indexing
	probs = np.abs(tensor) ** 2
	probs = np.transpose(probs, indices)
	probs = probs.reshape(-1)
	else:
	# Fancy indexing required
	meas_shape = tuple(qid_shape[i] for i in indices)
	probs = (
	np.abs(
	[
	tensor[
	linalg.slice_for_qubits_equal_to(
	indices, big_endian_qureg_value=b, qid_shape=qid_shape
	)
	]
	for b in range(np.prod(meas_shape, dtype=np.int64))
	]
	)
	** 2
	)
	probs = np.sum(probs, axis=tuple(range(1, len(probs.shape))))

	# To deal with rounding issues, ensure that the probabilities sum to 1.
	return probs / np.sum(probs)

Description

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions