Choosing a Simulation Method

Choosing a Simulation Method

Choose between statevector, matrix product state (MPS), and general tensor-network simulation based on the circuit structure, requested outputs, accuracy requirements, and available resources.

Version: FlagQuantum 0.3.1. All three complete scripts below were verified on macOS Apple Silicon with Python 3.12.14 and PyTorch 2.14.1 using local CPU execution. They require the core installation from Installation and Compatibility; the demonstrated native PyTorch paths do not require the optional JAX or cotengra packages.

Formula/code review — October 20, 2026: The equations are paired with code excerpts, and the local CPU examples were re-run with Python 3.12.14 and PyTorch 2.13.0. See the review record for executed commands, distributed checks, and untested hardware/provider paths. Earlier environment and output records are retained.

1. A Practical Starting Point

Method	Representation	Good starting scenarios	Main limitation
statevector	One complex amplitude for every computational-basis state	Small reference calculations, debugging, arbitrary circuit connectivity, and full-state inspection	State storage grows exponentially with qubit count.
mps	A chain of small tensors connected by bond indices	Circuits that maintain modest entanglement across cuts in the chosen qubit ordering, especially local chain-like interactions	Bond dimensions can grow exponentially; long-range interactions and decompositions add cost.
tensor_network	A circuit-derived network of tensors evaluated by contraction	Structured circuits with manageable contraction width, particularly when only a few outputs are needed	Contraction order and intermediate tensor sizes can make execution expensive or infeasible.

MPS is itself a tensor-network representation. The distinction here is between FlagQuantum's chain-state mps mode and its more general circuit-contraction tensor_network mode.

Start with statevector when the circuit fits comfortably in memory. Use it as a correctness reference on a reduced problem before moving to a compressed or contraction-based method. There is no universal qubit-count threshold or fastest method: depth, connectivity, entanglement, outputs, batching, and gradient requirements all matter.

2. Statevector: A Direct Reference Representation

An n-qubit pure state has 2**n complex amplitudes. Statevector simulation stores that dense vector and updates it as gates execute. It represents general pure states without imposing a low-entanglement structure.

When to Use It

• The full state fits in memory, including working space and any training graph.
• You want a small, straightforward reference for probabilities, expectations, or gradients.
• You need to inspect many amplitudes or the full probability distribution.

Limits

The amplitude array alone requires approximately 2**n * bytes_per_complex_value bytes per batch entry. For 30 qubits, that is 8 GiB with complex64 or 16 GiB with complex128, before buffers, output arrays, optimizer state, or backward storage. A single-device allocation estimate is not a guarantee that execution fits.

Requesting only one expectation can reduce output size, but it does not remove the dense statevector maintained by this method. These noiseless pure-state examples also do not establish mixed-state or noise support; density-matrix and trajectory workflows are separate capabilities.

3. MPS: Compress the State Along a Chain

An MPS factors amplitudes into a chain of site tensors. Neighboring tensors share indices whose sizes are called bond dimensions. A bond must be large enough to represent the state's correlations across the corresponding cut in the chain.

For fixed physical dimension and a maximum bond dimension chi, state storage scales roughly as O(n * chi**2) per batch entry. This is useful only while chi stays manageable. General states can require exponentially large bonds, and gate application or differentiation needs additional working memory.

When to Use It

• The evolving state has modest entanglement across the relevant chain cuts.
• Interactions are mostly local in a useful qubit ordering.
• You need observables or samples without constructing every amplitude.

Limits and Accuracy

MPS is not inherently approximate. With no bond cap and no truncation, it can represent an exact state up to numerical precision, while allowing bonds to grow. A low bond cap or nonzero truncation threshold can reduce cost by discarding information, introducing approximation error.

Low entanglement in a final state does not guarantee a cheap run: intermediate states may need larger bonds. Non-neighbor gates can require extra work, and changing the qubit ordering can change which cuts are difficult. Converting an arbitrary dense input state to MPS also starts from an exponentially large input description.

Use explicit accuracy settings and inspect truncation information when approximating. A bond cap alone is not a certified bound on an observable's error. Changes in rank and tensor factorizations can also affect numerical gradient behavior; compare gradients on a tractable reference problem.

