[BUG] BasisState on default.clifford triggers dense state allocation (fails at 64 qubits)

[JeremieGince]

Expected behavior

default.clifford should allow large-scale Clifford simulations (e.g., 64 qubits) when the circuit contains only Clifford operations and Pauli observables.

In particular, using a stabilizer-compatible state preparation such as qml.BasisState should not trigger construction of a full (2^n) state vector. The device should internally convert the basis state into a tableau representation and successfully return all requested Pauli expectation values.

Actual behavior

When qml.BasisState(np.zeros(64), wires=range(64)) is queued inside a QNode executed on default.clifford, PennyLane attempts to call state_vector() for the preparation operator. This leads to allocation of a dense tensor of size (2^{64}), which overflows memory and raises:

ValueError: array is too big; arr.size * arr.dtype.itemsize is larger than the maximum possible size.

As a result, Clifford expectation values cannot be computed for large qubit counts, even though the circuit and observables are fully Clifford-compatible and should be handled purely in stabilizer/tableau form.

Additional information

This function is important for the development of the MatchCake package: https://github.com/MatchCake/MatchCake.

Source code

import numpy as np
import pennylane as qml
from pennylane import qnode
from pennylane.operation import StatePrepBase
from pennylane.pauli import string_to_pauli_word


@qml.qnode(qml.device("default.clifford", tableau=True))
def GHZStatePrep(num_wires):
    # From: https://pennylane.ai/qml/demos/tutorial_clifford_circuit_simulations
    qml.Hadamard(wires=[0])
    for wire in range(num_wires):
        qml.CNOT(wires=[wire, wire + 1])
    return qml.expval(qml.Z(0) @ qml.Z(num_wires - 1))


expval = GHZStatePrep(num_wires=64)
print(f"GHZStatePrep: {expval}")


