Lanczos Algorithm#

lanczos_alg(H: BlockEncoding | QubitOperator, D: int, operand_prep: Callable[[...], Any], mes_kwargs: Dict[str, object] = {}, cutoff: float = 0.01, show_info: bool = False)[source]#

Estimate the ground state energy of a Hamiltonian using the Exact and efficient Lanczos method on a quantum computer.

This function implements the algorithm proposed in Kirby et al. by constructing a Krylov subspace using block-encodings of Chebyshev polynomials \(T_k(H)\), bypassing the need for real or imaginary time evolution.

The quantum Lanczos algorithm efficiently constructs a Krylov subspace by applying Chebyshev polynomials of the Hamiltonian \(T_k(H)\) to an initial state \(\ket{\psi_0}\). The Krylov space dimension \(D\) determines the accuracy of ground state energy estimation, with convergence guaranteed when the initial state has a sufficiently large overlap (\(|\gamma_0|=\Omega(1/\text{poly}(n))\) for \(n\) qubits) with the true ground state. The Chebyshev approach allows exact Krylov space construction (up to sample noise) without real or imaginary time evolution. This algorithm is motivated by the rapid convergence of the Lanczos method for estimating extremal eigenvalues, and its quantum version avoids the classical barrier of exponential cost in representing Krylov vectors.

The Krylov basis vectors are generated by applying Chebyshev polynomials of the Hamiltonian to \(\ket{\psi_0}\):

\[\ket{\psi_k}=T_k(H)\ket{\psi_0}, k=0,\dots,D-1,\]

which span the same Krylov space as \(\{H^k\ket{\psi_0}\}\). The generalized eigenvalue problem \(\mathbf{H}\vec{v}=\epsilon \mathbf{S}\vec{v}\) arises as a standard quantum subspace diagonalization task, with the projected Hamiltonian matrix \(\mathbf{H}\) and overlap matrix \(\mathbf{S}\) defined in this Chebyshev basis.

We can obtain \(\mathbf{S}\) by evaluating its matrix elements as

\[S_{ij}=\bra{\psi_0}T_i(H)T_j(H)\ket{\psi_0}=\langle T_i(H)T_j(H)\rangle_0,\]

where \(\langle\cdot\rangle_0\) denotes expectation value with respect to the initial state \(\ket{\psi_0}\).

Using the identity \(T_i(x)T_j(x)=\frac{1}{2}(T_{i+j}(x)+T_{|i-j|}(x))\), the overlap elements become

\[S_{ij}=\frac{1}{2}\langle T_{i+j}(H)+T_{|i-j|}(H)\rangle_0=\frac{1}{2}\bigg(\langle T_{i+j}(H)\rangle_0+\langle T_{|i-j|}(H)\rangle_0\bigg).\]

For \(\mathbf{H}\), we have \(H_{ij}=\langle T_i(H) H T_j(H)\rangle_0\). By using the fact that \(H=T_1(H)\) and applying the Chebyshev identity twice, we get

\[H_{ij}=\frac{1}{4}\bigg(\langle T_{i+j+1}(H)\rangle_0+\langle T_{|i+j-1|}(H)\rangle_0 + \langle T_{|i-j+1|}(H)\rangle_0 + \langle T_{|i-j-1|}(H)\rangle_0\bigg).\]

Because \(i,j=0,1,2,\dots, D-1\), all matrix elements \(S_{ij}\) and \(H_{ij}\) are expressed as linear combinations of expectation values of Chebyshev polynomials with respect to the initial state \(\langle T_k(H)\rangle_0\), where \(k=0,1,2,\dots, 2D-1\). The highest value \(2D-1\) comes from the first term for \(H_{ij}\) when \(i=j=D-1\). To construct \(\mathbf{H}\) and \(\mathbf{S}\) it is enough to estimate all of the expectation values \(\langle T_k(H)\rangle_0\).

The final step in this implementation is solving the generalized eigenvalue problem \(\mathbf{H}\vec{v}=\epsilon \mathbf{S}\vec{v}\). In practice, the overlap matrix \(\mathbf{S}\) can become ill-conditioned due to the linear dependencies in the Krylov basis or sampling noise in the expectation values \(T_k(H)\). To prevent numerical instability, the matrices are regularized via thresholding, a process involving discarding small eigenvalues of \(\mathbf{S}\) below a specified cutoff. The choice of this cutoff is vital; a value too small may fail to suppress noise, while a value too large may discard physically relevant information.

The entire approach can be summarized by the following steps:

Run quantum Lanczos subroutine lanczos_expvals() to obtain Chebyshev expectation values \(\langle T_k(H)\rangle_0\).
Build overlap and Hamiltonian subspace matrices \((\mathbf{S}, \mathbf{H})\).
Regularize overlap matrix \(\mathbf{S}\) and \(\mathbf{H}\) by projecting onto the subspace with well conditioned eigenvalues.
Solve generalized eigenvalue problem \(\mathbf{H}\vec{v}=\epsilon \mathbf{S}\vec{v}\).
Return lowest eigenvalue \(\epsilon_{\text{min}}\) as ground state energy estimate.