4. General Tensor Networks: Contract the Requested Quantity

A circuit can be represented as gate and state tensors joined by shared indices. Computing an amplitude or expectation involves summing over those internal indices, called tensor contraction.

The contraction path determines which tensors are combined first. Different paths can have very different intermediate sizes, memory requirements, and arithmetic cost, even for the same circuit.

When to Use It

• The circuit's structure permits small contraction intermediates.
• You need selected expectations, amplitudes, or marginals rather than the full state.
• You can inspect or control contraction planning for your workload.

Limits

Tensor networks do not make arbitrary circuits efficient. Highly connected or deep circuits can have expensive contractions. Path search itself can cost time. Slicing a contraction into pieces can reduce peak intermediate storage while increasing work; slicing alone does not imply multi-GPU or multi-node execution.

General tensor-network contraction can be exact up to floating-point precision. Approximation, when used by a particular workflow, is a separate choice. In FlagQuantum 0.3.1, the general tensor-network execution/training capability is labeled Experimental; noise-channel instructions in this mode are unsupported. Consult the release capability catalog for the exact scope rather than extrapolating from a successful small example.

5. Request Outputs That Match Your Goal

Goal	Suggested output strategy
Check a local observable	Request fq.expectation(...).
Inspect a few qubits	Request a marginal such as fq.probabilities(qubits=(0,)).
Compare methods on a tiny circuit	Compare outputs against an analytical result and a statevector reference.
Inspect every amplitude	Use statevector if practical; a dense conversion from another representation can erase its memory advantage.
Train a model	Validate both predictions and gradients for the actual execution path.

A full statevector or full probability array has 2**n entries regardless of how efficiently the internal state was represented. Computing a large dense output can therefore dominate a compressed calculation. Marginals also grow exponentially with the number of selected qubits.

The scripts below request small outputs and explicitly select each mode. They disable backend fallback and disallow approximate execution. These settings do not make an unsupported workload supported, impose a memory budget, or guarantee a particular numerical tolerance.

Compare Representations and Quantities in Code

Formula / quantity	Corresponding code excerpt	What the scripts verify
$N_{\mathrm{amplitudes}}=B2^n$	statevector_elements = batch_size * (2 ** n_qubits)	This is an allocation estimate, not a measured peak-memory value.
$M_{\mathrm{state}}=B2^n s$	statevector_elements * torch.tensor([], dtype=torch.complex64).element_size()	s is bytes per complex value; exclude work buffers and autograd storage.
$\lvert\psi\rangle=\cos(0.3)\lvert00\rangle+\sin(0.3)\lvert11\rangle$	fq.Circuit(2).ry(0, 0.6).cx(0, 1)	The same circuit is evaluated in all three modes.
$\langle Z_0\rangle=\cos(0.6),\quad\langle Z_0Z_1\rangle=1$	expected_z = torch.cos(torch.tensor(0.6)).item()	compare_simulation_methods.py checks the named expectations.
$d\langle Z_0\rangle/d\theta=-\sin\theta$	expected = -torch.sin(torch.tensor(0.6))	compare_simulation_gradients.py checks all three gradient paths, with atol=1e-5.
$\chi_{\mathrm{product}}=1,\quad\chi_{\mathrm{GHZ}}=2$	state.bond_dims, state.max_bond	inspect_mps_bonds.py checks normalization and the GHZ maximum bond.

MPS storage roughly follows $O(n\chi^2)$ for qubits with bounded bond dimension; this describes the representation, not the peak memory or runtime of a particular implementation. The memory expressions are estimates; the scripts below validate values, gradients, and bonds on small circuits.

Example 1: Compare the Same Circuit Across All Three Methods

Save as compare_simulation_methods.py:

import torch
import flagquantum as fq

def build_circuit():
    return fq.Circuit(2).ry(0, 0.6).cx(0, 1)