def majorana_to_pauli(i: int) -> qml.pauli.PauliWord:
    # From: https://github.com/MatchCake/MatchCake/blob/f277d226252e7733374743378e78620c68600b73/src/matchcake/utils/majorana.py#L123
    k = int(i // 2)
    gate = "X" if i % 2 == 0 else "Y"
    majorana_str = "".join(["Z"] * k + [gate])
    return string_to_pauli_word(majorana_str)


def compute_clifford_expvals(state_prep_op: StatePrepBase):
    # From: https://github.com/MatchCake/MatchCake/blob/f277d226252e7733374743378e78620c68600b73/src/matchcake/devices/expval_strategies/clifford_expval/clifford_expval_strategy.py#L24
    wires = state_prep_op.wires
    triu_indices = np.triu_indices(2 * len(wires), k=1)

    @qnode(qml.device("default.clifford", tableau=True))
    def clifford_circuit():
        state_prep_op.queue()
        return [qml.expval(majorana_to_pauli(mu) @ majorana_to_pauli(nu)) for mu, nu in zip(*triu_indices)]

    return clifford_circuit()


state_prep_op = qml.BasisState(np.zeros(64), wires=range(64))
expval = compute_clifford_expvals(state_prep_op)
print(f"Clifford Expvals: {expval}")

Tracebacks

E:\Github\MatchCake\.venv\Scripts\python.exe C:\Users\gince\AppData\Roaming\JetBrains\PyCharm2025.3\scratches\pennylane_clifford\issue.py 
GHZStatePrep: 1.0
Traceback (most recent call last):
  File "C:\Users\gince\AppData\Roaming\JetBrains\PyCharm2025.3\scratches\pennylane_clifford\issue.py", line 42, in <module>
    state_prep_op = qml.BasisState(np.zeros(64), wires=range(64))
  File "C:\Users\gince\AppData\Roaming\JetBrains\PyCharm2025.3\scratches\pennylane_clifford\issue.py", line 38, in compute_clifford_expvals
    
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\workflow\qnode.py", line 882, in __call__
    return self._impl_call(*args, **kwargs)
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\workflow\qnode.py", line 855, in _impl_call
    res = qml.execute(
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\workflow\execution.py", line 244, in execute
    results = run(tapes, device, config, inner_transform)
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\workflow\run.py", line 286, in run
    results = inner_execute(tapes)
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\workflow\run.py", line 251, in inner_execute
    results = device.execute(transformed_tapes, execution_config=execution_config)
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\devices\modifiers\simulator_tracking.py", line 28, in execute
    results = untracked_execute(self, circuits, execution_config)
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\devices\modifiers\single_tape_support.py", line 30, in execute
    results = batch_execute(self, circuits, execution_config)
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\devices\default_clifford.py", line 504, in execute
    return tuple(
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\devices\default_clifford.py", line 505, in <genexpr>
    self.simulate(c, seed=s, debugger=self._debugger) for c, s in zip(circuits, seeds)
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\devices\default_clifford.py", line 561, in simulate
    qml.math.reshape(prep.state_vector(wire_order=list(circuit.op_wires)), (1, -1))[0],
  File "E:\Github\MatchCake\.venv\lib\site-packages\pennylane\ops\qubit\state_preparation.py", line 182, in state_vector
    ket = math.zeros((2,) * num_wires)
  File "E:\Github\MatchCake\.venv\lib\site-packages\autoray\autoray.py", line 81, in do
    return func(*args, **kwargs)
ValueError: array is too big; `arr.size * arr.dtype.itemsize` is larger than the maximum possible size.

Process finished with exit code 1

System information

Backend tkagg is interactive backend. Turning interactive mode on.
Name: PennyLane
Version: 0.41.1
Summary: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Home-page: https://github.com/PennyLaneAI/pennylane
Author: 
Author-email: 
License: Apache License 2.0
Location: e:\github\matchcake\.venv\lib\site-packages
Requires: appdirs, autograd, autoray, cachetools, diastatic-malt, networkx, numpy, packaging, pennylane-lightning, requests, rustworkx, scipy, tomlkit, typing-extensions
Required-by: matchcake, PennyLane_Lightning
Platform info:           Windows-10-10.0.26100-SP0
Python version:          3.10.11
Numpy version:           2.2.5
Scipy version:           1.15.2
Installed devices:
- default.clifford (PennyLane-0.41.1)
- default.gaussian (PennyLane-0.41.1)
- default.mixed (PennyLane-0.41.1)
- default.qubit (PennyLane-0.41.1)
- default.qutrit (PennyLane-0.41.1)
- default.qutrit.mixed (PennyLane-0.41.1)
- default.tensor (PennyLane-0.41.1)
- null.qubit (PennyLane-0.41.1)
- reference.qubit (PennyLane-0.41.1)
- lightning.qubit (PennyLane_Lightning-0.41.0)

Existing GitHub issues

• I have searched existing GitHub issues to make sure the issue does not already exist.

[josh146]

Hi @JeremieGince, thanks for catching this! Indeed, StatePrep operations should be defaulting to an efficient Clifford operation where possible, this is something we should fix.

In the meantime, just wanted to ask if you are blocked by this bug, or if forcing the stateprep to be decomposed (e.g., BasisState -> List[PauliX]) is a reasonable workaround for MatchCake?

[JeremieGince]

Hi @josh146, thank you very much! I’m not sure why this hadn’t occurred to me earlier — it’s a simple and elegant solution that resolves my problem for now.

For reference, here is the code implementing the workaround:

def compute_clifford_expvals(state_prep_op: StatePrepBase):
    # From: https://github.com/MatchCake/MatchCake/blob/f277d226252e7733374743378e78620c68600b73/src/matchcake/devices/expval_strategies/clifford_expval/clifford_expval_strategy.py#L24
    wires = state_prep_op.wires
    triu_indices = np.triu_indices(2 * len(wires), k=1)

    @qnode(qml.device("default.clifford", tableau=True))
    def clifford_circuit():
        if isinstance(state_prep_op, BasisState):
            for op in state_prep_op.decomposition():
                op.queue()
        else:
            raise NotImplementedError("Only BasisState is implemented for Clifford subroutine.")
        return [qml.expval(majorana_to_pauli(mu) @ majorana_to_pauli(nu)) for mu, nu in zip(*triu_indices)]

    return clifford_circuit()


state_prep_op = qml.BasisState(np.random.choice([0, 1], size=64), wires=range(64))
expval = compute_clifford_expvals(state_prep_op)
print(f"Clifford Expvals: {expval}")

If you consider this sufficient to treat the issue as resolved, please feel free to close it.

[josh146]

Nice, glad to see that it worked @JeremieGince! Let us know if anything else comes up that we can help with.

Regarding this issue, let's keep it open, as it is a bug we would like to fix and address.

[OmkarSarkar204]

Hi! @josh146 I’d like to work on this issue. I’ve started looking into it and will open a PR once I have a working fix.

[CatalinaAlbornoz]

Hi @OmkarSarkar204 ,

Thanks for offering to work on this! It's great to see the enthusiasm. In this case however, it looks like the solution to this issue might not be obvious. Do you you already have a solution in mind?

Given that this is a tricky issue, I would recommend that you first share the idea for your solution here before spending a lot of time opening a PR with a fix. This way we can review at a high level and identify whether the solution is likely to work or not.

[CatalinaAlbornoz]

Hi @OmkarSarkar204 ,

My colleague Utkarsh had a suggestion that you may find useful!

The most robust way to handle this is to use the Basis state decomposition to build and update the Tableau instead of the state vector representation of the prep operation.

stim_tableau = stim.Tableau.from_state_vector(math.reshape(prep.state_vector(wire_order=list(circuit.op_wires)), (1, -1))[0], endian="big")

[OmkarSarkar204]

Hi @CatalinaAlbornoz , sorry for the delay, I am currently traveling.

Thanks for the suggestion from Utkarsh. That aligns with my plan to handle BasisState using its decomposition to update the tableau directly, bypassing the state_vector() call and avoiding dense state allocation.

I haven’t dived into the code yet, but I plan to reproduce the issue first this weekend once I reach home. After that, I’ll finalize my approach and share the details here before starting the implementation.

[CatalinaAlbornoz]

Sounds good! Enjoy your travel @OmkarSarkar204 !

[OmkarSarkar204]

Hi @CatalinaAlbornoz ,

I was able to reproduce the crash,

ValueError: array is too big; `arr.size * arr.dtype.itemsize` is larger than the maximum possible size.

and it appears to be caused by an attempt to build the full state vector, which runs out of memory.

I located the section handling state preparation (and the TODO note left there # TODO: Add a method to prepare directly from a Tableau ). I’ve implemented a fix that intercepts BasisState specifically. Instead of letting it trigger the dense state vector creation, I now iterate through its decomposition and append the resulting gates directly to the Stim circuit.

Does this approach look reasonable, and would it be ok for me to proceed with a PR based on this?

[CatalinaAlbornoz]

Hi @OmkarSarkar204 , thanks for proceeding with the PR! Let's continue the conversation there.