rqutils.sqd module

Sample-based quantum diagonalization of general Pauli-sum Hamiltonians

Overview

SQD is an algorithm for finding an approximate ground eigenpair of a very large (computationally intractable) Hamiltonian by projecting it onto a subspace. In our case, the Hamiltonian is expressed as a linear combination of Pauli strings, and the subspace is identified from a (possibly redundant) set of bitstrings (states).

Typically, the projected Hamiltonian is itself still too large to be stored in memory as a matrix, requiring some matrix-free method to solve the eigenvalue problem. The technical challenge is then threefold:

Extracting the set \(S\) of unique bitstrings from the input list.
Computing (on the fly) the matrix elements \(\langle j | Q | k \rangle\) for all \(j, k \in S\) for each term \(Q\) of the Hamiltonian.
Solving the eigenvalue problem.

The first point will have to be done with np.unique or an equivalent accelerated function. For the last point, we can use the ground_locg solver provided in this package, which takes a matrix-vector application function and an initial-guess vector as inputs. The central task here is thus providing the matvec function that covers the second point.

Usage examples can be found at examples/sqd.py.

Algorithm

The algorithm takes advantage of the symplectic representation of Pauli strings. A term in the Hamiltonian is a product of a real coefficient \(\alpha\) and a Pauli string \(Q\). In the symplectic representation,

\[\alpha Q = \alpha (-i)^{xz} \left(Z^{z_{n-1}} \otimes \cdots Z^{z_{0}}\right) \left(X^{x_{n-1}} \otimes \cdots X^{x_{0}}\right)\]

where \(x\) (X signature) and \(z\) (Z signature) are binary vectors of length \(n\) (number of qubits) and \(xz\) represents their inner product. For a given bitstring \(s = [s_{n-1}, \dots, s_{0}]\), the X signature of the Hamiltonian term determines the existence and location of the matrix element, and the Z signature gives the sign.

In the preparation stage of the algorithm, the terms of the input Hamiltonian are grouped by the X signature (for example, XIZ and YZI will belong to the same group). For each X signature, there will be multiple Z signatures and the corresponding phased coefficients (\(\alpha (-i)^{xz}\) above). The X and Z signatures are bit-packed into arrays of 8-bit integers.

The input states are then sorted and similarly bit-packed to allow bitwise operations between the states and the X/Z signatures. The resulting array \(S = [s^{0}, \dots, s^{N-1}]\) is the basis on which the Hamiltonian is projected. We then define the initial one-hot vector of length \(N\) as the input to the LOBPCG function.

Let \(J\) be the number of distinct X signatures in the Hamiltonian, and \(K^{(j)}\) be the number of Z signatures and coefficients associated with the \(j\) th X signature. The matrix-vector operation to be passed to the solver acts on the length-\(N\) vector \(v\) of coefficients as