expected_z = torch.cos(torch.tensor(0.6)).item()
for mode in ("statevector", "mps", "tensor_network"):
    result = fq.run(
        build_circuit(),
        outputs=[
            fq.expectation(fq.Z(0), name="z0"),
            fq.expectation(fq.Z(0) @ fq.Z(1), name="zz"),
            fq.probabilities(qubits=(0,)),
        ],
        options=fq.ExecutionOptions(
            mode=mode,
            backend="pytorch",
            device="cpu",
            allow_backend_fallback=False,
            allow_approximate=False,
        ),
    )
    z0 = result.expectation("z0").item()
    zz = result.expectation("zz").item()
    probabilities = result.probabilities[0]
    print(f"{mode}: <Z0>={z0:.6f}, <Z0 Z1>={zz:.6f}")
    print("  Qubit 0 probabilities:", [round(p.item(), 6) for p in probabilities])
    assert abs(z0 - expected_z) < 1e-6
    assert abs(zz - 1.0) < 1e-6
    assert abs(probabilities.sum().item() - 1.0) < 1e-6

print("Three-method output checks passed.")

Run python compare_simulation_methods.py. Expected output:

statevector: <Z0>=0.825336, <Z0 Z1>=1.000000
  Qubit 0 probabilities: [0.912668, 0.087332]
mps: <Z0>=0.825336, <Z0 Z1>=1.000000
  Qubit 0 probabilities: [0.912668, 0.087332]
tensor_network: <Z0>=0.825336, <Z0 Z1>=1.000000
  Qubit 0 probabilities: [0.912668, 0.087332]
Three-method output checks passed.

The circuit prepares cos(0.3)|00> + sin(0.3)|11>. Consequently, <Z0> = cos(0.6) and <Z0 Z1> = 1. The displayed probabilities concern only qubit 0; no full-register probability vector is requested.

This is a correctness comparison, not a performance benchmark. For a two-qubit circuit, execution overhead can dominate, so its runtime would not establish which method is best for a larger workload.

6. Inspect MPS Bonds on Simple States

A product state needs bond dimension 1. A GHZ state is entangled, but needs only bond dimension 2 across each nontrivial chain cut. Thus, the presence of entanglement alone does not imply that an MPS will be large; the required ranks across cuts matter.

The following script uses the documented native MPS entry point to inspect the representation. Its max_bond and cutoff arguments belong to run_mps; they are not fields of the generic fq.ExecutionOptions object in 0.3.1.

Example 2: Compare Product-state and GHZ Bonds

Save as inspect_mps_bonds.py:

import torch
import flagquantum as fq
from flagquantum.simulation.mps import run_mps

product = fq.Circuit(6)
ghz = fq.Circuit(6).h(0)
for qubit in range(1, 6):
    ghz.cx(qubit - 1, qubit)

for label, circuit in (("Product", product), ("GHZ", ghz)):
    state = run_mps(circuit, max_bond=None, cutoff=0.0)
    print(f"{label} bond dimensions:", state.bond_dims)
    print(f"{label} maximum bond:", state.max_bond)
    assert torch.allclose(state.state_norm(), torch.ones(1), atol=1e-6)

assert run_mps(ghz, max_bond=None, cutoff=0.0).max_bond == 2
print("MPS bond checks passed.")

Run python inspect_mps_bonds.py. Expected output:

Product bond dimensions: (1, 1, 1, 1, 1)
Product maximum bond: 1
GHZ bond dimensions: (2, 2, 2, 2, 2)
GHZ maximum bond: 2
MPS bond checks passed.

Here max_bond=None and cutoff=0.0 request untruncated execution. The native result is an MPSState, not the ExecutionResult returned by fq.run. The norm checks verify normalization for these examples; they do not establish an error guarantee for an approximate run.

7. Check Gradients Before Switching a Training Workload

Matching forward values is necessary but not sufficient for training. The gradient path, numerical stability, and supported observable must also match your requirements. Test the derivative on a small problem before training a larger model.

Example 3: Compare a Parameter Gradient

Save as compare_simulation_gradients.py:

import torch
import flagquantum as fq

def build_circuit(parameters):
    return fq.Circuit(2).ry(0, parameters[0]).cx(0, 1)

