cirq.sample_state_vector fails when the number of qubits > 32

[rht]

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.

[verult]

From Cirq sync:

• Need to investigate feasibility of alternative solutions that do not have the limit of 32

[NoureldinYosri]

The 32 limit seems interinsic to numpy, but it seems there are multiple scientific libraries running into this problem so numpy might pump up that limit numpy/numpy#5744 (comment) , however it probably won't be easy for them numpy/numpy#5744 (comment) and will take time.

The code in this issue can be rewritten in a way that doesn't require reshaping the state array, however this is not an isolated issue and we will get more occurrences of it in other places as people build larger circuits/run simulations on large GPUs.

For now I think it will be enough to write helper functions that allow us to do the needed operations but on a flattened array instead of a reshaped array. For a start a function for indexing and trasposing will be needed to fix this issue.

---

An alternative would be to increase the limit yourself, build a local version of numpy (and hope it works fine) as mentioned in https://stackoverflow.com/a/35192754

---

UPDATE: cupy also has the same limit https://github.com/cupy/cupy/blob/cb71365934dd5fa8197578ba55b74370cf07951a/cupy/_core/core.pyx#L197

[rht]

Update: I found another maximum supported dimension issue in

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

Line 79 in 16ed00d

state_vector = initial_state.reshape(qid_shape)

.

Traceback (most recent call last):
  File "/home/rht/oom.py", line 11, in <module>
    sim.simulate(qc_cirq)
  File "/home/rht/.local/lib/python3.9/site-packages/cirq/sim/simulator.py", line 495, in simulate
    return self.simulate_sweep(
  File "/home/rht/.local/lib/python3.9/site-packages/cirq/sim/simulator.py", line 510, in simulate_sweep
    return list(self.simulate_sweep_iter(program, params, qubit_order, initial_state))
  File "/home/rht/.local/lib/python3.9/site-packages/qsimcirq/qsim_simulator.py", line 529, in simulate_sweep_iter
    final_state = cirq.StateVectorSimulationState(
  File "/home/rht/.local/lib/python3.9/site-packages/cirq/sim/state_vector_simulation_state.py", line 350, in __init__
    state = _BufferedStateVector.create(
  File "/home/rht/.local/lib/python3.9/site-packages/cirq/sim/state_vector_simulation_state.py", line 79, in create
    state_vector = initial_state.reshape(qid_shape)
ValueError: maximum supported dimension for an ndarray is 32, found 34

The code I used to reproduce this is the one in quantumlib/qsim#612 (comment), but with num_qubits of 34. This time, I'm using Cirq 1.1.0. I suppose this initial_state.reshape(qid_shape) is not present in Cirq 0.14, and hence why I didn't encounter this error when posting this issue (#6031).

[rht]

I suppose cirq.sample_state_vector is not commonly used, but the new error I wrote in #6031 (comment) has a code pattern that is very commonly used. As such, this affects cuQuantum Appliance 23.03 (which uses Cirq 1.0.0), when using the cuStateVec backend (cuTensorNet should be unaffected); cc: @leofang @ahehn-nv @tlubowe.

[NoureldinYosri]

from a first glance it seems that the reshape is not actually needed ... the qid_shape is passed anyway to other functions and isn't actually used except during the calculation of probabilities which I fix in #6090 so after that PR is merged I think we can remove the reshape line and it will still work -maybe with a few tweeks-.

[leofang]

It does not make any sense to be limited by NumPy's arbitrary choice (32 or whatever). For single-node workloads the only true limitation is the 64-bit unsigned integer limit, which translates to 64 qubits, and for multi-node workloads, unlimited in practice since the state vector can be processed/stored in a distributed manner (cuQuantum has that). Also, other simulation methods like tensor network allow exceeding 64 qubits very easily, so I would advice for all QC frameworks to not make any assumptions in the max number of qubits.

[github-actions]

This issue is stale because it has been open 30 days with no activity. Remove stale label or comment or this will be closed in 30 days

[hthayko]

Hi there - any updates on the cirq 1.1.0 issue?
Basically the new cirq doesn't allow >32 qubit simulations. I would say this is a pretty serious issue.
This is preventing us from upgrading our cirq 0.14 to cirq 1.0

[NoureldinYosri]

@hthayko we just merged #6090 which fixes the original issue. The new issue from #6031 (comment) turned out to be trickier than original anticipated since all the operations of _BufferedStateVector assume that the state vector is in the correct shape.

---

Since v0.14 The StateVectorSimulationState class acquired new functionalities that were all optimized for performance. all of these optimizations assumed the state vector is in a specific shape.

If #6090 is enough fixes your issue then you can upgrade cirq 1.0. otherwise you will need to wait until the numpy issue is resolved. luckly this seems to be close see numpy/numpy#5744 (comment) and my suggested solution numpy/numpy#5744 (comment)

If you don't want to remain blocked on numpy you can try the hack in numpy/numpy#5744 (comment) until the issue is resolved.

[NoureldinYosri]

The https://github.com/quantumlib/qsim simulator has no limit on the number of qubits of the simulation. As long as you have enough RAM the qsim simulator will run the simulation with no problem.

With this I will close this issue since we have an alternative simulator that handles more than 32 qubits and the limitation imposed on the state vector simulator comes from numpy and is not intrinsic to the state vector simulator.

[rht]

You should check the code in #6031 (comment). I used qsimcirq there.

[NoureldinYosri]

@rht As said before this is a numpy issue and isn't internsic to the simulators ... if we do it in a different way then we are sacrificing efficiency and by that I mean multiplying the run time by several orders of magnitude. not to mention filling our code base with throwaway code.

When numpy starts supporting more than 32 dimensions then the state vector simulation will work out of the box. If you really need to do this kind of simulation then I suggest keeping an eye on the numpy issue numpy/numpy#5744 where they indicated that they don't have an issue with this but are lacking people to work on the fixing blockers for this, mainly the old iterators numpy/numpy#5744 (comment)

[NoureldinYosri]

@rht at your risk and as a work around you can compile your own numpy after updating the limit on number of dimensions numpy/numpy#5744 (comment)

[tanujkhattar]

@NoureldinYosri so what's your recommendation for users? That they use the c++ based interface of qsim? Can we make qsimcirq not fail for > 32 qubits ?

[NoureldinYosri]

actually after taking another look. it looks like for qsim the breakage happens during returning the result rather than while computing it. This is why I didn't notice the problem in #6031 (comment) since I didn't wait for the simulation to finish.

In that case we can create a simplified version of the StateVectorSimulationState for reporting the results.

so yea I guess we aren't blocked on numpy for this after all.

---

So next I will implement a simplified version of StateVectorSimulationState for reporting the simulation results that works within the dimensions limit.

[NoureldinYosri]

@rht does quantumlib/qsim#623 work for you?

[rht]

Yes, it looks promising, that I can finally use sim.simulate(circuit) for >32 qubits. A 1D array result works for me.

[NoureldinYosri]

@rht qsim simulator now has a new method simulate_into_1d_array which you can use to run your simulation and get the state vector as a 1D array.

[rht]

Yes, thank you so much for the simulate_into_1d_array!