Parameters:

HQubitOperator or BlockEncoding: Hamiltonian for which to estimate the ground-state energy. If a QubitOperator is provided, it is automatically converted to a Pauli block-encoding.
Dint: Krylov space dimension.
operand_prepcallable: Function returning the (operand) QuantumVariables in the initial system state \(\ket{\psi_0}\), i.e., operands=operand_prep(). Must return a QuantumVariable or a tuple of QuantumVariables.
mes_kwargsdict: The keyword arguments for the measurement function. By default, 100_000 shots are executed for measuring each expectation value.
cutofffloat: Regularization cutoff threshold for overlap matrix \(\mathbf{S}\). The default is 1e-2.
show_infobool, optional: If True, a dictionary with detailed information is returned. The default is False.

Returns:

energyfloat

Estimated ground state energy of the Hamiltonian H.

infodict, optional

Full details including:

‘energy’float
Ground-state energy estimate
‘eigvals’ArrayLike, shape (D,)
Eigenvalues of regularized problem
‘eigvecs’ArrayLike, shape (D, D)
Eigenvectors of regularized problem
‘H’ArrayLike, shape (D, D)
The Hamiltonian matrix
‘H_reg’ArrayLike, shape (D, D)
Regularized Hamiltonian matrix
‘S’ArrayLike, shape (D, D)
The overlap matrix
‘S_reg’ArrayLike, shape (D, D)
Regularized overlap matrix
‘Tk_expvals’ArrayLike, shape (2D,)
Chebyshev expectation values

Examples

Example 1: Jasp Mode (Dynamic Execution)

This mode uses Qrisp’s Jasp framework for JIT-compilation and tracing, ideal for high-performance execution.

from qrisp import QuantumVariable
from qrisp.algorithms.lanczos import lanczos_alg
from qrisp.operators import X, Y, Z
from qrisp.vqe.problems.heisenberg import create_heisenberg_init_function
from qrisp.jasp import jaspify
import networkx as nx

# Define a 1D Heisenberg model
L = 6
G = nx.cycle_graph(L)
H = (1/4)*sum((X(i)*X(j) + Y(i)*Y(j) + 0.5*Z(i)*Z(j)) for i,j in G.edges())

# Prepare initial state function (tensor product of singlets)
M = nx.maximal_matching(G)
U_singlet = create_heisenberg_init_function(M)

def operand_prep():
    qv = QuantumVariable(H.find_minimal_qubit_amount())
    U_singlet(qv)
    return qv

D = 6  # Krylov dimension

@jaspify(terminal_sampling=True)
def main():
    return lanczos_alg(H, D, operand_prep, show_info=True)

energy, info = main()
print(f"Ground state energy estimate: {energy}")

We can compare the results obtained by classical calculation:

print(f"Ground state energy: {H.ground_state_energy()}")

Example 2: Standard Mode (Static Execution)

The standard mode is useful for simple scripts and direct hardware interface without JAX overhead.

# Using the same Hamiltonian and prep function from the previous example
energy = lanczos_alg(H, D=6, operand_prep=operand_prep)
print(f"Ground state energy: {energy}")
print(f"Ground state energy: {H.ground_state_energy()}")

`lanczos_even`(BE, k, operand_prep)	This function implements the Krylov space construction via block-encodings of Chebyshev polynomials \(T_k(H)\), following the layout in Figure 1(a) of Kirby et al..
`lanczos_odd`(BE, k, operand_prep)	This function implements the Krylov space construction via block-encodings of Chebyshev polynomials \(T_k(H)\), following the layout in Figure 1(b) of Kirby et al..
`lanczos_expvals`(H, D, operand_prep[, mes_kwargs])	Estimate the expectation values of Chebyshev polynomials \(\langle T_k(H) \rangle_0\) for the exact and efficient Quantum Lanczos method.
`build_S_H_from_Tk`(expvals)	Construct the overlap matrix \(\mathbf{S}\) and the Krylov Hamiltonian matrix \(\mathbf{H}\) from Chebyshev polynomial expectation values.
`regularize_S_H`(S, H_mat[, cutoff])	Regularize the overlap matrix \(\mathbf{S}\) by retaining only eigenvectors with sufficiently large eigenvalues and project the Hamiltonian matrix \(\mathbf{H}\) accordingly.
`generalized_eigh`(A, B)	Solves the generalized eigenvalue problem \(A v = \lambda B v\) for a complex Hermitian or real symmetric matrix \(A\) and a real symmetric positive-definite matrix \(B\).