\[v' = \sum_{j=1}^{J} \left( \sum_{k=1}^{K^{(j)}} \alpha^{(j,k)} (-i)^{x^{(j)}z^{(j,k)}} D[z^{(j,k)}] \right) \circ B[x^{(j)}](v).\]

The operation

\[w = B[x](v)\]

consists of the following steps:

Compute the source state \(t^{i} \leftarrow s^{i} \oplus x\) of \(w^{i}\).
If a source index \(j^{i}\) exists such that \(s^{j^{i}} = t^{i}\), \(w^{i} \leftarrow v^{j^{i}}\). Otherwise \(w^{i} \leftarrow 0\).

The operation \(D[z]\) is a diagonal operation that applies a sign factor to each vector entry:

\[D[z](w^{i}) = (-1)^{zs^{i}} w^{i}.\]

Caching

In the expressions above, source indices \([j^{i}]\) and sign factors \([(-1)^{zs^{i}}]\) do not depend on the coefficient vector \(v\) and can be determined once \(S\) is given. In fact, the composition of the sign factors with the coefficients

\[C^{(j)} = \sum_{k=1}^{K^{(j)}} \alpha^{(j,k)} (-i)^{x^{(j)}z^{(j,k)}} [(-1)^{z^{(j,k)}s^{i}}]\]

is entirely static in the same way. It is therefore possible to consider caching these vectors and reusing them in the repeated call to the matrix-vector function. There is however a tradeoff between the compute time and memory footprint, as is always the case with caching.

Concretely, caching the source indices \([j^{i}]\) requires \(4 J N\) bytes of memory, assuming \(N \leq 2^{31}\) (which is actually the hard limit set by other constraints; see the next section) and therefore 32-bit (4-byte) integers are used for vector indexing. Caching the sign bits will require \(\kappa N\) bytes, where

\[\kappa = \lceil \frac{\sum_{j} K^{(j)}}{8} \rceil\]

is the number of bytes required to pack the sign bits for each vector entry. This “dense” packing however may cause inefficiencies in computation that defeats the purpose of caching. For a more straightforward caching, the memory requirement is inflated to \(\bar{\kappa} J N\) bytes, where

\[\bar{\kappa} = \lceil \frac{\max_{j} K^{(j)}}{8} \rceil.\]

If the entire composition of the coefficients are cached instead of the sign bits, \(8 J N\) or \(16 J N\) bytes are used, depending on whether the vector is real or complex (i.e., if there are terms in the Hamiltonian with odd numbers of Ys).

Given that identifying \([j^{i}]\) is an expensive operation, while the other diagonal operations are not, the default behavior of this function is to cache only the source indices. Note however that there are cases when further caching can actually be advantageous also in terms of memory. This is because \(S\) will not be used after caching both the source indices and sign bits / coefficient sums. Since \(S\) occupies \(\lceil n/8 \rceil N\) bytes of memory, caching setting should be adjusted according to the values of \(n\) and \(\{K^{(j)}\}_j\).

Distributed arrays and scaling limits

When the SQD function is called within a context where the global mesh is set via jax.set_mesh(mesh), the state vector is distributed (sharded) among the devices in the mesh, and accordingly all arrays with an axis with size \(N\) follow the same sharding. Even the most aggressive caching strategy described above will be possible this way.

However, there is a limit to scaling in \(N\) (SQD subspace dimension) imposed by the need to sort the states list during the initial uniquification, and also whenever the source indices for an X signature is computed. At the moment, sorting must take place within a single device, with at most \(2^32\) elements involved. Furthermore, source indices identification sorts through a stack of two state lists. Therefore, the maximum achievable \(N\) is \(2^31\). A comparable limit is set by the GPU memory, which is at most O(100)GB per device as of mid-2026.

When the source indices are cached but neither the sign bits nor the diagonals are, the state list \(S\) will also be sharded after the computation of the source indices are done.

SQD API

rqutils.sqd.sqd(hamiltonian: HamiltonianInput, states: StateList, states_size: int | None = None, return_eigvec: bool = True, cache_level: tuple[int, int] = (1, 0)) → float | tuple[float, Vector, StateList]

Perform a sample-based quantum diagonalization of the Hamiltonian.

The Hamiltonian can be given in three different forms:

A tuple of two lists, where the first list enumerates the Pauli strings \(Q\) as strs and the second contains the coefficients \(\alpha\).
Qiskit SparsePauliOp
PauliSumXZ (From rqutils.paulis.symplectic)

States must have binary values and can be passed as an array of integers or booleans.

Internally, the states are bit-packed and represented by \(\lceil (n+1)/8 \rceil\) uint8 s, where the extra bit is placed at position 0 and serves as the indicator for spurious (fill-in) entries. Accordingly, the PauliSumXZ representation of the Hamiltonian must be made with add_padding=True.

Cache level is a 2-tuple where the first element specifies the caching of the source indices (0=no caching, 1=cached) and the second specifies the caching of the diagonal elements (0=no caching, 1=cache sign bits, 2=cache diagonals).

Parameters:

hamiltonian – Hamiltonian to be projected and diagonalized.
states – Binary array of computational basis states to project the Hamiltonian onto. Shape (subspace_dim, num_qubits).
states_size – Fix the size of the states array used in computation to the specified value so that compilation is not triggered at each call with slightly different array sizes.
return_eigvec – Whether to return the eigenvector (coefficients and unique state bitstrings).
cache_level – Switches for caching the results of source indices and sign bits / diagonals. See the module documentation for the detailed discussion of the resource tradeoff involved.

Returns:

Calculated ground state energy, or a tuple of energy, ground state vector, and sorted uniquified states (if return_eigvec=True).

rqutils.sqd.hproj(hamiltonian: HamiltonianInput, states: StateList, unique_states: bool = False) → csr_array

Return the Hamiltonian projected onto the given subspace.

The Hamiltonian can be given in three different forms:

A tuple of two lists, where the first list enumerates the Pauli strings \(Q\) as strs and the second contains the coefficients \(lpha\).
Qiskit SparsePauliOp
PauliSumXZ (From rqutils.paulis.symplectic)

States must have binary values and can be passed as an array of integers or booleans.

Parameters:

hamiltonian – Hamiltonian to be projected and diagonalized.
states – Binary array of computational basis states to project the Hamiltonian onto. Shape (subspace_dim, num_qubits).
unique_states – Whether the states can be assumed to be already uniquified and sorted.

Returns:

The projected Hamiltonian as a sparse matrix.