for mode in ("statevector", "mps", "tensor_network"):
    model = fq.Module(
        build_circuit,
        n_parameters=1,
        init=torch.tensor([0.6]),
        device="cpu",
        policy=fq.RuntimePolicy(
            execution_options=fq.ExecutionOptions(
                mode=mode,
                backend="pytorch",
                device="cpu",
                allow_backend_fallback=False,
                allow_approximate=False,
            ),
            observable="z",
            observable_qubits=(0,),
        ),
    )
    loss = model().mean()
    loss.backward()
    gradient = next(model.parameters()).grad[0]
    expected = -torch.sin(torch.tensor(0.6))
    assert torch.isfinite(gradient)
    assert torch.allclose(gradient, expected, atol=1e-5)
    print(f"{mode}: <Z0>={loss.item():.6f}, gradient={gradient.item():.6f}")

print("Three-method gradient checks passed.")

Run python compare_simulation_gradients.py. Expected output:

statevector: <Z0>=0.825336, gradient=-0.564642
mps: <Z0>=0.825336, gradient=-0.564642
tensor_network: <Z0>=0.825336, gradient=-0.564642
Three-method gradient checks passed.

Each iteration creates a fresh module with the same initial angle. The analytical derivative is -sin(0.6). The script checks that each gradient is finite and agrees with this reference within 1e-5.

Agreement for this circuit does not guarantee gradients for every gate, rank configuration, output type, noise model, or backend. In particular, the release catalog describes the general tensor-network execution/training surface as experimental even though selected exact gradient workflows have their own scoped support statements.

8. Make a Method-selection Decision

1. Specify the result. Decide whether you need a few observables, samples, marginals, or the entire state, and whether backward computation is required.
2. Estimate the dense baseline. Include batch size, precision, working space, and gradients rather than just the amplitude array.
3. Inspect structure. Consider MPS bond growth and interaction locality, or tensor-network contraction intermediates and output size.
4. Validate a reduced case. Compare values and, when needed, gradients with analytical expectations or a statevector reference.
5. Measure the real workload. Record peak memory and forward/backward time on the intended hardware. Increase circuit size gradually and keep accuracy settings explicit.

If statevector fits and meets your needs, it is a useful initial choice. If it does not fit, examine MPS and tensor-network structure before assuming either method will solve the capacity problem. Qubit count alone is insufficient to predict performance.

9. Release Scope and Common Pitfalls

The 0.3.1 capability catalog distinguishes Production supported local statevector simulation/training, Development evidence for the combined differentiable/sharded MPS training entry, and Experimental general tensor-network execution/training. Other specifically qualified gradient or multi-node workflows have separate entries; their maturity does not automatically extend to every use of a representation.

Pitfall	Better approach
Treating MPS as always approximate	Check whether truncation is enabled and whether bonds can grow.
Assuming a small bond cap has negligible error	Inspect truncation evidence and compare observables against a reference.
Assuming tensor networks are always faster	Evaluate planning cost, contraction intermediates, and the requested outputs.
Asking a compressed method for an enormous dense state	Request the quantities you actually need.
Using device="mps" to select matrix product states	Use mode="mps". Apple's MPS device naming is a separate concept and does not establish device support here.
Treating local tensor slicing as distributed execution	Verify the separately supported distribution semantics and transport path.
Assuming a CPU demonstration certifies GPU behavior	Validate the intended device and workload separately.
Silently allowing a different backend during a comparison	Disable fallback and report any explicit representation conversions.

Continue Learning and References

• Distributed Simulation and Training: sharding, communication preflight, and a verified two-worker training workflow.
• Accelerator Setup: runtime detection, explicit device execution, and gradient/optimizer checks.
• Building and Simulating Circuits: gates and result requests.
• Training Parameterized Circuits: objectives and optimizer updates.
• Hybrid Classical-Quantum Models: differentiable classical-quantum composition.
• Qubit Ordering and Result Conventions: probability-vector ordering and marginal selections.
• Wiki Directory: published guides.
• Release MPS entry points, MPS state, tensor-network entry points, and execution options.