Sample-Based Krylov Quantum Diagonalization (SKQD)

Build a Krylov subspace by repeatedly time-evolving and sampling a state, then diagonalize classically — scaling to many qubits with sparse linear algebra.

The problem it solves

Finding the ground-state energy of a many-body Hamiltonian means extracting the lowest eigenvalue of a matrix whose dimension is 2^N — hopelessly large to diagonalize directly. The classical workhorse for this is the Krylov subspace method (Lanczos): rather than touch the full matrix, build a small subspace {|ψ⟩, H|ψ⟩, H²|ψ⟩, …} from a reference state and diagonalize H restricted to it. Krylov spaces converge to the extremal (lowest and highest) eigenvalues remarkably fast, so a handful of basis vectors usually suffices.

The catch is that applying powers of H — Hk|ψ⟩ — is not something a quantum computer does naturally. What it does do naturally is real-time evolution e−iHt, implemented by Trotterization. The two are closely related: e−iHt expands as a power series in H, so the time-evolved states {e−iH·k·dt|ψ⟩} at several times capture the spectral information that the Hk|ψ⟩ exploration probes — but one does not explicitly construct each Hk|ψ⟩ term by term, and the numerical stability and the correlation between basis vectors depend on the time step and the chosen sampling times. So you can build a Krylov-type basis purely from circuits the device can run.

Intuition — a Krylov basis from time evolution and sampling

Each Trotterized evolution nudges the reference state further along its own dynamics, exploring different parts of the Hilbert space. Measure each evolved state in the computational basis and you collect the configurations (bitstrings) that carry significant amplitude at that evolution time. Taking the union of these dominant configurations across all the sampled times yields a compact set that spans a subspace rich in the ground state — the quantum device proposes the basis, and nothing more is asked of it.

Generating Krylov states

The loop below evolves the reference state for k Trotter steps, samples it, and adds the observed bitstrings to a growing pool of basis configurations.

import qalgora

# repeatedly evolve and sample, keeping only the dominant, valid configurations
configs = set()
for k in range(krylov_dim):
    counts = qalgora.sample(trotter_evolve, k, dt, shots_count=5000)
    # keep the top-k most frequent bitstrings above a frequency threshold,
    # restricted to the valid symmetry sector (e.g. fixed particle number)
    ranked = sorted(counts.items(), key=lambda kv: kv[1], reverse=True)
    for bitstring, freq in ranked[:top_k]:
        if freq >= freq_threshold and in_valid_sector(bitstring):
            configs.add(bitstring)
    if len(configs) >= max_subspace_dim:   # cap the subspace size
        break

Solving the projected eigenproblem

It helps to separate two distinct subspaces. (A) The time-evolved Krylov basis {e−iH·k·dt|ψ0⟩} is in general non-orthogonal, so one solves a generalized eigenproblem H c = E S c; here both H and the overlap matrix S must be obtained from quantum overlap / matrix-element measurements. (B) The sampled-bitstring basis — the union of computational-basis configurations collected at several evolution times — is treated as an orthogonal computational basis, so S = I and the problem collapses to an ordinary sparse eigenvalue problem; this variant is closer to a Krylov-guided selected CI, and its matrix elements Hij are constructible classically from the bitstrings. The code below uses variant (B), solving the S = I problem with a Lanczos-type routine such as eigsh.

Sparse classical diagonalization

from scipy.sparse.linalg import eigsh

H_sparse = build_sparse_hamiltonian(list(configs))   # CSR matrix on GPU
ground_energy = eigsh(H_sparse, k=1, which="SA")[0][0]
print("SKQD ground-state energy:", ground_energy)

Conditioning and noise — the central caveat

The price of fast Krylov convergence is ill-conditioning. As the Krylov dimension grows, the time-evolved states become nearly linearly dependent — they start to point in almost the same direction — so the overlap matrix S becomes badly conditioned. Once you estimate the matrix elements by sampling, statistical noise on those entries gets amplified by the inverse of a near-singular S, which can throw up spurious eigenvalues below the true ground state — these below-ground-state pseudo-eigenvalues are an artifact of a noisy, ill-conditioned S, not the result of an exact variational diagonalization over the subspace. Practical SKQD therefore caps the Krylov dimension, regularizes or thresholds the overlap matrix, and takes enough shots to keep the matrix-element error small. Trotter error in the evolution adds a second, separate source of bias.

Seeing it in code

Read the two snippets against the theory. In the first, the loop over k performs the repeated time-evolve-and-sample step, and the filtering — top-k by frequency, a frequency threshold, a valid-sector check, and a cap on the subspace size — selects the dominant configurations that define the subspace, rather than blindly accumulating every sampled bitstring. In the second, build_sparse_hamiltonian(list(configs)) projects H onto that subspace, and eigsh(..., k=1, which="SA") is the sparse Lanczos solve that returns the lowest ("smallest algebraic") eigenvalue.