From d9dd3e4f58f71e1a225e79c937ddae3424ba8a23 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Wed, 6 Aug 2025 17:42:24 +0200 Subject: [PATCH 01/32] Implemented the quantum part of the Lanczos method without the postprocessing to construct H and S DxD matrices --- src/qrisp/algorithms/lanczos.py | 47 +++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 src/qrisp/algorithms/lanczos.py diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py new file mode 100644 index 000000000..3045fe64b --- /dev/null +++ b/src/qrisp/algorithms/lanczos.py @@ -0,0 +1,47 @@ +from qrisp.operators import * +from qrisp import * +from qrisp.algorithms.grover import* +import numpy as np + +def Lanczos(H, D): + + unitaries, coeffs = H.unitaries() + + num_unitaries = len(unitaries) + + def state_prep(case): + prepare(case, np.sqrt(coeffs)) + + n = np.int64(np.ceil(np.log2(num_unitaries))) + + + + def UR(case_indicator, operand, unitaries): + qswitch(operand, case_indicator, unitaries) + diffuser(case_indicator, state_function = state_prep) + for k in jrange(0, 2*D): + case_indicator = QuantumFloat(n) + operand = QuantumVariable(2) + if k % 2 == 0: + with conjugate(state_prep)(case_indicator): + for _ in jrange(k//2): + UR(case_indicator, operand, unitaries) + print(case_indicator.qs) + multi_measurement([case_indicator, operand]) + + else: + state_prep(case_indicator) + for _ in jrange(k//2): + UR(case_indicator, operand, unitaries) + qv = QuantumVariable(1) + h(qv) + with control(qv[0]): + qswitch(operand, case_indicator, unitaries) + h(qv) + print(case_indicator.qs) + qv.get_measurement() + qv.delete() + + case_indicator.delete() + operand.delete() + From ae61fe864b6eeef7b58473e03809d3e023c5b58c Mon Sep 17 00:00:00 2001 From: mat70593 Date: Wed, 6 Aug 2025 17:46:22 +0200 Subject: [PATCH 02/32] removed multimeasurement for k=even based on the paper --- src/qrisp/algorithms/lanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 3045fe64b..b64867061 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -27,7 +27,7 @@ def UR(case_indicator, operand, unitaries): for _ in jrange(k//2): UR(case_indicator, operand, unitaries) print(case_indicator.qs) - multi_measurement([case_indicator, operand]) + case_indicator.get_measurement() else: state_prep(case_indicator) From b7157844201613a013e1a649a1e6b0aca5869abc Mon Sep 17 00:00:00 2001 From: mat70593 Date: Thu, 21 Aug 2025 16:06:25 +0200 Subject: [PATCH 03/32] added docstrings and classical post-processing functions, as well as create the lanczos_alg function executing the subroutine, returning the ground state energy --- src/qrisp/algorithms/lanczos.py | 248 ++++++++++++++++++++++++++++---- 1 file changed, 223 insertions(+), 25 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index b64867061..8b37b06f9 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -2,46 +2,244 @@ from qrisp import * from qrisp.algorithms.grover import* import numpy as np +import scipy -def Lanczos(H, D): +def inner_lanczos(H, D, state_prep_func): + """ + Perform the quantum subroutine of the exact and efficient Lanczos method to estimate Chebyshev polynomials of a Hamiltonian. - unitaries, coeffs = H.unitaries() + This function implements the Krylov space construction via block-encodings + of Chebyshev polynomials T_k(H), as described in "Exact and efficient Lanczos method + on a quantum computer" (arXiv:2208.00567). + + At each order k = 0, ..., 2D-1 it prepares and measures circuits corresponding + either to $\bra{\psi\lfloor k/2\rfloor}R\ket{\psi\lfloor k/2\rfloor}$ for even k, or + $\bra{\psi\lfloor k/2\rfloor U\ket{\psi\lfoor k/2\rfloor}$ for odd k. + The measured statistics encode the expectation values $⟨T_k(H)⟩$. - num_unitaries = len(unitaries) + Parameters + ---------- + H : Hamiltonian (QubitOperator) + Hamiltonian represented as a Pauli sum operator with block-encoding accessible via `unitaries()`. + D : int + Krylov space dimension. Determines maximum Chebyshev order (2D-1). + state_prep_func : callable + Function preparing the initial system state (operand) $\ket{\psi_0}$ on a QuantumVariable. - def state_prep(case): - prepare(case, np.sqrt(coeffs)) + Returns + ------- + meas_res : dictionary + Measurement results for each Chebyshev polynomial order, suitable for conversion into expectation values. + """ + # Extract unitaries and size of the case_indicator QuantumFloat. + unitaries, coeffs = H.unitaries() + num_unitaries = len(unitaries) n = np.int64(np.ceil(np.log2(num_unitaries))) - + # Prepare $\ket{G} = \sum_i \sqrt{\alpha_i}\ket{i}$ on the auxiliary register (case_indicator) (paper Def. 1, Eq. (4)) + def state_prep(case): + prepare(case, np.sqrt(coeffs)) def UR(case_indicator, operand, unitaries): - qswitch(operand, case_indicator, unitaries) - diffuser(case_indicator, state_function = state_prep) + qswitch(operand, case_indicator, unitaries) #qswitch applies $U = \sum_i\ket{i}\bra{i}\otimes P_i$. + diffuser(case_indicator, state_function=state_prep) # diffuser implements the reflection operator R about $\ket{G}$. + + meas_res = {} + for k in jrange(0, 2*D): case_indicator = QuantumFloat(n) - operand = QuantumVariable(2) + operand = QuantumVariable(H.find_minimal_qubit_amount()) + state_prep_func(operand) # prepare operand $\ket{\psi_0}$ + if k % 2 == 0: + # EVEN k: Figure 1 top with conjugate(state_prep)(case_indicator): for _ in jrange(k//2): UR(case_indicator, operand, unitaries) - print(case_indicator.qs) - case_indicator.get_measurement() + + meas = case_indicator.get_measurement() + meas_res[k] = meas + + else: + # ODD k: Figure 1 bottom + state_prep(case_indicator) + for _ in jrange(k//2): + UR(case_indicator, operand, unitaries) + qv = QuantumVariable(1) + h(qv) # Hadamard test for + with control(qv[0]): + qswitch(operand, case_indicator, unitaries) # control-U on the case_indicator QuantumFloat + h(qv) # Hadamard test for + meas = qv.get_measurement() + meas_res[k] = meas - else: - state_prep(case_indicator) - for _ in jrange(k//2): - UR(case_indicator, operand, unitaries) - qv = QuantumVariable(1) - h(qv) - with control(qv[0]): - qswitch(operand, case_indicator, unitaries) - h(qv) - print(case_indicator.qs) - qv.get_measurement() - qv.delete() + case_indicator.delete() + operand.delete() + + return meas_res + +def compute_expectation(counts): + """ + Convert measurement counts into an expectation value. + + Assumes measurement outcomes correspond to ±1 eigenvalues of observables + (reflection R or block-encoding unitary U). + + For even k we measure the auxilary case_indicator QuantumFloat in the computational basis. We then map: + +1 if all-zeroes outcome (projector 2|0><0| - I), else -1 (paper step 2). + + For odd k we perform the hadamard test for the unitary on the controlled unitary and then measure the + auxilary case_indicator QuantumFloat in the Z basis + + Parameters + ---------- + counts : dictionary + Measurement counts or probabilities, keyed by outcomes. + + Returns + ------- + expval : float + Expectation value of the measured observable. + """ + expval = 0.0 + for outcome, prob in counts.items(): + + if int(outcome) == 0: + expval += prob * 1 + else: + expval += prob * (-1) + return expval + +def build_S_H_from_Tk(Tk_expectation, D): + """ + Construct Lanczos overlap matrix S and Hamiltonian matrix H from Chebyshev polynomials. + + Uses Chebyshev recurrence identities to compute matrix elements + (Eq. (17), (19) in the reference paper) from measured expectations. + + Parameters + ---------- + Tk_expectation : dictionary + Mapping k to ⟨T_k(H)⟩. + + D : integer + Krylov space dimension. + + Returns + ------- + S : ndarray + overlap (Gram) matrix S or Krylov states. + + H_mat : ndarray + Hamiltonian matrix in Krylov subspace. + """ + def Tk(k): + k = abs(k) + return Tk_expectation.get(k, 0) + + S = np.zeros((D, D)) + H_mat = np.zeros((D, D)) + + for i in range(D): + for j in range(D): + S[i, j] = 0.5 * (Tk(i + j) + Tk(abs(i - j))) + H_mat[i, j] = 0.25 * ( + Tk(i + j + 1) + + Tk(abs(i + j - 1)) + + Tk(abs(i - j + 1)) + + Tk(abs(i - j - 1)) + ) + return S, H_mat + +def regularize_S_H(S, H_mat, cutoff=1e-3): + """ + Regularize overlap matrix S by thresholding eigenvalues below cutoff * max_eigenvalue. + Project both S and H onto subspace defined by large eigenvalues. + + Parameters + ---------- + S : Overlap matrix + H_mat : Hamiltonian matrix + cutoff : float + Eigenvalue threshold for regularizing S + + Returns + ------- + S_reg : ndarray + regularized overlap matrix S. + H_reg : ndarray + regularized Hamiltonian matrix in Krylov subspace. + """ + eigvals, eigvecs = np.linalg.eigh(S) + max_eigval = eigvals.max() + threshold = cutoff * max_eigval + mask = eigvals > threshold + U = eigvecs[:, mask] + S_reg = U.T @ S @ U + H_reg = U.T @ H @ U + return S_reg, H_reg + +def lanczos_alg(H, D, state_prep_func, cutoff=1e-2): + """ + Exact and efficient Lanczos method on a quantum computer for ground state energy estimation. + + Implements the following steps: + 1. Run quantum Lanczos subroutine to obtain Chebyshev expectation values. + 2. Build overlap and Hamiltonian subspace matrices (S, H). + 3. Regularize overlap matrix S and H. + 4. Solve generalized eigenvalue problem Hv=ESv. + 5. Return lowest eigenvalue as ground state energy estimate. + + Parameters + ---------- + H : Hamiltonian (QuantumOperator) + Hamiltonian operator expressed as the Pauli sum. + D : int + Krylov space dimension. + state_prep_func : callable + Function preparing the initial system state (operand) $\ket{\psi_0}$ on a QuantumVariable. + cutoff : float + Regularization cutoff threshold for overlap matrix S. + + Returns + ------- + energy : float + Estimated ground state energy of the Hamiltonian H. + results : dictionary + Full details including: + - 'Tk_expvals': dict of Chebyshev expectation values + - 'energy': ground state estimate + - 'eigvals': eigenvalues of regularized problem + - 'eigvecs': eigenvectors + - 'S_reg': regularized overlap matrix + - 'H_reg': regularized Hamiltonian matrix + """ + unitaries, coeffs = H.unitaries() - case_indicator.delete() - operand.delete() + # Step 1: Quantum Lanczos: Get expectation values of Chebyshev polynomials + meas_counts = inner_lanczos(H, D, state_prep_func) + + # Step 2: Convert counts to expectation values + Tk_expvals = {k: compute_expectation(counts) for k, counts in meas_counts.items()} + + # Step 3: Build matrices S and H + S, H_mat = build_S_H_from_Tk(Tk_expvals, D) + + # Step 4: Regularize matrices via thresholding + S_reg, H_reg = regularize_S_H(S, H_mat, cutoff=cutoff) + + # Step 5: Solve generalized eigenvalue problem Hv = ESv + evals, evecs = scipy.linalg.eigh(H_reg, S_reg) + + ground_state_energy = np.min(evals)*sum(coeffs) + results = { + 'Tk_expvals': Tk_expvals, + 'energy': ground_state_energy, + 'eigvals': evals, + 'eigvecs': evecs, + 'S_reg': S_reg, + 'H_reg': H_reg, + } + return ground_state_energy, results \ No newline at end of file From f5d30e7151db8a33f765107bb83aba425a38e820 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Thu, 21 Aug 2025 16:16:53 +0200 Subject: [PATCH 04/32] added example to test the lanczos_alg functionality --- src/qrisp/algorithms/lanczos.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 8b37b06f9..0d5a6837d 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -214,6 +214,32 @@ def lanczos_alg(H, D, state_prep_func, cutoff=1e-2): - 'eigvecs': eigenvectors - 'S_reg': regularized overlap matrix - 'H_reg': regularized Hamiltonian matrix + + Examples + -------- + + :: + + from qrisp.vqe.problems.heisenberg import create_heisenberg_init_function + import networkx as nx + + L = 6 + G = nx.Graph() + G.add_edges_from([(k, (k+1) % L) for k in range(L - 1)]) + + # Define Hamiltonian e.g. Heisenberg with custom couplings + 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) + + D = 6 # Krylov dimension + energy, info = lanczos_alg(H, D, U_singlet) + + print(f"Ground state energy estimate: {energy}") + + """ unitaries, coeffs = H.unitaries() From 06ec68bf884430760f457371adcb747cfbb57a11 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Fri, 22 Aug 2025 14:23:29 +0200 Subject: [PATCH 05/32] Improved docstrings --- src/qrisp/algorithms/lanczos.py | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 0d5a6837d..0a57c5a6c 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -184,12 +184,25 @@ def lanczos_alg(H, D, state_prep_func, cutoff=1e-2): """ Exact and efficient Lanczos method on a quantum computer for ground state energy estimation. + This function implements the Lanczos method on a quantum computer using block-encodings of Chebyshev + polynomials T_k(H), closely following the algorithm proposed in + "Exact and efficient Lanczos method on a quantum computer" (arXiv:2208.00567). + + The quantum Lanczos algorithm efficiently constructs a Krylov subspace by applying Chebyshev polynomials + of the Hamiltonian to an initial state. 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. + Implements the following steps: - 1. Run quantum Lanczos subroutine to obtain Chebyshev expectation values. + 1. Run quantum Lanczos subroutine to obtain Chebyshev expectation values $\langle T_k(H)\rangle$. 2. Build overlap and Hamiltonian subspace matrices (S, H). - 3. Regularize overlap matrix S and H. - 4. Solve generalized eigenvalue problem Hv=ESv. - 5. Return lowest eigenvalue as ground state energy estimate. + 3. Regularize overlap matrix S and H by projecting onto the subspace with well conditioned eigenvalues. + 4. Solve generalized eigenvalue problem $\mathbf{H}\vec{v}=\epsilon\mathbf{S}\vec{v}$. + 5. Return lowest eigenvalue $\epsilon_{\text{min}}$ as ground state energy estimate. Parameters ---------- @@ -239,23 +252,23 @@ def lanczos_alg(H, D, state_prep_func, cutoff=1e-2): print(f"Ground state energy estimate: {energy}") - + """ unitaries, coeffs = H.unitaries() # Step 1: Quantum Lanczos: Get expectation values of Chebyshev polynomials meas_counts = inner_lanczos(H, D, state_prep_func) - # Step 2: Convert counts to expectation values + # Convert counts to expectation values Tk_expvals = {k: compute_expectation(counts) for k, counts in meas_counts.items()} - # Step 3: Build matrices S and H + # Step 2: Build matrices S and H S, H_mat = build_S_H_from_Tk(Tk_expvals, D) - # Step 4: Regularize matrices via thresholding + # Step 3: Regularize matrices via thresholding S_reg, H_reg = regularize_S_H(S, H_mat, cutoff=cutoff) - # Step 5: Solve generalized eigenvalue problem Hv = ESv + # Step 4: Solve generalized eigenvalue problem $\mathbf{H}\vec{v}=\epsilon\mathbf{S}\vec{v}$ evals, evecs = scipy.linalg.eigh(H_reg, S_reg) ground_state_energy = np.min(evals)*sum(coeffs) From 2ef426b2ca48c32bec3a7ba46be563ec3746ad4a Mon Sep 17 00:00:00 2001 From: mat70593 Date: Fri, 22 Aug 2025 14:26:27 +0200 Subject: [PATCH 06/32] Added mes_kwargs as a keyword argument to inner_lanczos and lanczos_alg --- src/qrisp/algorithms/lanczos.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 0a57c5a6c..babc1b960 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -4,7 +4,7 @@ import numpy as np import scipy -def inner_lanczos(H, D, state_prep_func): +def inner_lanczos(H, D, state_prep_func, mes_kwargs): """ Perform the quantum subroutine of the exact and efficient Lanczos method to estimate Chebyshev polynomials of a Hamiltonian. @@ -25,6 +25,8 @@ def inner_lanczos(H, D, state_prep_func): Krylov space dimension. Determines maximum Chebyshev order (2D-1). state_prep_func : callable Function preparing the initial system state (operand) $\ket{\psi_0}$ on a QuantumVariable. + mes_kwargs : dictionary + The keyword arguments for the measurement function. Returns ------- @@ -58,7 +60,7 @@ def UR(case_indicator, operand, unitaries): for _ in jrange(k//2): UR(case_indicator, operand, unitaries) - meas = case_indicator.get_measurement() + meas = case_indicator.get_measurement(**mes_kwargs) meas_res[k] = meas else: @@ -71,7 +73,7 @@ def UR(case_indicator, operand, unitaries): with control(qv[0]): qswitch(operand, case_indicator, unitaries) # control-U on the case_indicator QuantumFloat h(qv) # Hadamard test for - meas = qv.get_measurement() + meas = qv.get_measurement(**mes_kwargs) meas_res[k] = meas case_indicator.delete() @@ -177,10 +179,10 @@ def regularize_S_H(S, H_mat, cutoff=1e-3): mask = eigvals > threshold U = eigvecs[:, mask] S_reg = U.T @ S @ U - H_reg = U.T @ H @ U + H_reg = U.T @ H_mat @ U return S_reg, H_reg -def lanczos_alg(H, D, state_prep_func, cutoff=1e-2): +def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): """ Exact and efficient Lanczos method on a quantum computer for ground state energy estimation. @@ -212,6 +214,8 @@ def lanczos_alg(H, D, state_prep_func, cutoff=1e-2): Krylov space dimension. state_prep_func : callable Function preparing the initial system state (operand) $\ket{\psi_0}$ on a QuantumVariable. + mes_kwargs : dictionary + The keyword arguments for the measurement function. cutoff : float Regularization cutoff threshold for overlap matrix S. @@ -257,7 +261,7 @@ def lanczos_alg(H, D, state_prep_func, cutoff=1e-2): unitaries, coeffs = H.unitaries() # Step 1: Quantum Lanczos: Get expectation values of Chebyshev polynomials - meas_counts = inner_lanczos(H, D, state_prep_func) + meas_counts = inner_lanczos(H, D, state_prep_func, mes_kwargs) # Convert counts to expectation values Tk_expvals = {k: compute_expectation(counts) for k, counts in meas_counts.items()} From 601b7b0d73d3e0148b8678f034edd02d63023dd4 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Fri, 22 Aug 2025 15:01:22 +0200 Subject: [PATCH 07/32] added lanczos to the documentation --- documentation/source/reference/Algorithms/Lanczos.rst | 11 +++++++++++ documentation/source/reference/Algorithms/index.rst | 1 + src/qrisp/__init__.py | 1 + src/qrisp/algorithms/__init__.py | 1 + 4 files changed, 14 insertions(+) create mode 100644 documentation/source/reference/Algorithms/Lanczos.rst diff --git a/documentation/source/reference/Algorithms/Lanczos.rst b/documentation/source/reference/Algorithms/Lanczos.rst new file mode 100644 index 000000000..28eda8fd8 --- /dev/null +++ b/documentation/source/reference/Algorithms/Lanczos.rst @@ -0,0 +1,11 @@ +.. _lanczos_alg: + + +Lanczos method +============== + +.. currentmodule:: qrisp.lanzcos + +.. autofunction:: lanczos_alg + +.. autofunction:: inner_lanczos \ No newline at end of file diff --git a/documentation/source/reference/Algorithms/index.rst b/documentation/source/reference/Algorithms/index.rst index e56cce031..e0059c1c5 100644 --- a/documentation/source/reference/Algorithms/index.rst +++ b/documentation/source/reference/Algorithms/index.rst @@ -43,6 +43,7 @@ We encourage you to explore these algorithms, delve into their documentation, an QITE Shor Grover + Lanczos QuantumBacktrackingTree quantum_counting QMCI diff --git a/src/qrisp/__init__.py b/src/qrisp/__init__.py index c4221877f..56c2ff4f8 100644 --- a/src/qrisp/__init__.py +++ b/src/qrisp/__init__.py @@ -45,6 +45,7 @@ "vqe", "qite", "qmci", + "lanczos" ]: sys.modules["qrisp." + i] = sys.modules["qrisp.algorithms." + i] diff --git a/src/qrisp/algorithms/__init__.py b/src/qrisp/algorithms/__init__.py index d09f2a103..57f90fc02 100644 --- a/src/qrisp/algorithms/__init__.py +++ b/src/qrisp/algorithms/__init__.py @@ -25,3 +25,4 @@ import qrisp.algorithms.vqe as vqe import qrisp.algorithms.qite as qite import qrisp.algorithms.qmci as qmci +import qrisp.algorithms.lanczos as lanczos From 66f4525b3cf700b4aee6b615919d14e902c097b6 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Fri, 22 Aug 2025 15:08:38 +0200 Subject: [PATCH 08/32] minor changes to the docstrings --- src/qrisp/algorithms/lanczos.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index babc1b960..7b8d141a6 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -9,20 +9,20 @@ def inner_lanczos(H, D, state_prep_func, mes_kwargs): Perform the quantum subroutine of the exact and efficient Lanczos method to estimate Chebyshev polynomials of a Hamiltonian. This function implements the Krylov space construction via block-encodings - of Chebyshev polynomials T_k(H), as described in "Exact and efficient Lanczos method + of Chebyshev polynomials $T_k(H)$, as described in "Exact and efficient Lanczos method on a quantum computer" (arXiv:2208.00567). - At each order k = 0, ..., 2D-1 it prepares and measures circuits corresponding + At each order $k = 0, \dots, 2D-1$ it prepares and measures circuits corresponding either to $\bra{\psi\lfloor k/2\rfloor}R\ket{\psi\lfloor k/2\rfloor}$ for even k, or $\bra{\psi\lfloor k/2\rfloor U\ket{\psi\lfoor k/2\rfloor}$ for odd k. - The measured statistics encode the expectation values $⟨T_k(H)⟩$. + The measured statistics encode the expectation values $\angleT_k(H)\rangle$. Parameters ---------- H : Hamiltonian (QubitOperator) Hamiltonian represented as a Pauli sum operator with block-encoding accessible via `unitaries()`. D : int - Krylov space dimension. Determines maximum Chebyshev order (2D-1). + Krylov space dimension. Determines maximum Chebyshev order $(2D-1)$. state_prep_func : callable Function preparing the initial system state (operand) $\ket{\psi_0}$ on a QuantumVariable. mes_kwargs : dictionary @@ -85,11 +85,11 @@ def compute_expectation(counts): """ Convert measurement counts into an expectation value. - Assumes measurement outcomes correspond to ±1 eigenvalues of observables + Assumes measurement outcomes correspond to $\pm 1$ eigenvalues of observables (reflection R or block-encoding unitary U). For even k we measure the auxilary case_indicator QuantumFloat in the computational basis. We then map: - +1 if all-zeroes outcome (projector 2|0><0| - I), else -1 (paper step 2). + $+1$ if all-zeroes outcome (projector $2\ket{0}\bra{0} - I), else $-1$ (paper step 2). For odd k we perform the hadamard test for the unitary on the controlled unitary and then measure the auxilary case_indicator QuantumFloat in the Z basis @@ -118,7 +118,7 @@ def build_S_H_from_Tk(Tk_expectation, D): Construct Lanczos overlap matrix S and Hamiltonian matrix H from Chebyshev polynomials. Uses Chebyshev recurrence identities to compute matrix elements - (Eq. (17), (19) in the reference paper) from measured expectations. + (Equations (17), (19) in the reference paper) from measured expectations. Parameters ---------- @@ -187,7 +187,7 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): Exact and efficient Lanczos method on a quantum computer for ground state energy estimation. This function implements the Lanczos method on a quantum computer using block-encodings of Chebyshev - polynomials T_k(H), closely following the algorithm proposed in + polynomials $T_k(H)$, closely following the algorithm proposed in "Exact and efficient Lanczos method on a quantum computer" (arXiv:2208.00567). The quantum Lanczos algorithm efficiently constructs a Krylov subspace by applying Chebyshev polynomials From 4a790e3602e388e2013e46e2b2408633a029e435 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Fri, 22 Aug 2025 15:43:52 +0200 Subject: [PATCH 09/32] fixed typo preventing the documentation to be compiled. Thanks @renezander90 ! --- documentation/source/reference/Algorithms/Lanczos.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/documentation/source/reference/Algorithms/Lanczos.rst b/documentation/source/reference/Algorithms/Lanczos.rst index 28eda8fd8..44c64327b 100644 --- a/documentation/source/reference/Algorithms/Lanczos.rst +++ b/documentation/source/reference/Algorithms/Lanczos.rst @@ -4,7 +4,7 @@ Lanczos method ============== -.. currentmodule:: qrisp.lanzcos +.. currentmodule:: qrisp.lanczos .. autofunction:: lanczos_alg From 8488c0406f6046d9ccf618f084f5996b4478e638 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Fri, 22 Aug 2025 16:10:33 +0200 Subject: [PATCH 10/32] extended documentation for the helper functions --- .../source/reference/Algorithms/Lanczos.rst | 8 ++++- .../qrisp.lanczos.build_S_H_from_Tk.rst | 6 ++++ .../qrisp.lanczos.compute_expectation.rst | 6 ++++ .../generated/qrisp.lanczos.inner_lanczos.rst | 6 ++++ .../qrisp.lanczos.regularize_S_H.rst | 6 ++++ src/qrisp/algorithms/lanczos.py | 30 +++++++++++-------- 6 files changed, 49 insertions(+), 13 deletions(-) create mode 100644 documentation/source/reference/Algorithms/generated/qrisp.lanczos.build_S_H_from_Tk.rst create mode 100644 documentation/source/reference/Algorithms/generated/qrisp.lanczos.compute_expectation.rst create mode 100644 documentation/source/reference/Algorithms/generated/qrisp.lanczos.inner_lanczos.rst create mode 100644 documentation/source/reference/Algorithms/generated/qrisp.lanczos.regularize_S_H.rst diff --git a/documentation/source/reference/Algorithms/Lanczos.rst b/documentation/source/reference/Algorithms/Lanczos.rst index 44c64327b..0b46e7c9a 100644 --- a/documentation/source/reference/Algorithms/Lanczos.rst +++ b/documentation/source/reference/Algorithms/Lanczos.rst @@ -8,4 +8,10 @@ Lanczos method .. autofunction:: lanczos_alg -.. autofunction:: inner_lanczos \ No newline at end of file +.. autosummary:: + :toctree: generated/ + + inner_lanczos + compute_expectation + build_S_H_from_Tk + regularize_S_H \ No newline at end of file diff --git a/documentation/source/reference/Algorithms/generated/qrisp.lanczos.build_S_H_from_Tk.rst b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.build_S_H_from_Tk.rst new file mode 100644 index 000000000..dbb61baca --- /dev/null +++ b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.build_S_H_from_Tk.rst @@ -0,0 +1,6 @@ +qrisp.lanczos.build\_S\_H\_from\_Tk +=================================== + +.. currentmodule:: qrisp.lanczos + +.. autofunction:: build_S_H_from_Tk \ No newline at end of file diff --git a/documentation/source/reference/Algorithms/generated/qrisp.lanczos.compute_expectation.rst b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.compute_expectation.rst new file mode 100644 index 000000000..ad4c215a0 --- /dev/null +++ b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.compute_expectation.rst @@ -0,0 +1,6 @@ +qrisp.lanczos.compute\_expectation +================================== + +.. currentmodule:: qrisp.lanczos + +.. autofunction:: compute_expectation \ No newline at end of file diff --git a/documentation/source/reference/Algorithms/generated/qrisp.lanczos.inner_lanczos.rst b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.inner_lanczos.rst new file mode 100644 index 000000000..99a74022d --- /dev/null +++ b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.inner_lanczos.rst @@ -0,0 +1,6 @@ +qrisp.lanczos.inner\_lanczos +============================ + +.. currentmodule:: qrisp.lanczos + +.. autofunction:: inner_lanczos \ No newline at end of file diff --git a/documentation/source/reference/Algorithms/generated/qrisp.lanczos.regularize_S_H.rst b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.regularize_S_H.rst new file mode 100644 index 000000000..90c1766f4 --- /dev/null +++ b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.regularize_S_H.rst @@ -0,0 +1,6 @@ +qrisp.lanczos.regularize\_S\_H +============================== + +.. currentmodule:: qrisp.lanczos + +.. autofunction:: regularize_S_H \ No newline at end of file diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 7b8d141a6..543fecf5c 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -5,7 +5,7 @@ import scipy def inner_lanczos(H, D, state_prep_func, mes_kwargs): - """ + r""" Perform the quantum subroutine of the exact and efficient Lanczos method to estimate Chebyshev polynomials of a Hamiltonian. This function implements the Krylov space construction via block-encodings @@ -14,8 +14,8 @@ def inner_lanczos(H, D, state_prep_func, mes_kwargs): At each order $k = 0, \dots, 2D-1$ it prepares and measures circuits corresponding either to $\bra{\psi\lfloor k/2\rfloor}R\ket{\psi\lfloor k/2\rfloor}$ for even k, or - $\bra{\psi\lfloor k/2\rfloor U\ket{\psi\lfoor k/2\rfloor}$ for odd k. - The measured statistics encode the expectation values $\angleT_k(H)\rangle$. + $\bra{\psi\lfloor k/2\rfloor}U\ket{\psi\lfloor k/2\rfloor}$ for odd k. + The measured statistics encode the expectation values $\langle T_k(H)\rangle$. Parameters ---------- @@ -82,7 +82,7 @@ def UR(case_indicator, operand, unitaries): return meas_res def compute_expectation(counts): - """ + r""" Convert measurement counts into an expectation value. Assumes measurement outcomes correspond to $\pm 1$ eigenvalues of observables @@ -114,11 +114,13 @@ def compute_expectation(counts): return expval def build_S_H_from_Tk(Tk_expectation, D): - """ - Construct Lanczos overlap matrix S and Hamiltonian matrix H from Chebyshev polynomials. + r""" + Construct the overlap matrix S and the Krylov Hamiltonian matrix H from Chebyshev polynomial expectation values. - Uses Chebyshev recurrence identities to compute matrix elements - (Equations (17), (19) in the reference paper) from measured expectations. + Using Chebyshev recurrence relations, this function generates the matrix elements for + both the overlap matrix (S) and the Hamiltonian matrix (H) in the Krylov subspace. + The approach follows Equations (17) and (19) in "Exact and efficient Lanczos method + on a quantum computer" (arXiv:2208.00567). Parameters ---------- @@ -155,10 +157,14 @@ def Tk(k): return S, H_mat def regularize_S_H(S, H_mat, cutoff=1e-3): - """ - Regularize overlap matrix S by thresholding eigenvalues below cutoff * max_eigenvalue. - Project both S and H onto subspace defined by large eigenvalues. + r""" + Regularize the overlap matrix S by retaining onlz eigenvectors with sufficientlz large eigenvalues and projects H_mat accordinglz. + This function applies a spectral cutoff: only directions in the Krylov subspace with eigenvalues + above ``cutoff * max_eigenvalue`` are kept. Both the overlap matrix (S) and the Hamiltonian matrix (H_mat) + are projected onto this reduced subspace, ensuring numerical stability for subsequent + generalized eigenvalue calculations. + Parameters ---------- S : Overlap matrix @@ -183,7 +189,7 @@ def regularize_S_H(S, H_mat, cutoff=1e-3): return S_reg, H_reg def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): - """ + r""" Exact and efficient Lanczos method on a quantum computer for ground state energy estimation. This function implements the Lanczos method on a quantum computer using block-encodings of Chebyshev From 767ea168ec17be8c585b2b524cc03d36c1bc010f Mon Sep 17 00:00:00 2001 From: mat70593 Date: Fri, 22 Aug 2025 16:13:46 +0200 Subject: [PATCH 11/32] added Lanczos method to the index.rst table in the algorithms folder --- documentation/source/reference/Algorithms/index.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/documentation/source/reference/Algorithms/index.rst b/documentation/source/reference/Algorithms/index.rst index e0059c1c5..e21da9201 100644 --- a/documentation/source/reference/Algorithms/index.rst +++ b/documentation/source/reference/Algorithms/index.rst @@ -27,6 +27,8 @@ This algorithms submodule of Qrisp provides a collection of commonly used quantu - estimating the amount of solutions for a given Grover oracle * - :ref:`Quantum Monte Carlo Integration ` - numerical integration + * - :ref:`Lanczos Method ` + - finding the ground state energy of a Hamiltonian We encourage you to explore these algorithms, delve into their documentation, and experiment with their implementations. From 3bd612b438b77f42f41af7bd3d7f86e358c19d23 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Tue, 26 Aug 2025 14:13:10 +0200 Subject: [PATCH 12/32] Review Lanczos implementation: minor improvements to documentation; fixed imports --- src/qrisp/algorithms/lanczos.py | 146 ++++++++++++++++++++------------ 1 file changed, 90 insertions(+), 56 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 543fecf5c..2f0866fa1 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -1,37 +1,60 @@ -from qrisp.operators import * -from qrisp import * -from qrisp.algorithms.grover import* +""" +******************************************************************************** +* Copyright (c) 2025 the Qrisp authors +* +* This program and the accompanying materials are made available under the +* terms of the Eclipse Public License 2.0 which is available at +* http://www.eclipse.org/legal/epl-2.0. +* +* This Source Code may also be made available under the following Secondary +* Licenses when the conditions for such availability set forth in the Eclipse +* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 +* with the GNU Classpath Exception which is +* available at https://www.gnu.org/software/classpath/license.html. +* +* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 +******************************************************************************** +""" + +#from qrisp.operators import * +from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate +from qrisp.alg_primitives.prepare import prepare +from qrisp.alg_primitives.switch_case import qswitch +from qrisp.algorithms.grover import diffuser +from qrisp.jasp import jrange import numpy as np import scipy + def inner_lanczos(H, D, state_prep_func, mes_kwargs): r""" - Perform the quantum subroutine of the exact and efficient Lanczos method to estimate Chebyshev polynomials of a Hamiltonian. + Perform the quantum subroutine of the exact and efficient Lanczos method to estimate expectation values of Chebyshev polynomials of a Hamiltonian. This function implements the Krylov space construction via block-encodings - of Chebyshev polynomials $T_k(H)$, as described in "Exact and efficient Lanczos method - on a quantum computer" (arXiv:2208.00567). + of Chebyshev polynomials $T_k(H)$, as described in + `"Exact and efficient Lanczos method on a quantum computer" `_. At each order $k = 0, \dots, 2D-1$ it prepares and measures circuits corresponding - either to $\bra{\psi\lfloor k/2\rfloor}R\ket{\psi\lfloor k/2\rfloor}$ for even k, or - $\bra{\psi\lfloor k/2\rfloor}U\ket{\psi\lfloor k/2\rfloor}$ for odd k. + either to $\bra{\psi\lfloor k/2\rfloor}R\ket{\psi\lfloor k/2\rfloor}$ for even $k$, or + $\bra{\psi\lfloor k/2\rfloor}U\ket{\psi\lfloor k/2\rfloor}$ for odd $k$. The measured statistics encode the expectation values $\langle T_k(H)\rangle$. Parameters ---------- - H : Hamiltonian (QubitOperator) - Hamiltonian represented as a Pauli sum operator with block-encoding accessible via `unitaries()`. + H : QubitOperator + Hamiltonian for which to estimate the ground-state energy. D : int Krylov space dimension. Determines maximum Chebyshev order $(2D-1)$. state_prep_func : callable - Function preparing the initial system state (operand) $\ket{\psi_0}$ on a QuantumVariable. - mes_kwargs : dictionary + Function preparing the initial system state $\ket{\psi_0}$ on a (operand) QuantumVariable. + mes_kwargs : dict The keyword arguments for the measurement function. Returns ------- - meas_res : dictionary - Measurement results for each Chebyshev polynomial order, suitable for conversion into expectation values. + meas_res : dict + Measurement results for each Chebyshev polynomial order $k$, suitable for calculation of expectation values $\langle T_k(H)\rangle$. + """ # Extract unitaries and size of the case_indicator QuantumFloat. @@ -81,31 +104,33 @@ def UR(case_indicator, operand, unitaries): return meas_res -def compute_expectation(counts): + +def compute_expectation(meas_res): r""" - Convert measurement counts into an expectation value. + Convert measurement results into an expectation value. Assumes measurement outcomes correspond to $\pm 1$ eigenvalues of observables - (reflection R or block-encoding unitary U). + (reflection $R$ or block-encoding unitary $U$). - For even k we measure the auxilary case_indicator QuantumFloat in the computational basis. We then map: - $+1$ if all-zeroes outcome (projector $2\ket{0}\bra{0} - I), else $-1$ (paper step 2). + For even $k$ we measure the auxilary case_indicator QuantumFloat in the computational basis. We then map: + $+1$ if all-zeros outcome (projector $2\ket{0}\bra{0} - \mathbb{I}$), else $-1$ (paper step 2). - For odd k we perform the hadamard test for the unitary on the controlled unitary and then measure the - auxilary case_indicator QuantumFloat in the Z basis + For odd $k$ we perform the Hadamard test for the block-encoding unitary $U$ and then measure the + auxilary QuantumVariable in the $Z$ basis. Parameters ---------- - counts : dictionary - Measurement counts or probabilities, keyed by outcomes. + meas_res : dict + A dictionary of values and their corresponding measurement probabilities. Returns ------- expval : float Expectation value of the measured observable. + """ expval = 0.0 - for outcome, prob in counts.items(): + for outcome, prob in meas_res.items(): if int(outcome) == 0: expval += prob * 1 @@ -113,30 +138,30 @@ def compute_expectation(counts): expval += prob * (-1) return expval + def build_S_H_from_Tk(Tk_expectation, D): r""" - Construct the overlap matrix S and the Krylov Hamiltonian matrix H from Chebyshev polynomial expectation values. + Construct the overlap matrix $\mathbf{S}$ and the Krylov Hamiltonian matrix $\mathbf{H}$ from Chebyshev polynomial expectation values. Using Chebyshev recurrence relations, this function generates the matrix elements for - both the overlap matrix (S) and the Hamiltonian matrix (H) in the Krylov subspace. - The approach follows Equations (17) and (19) in "Exact and efficient Lanczos method - on a quantum computer" (arXiv:2208.00567). + both the overlap matrix ($\mathbf{S}$) and the Hamiltonian matrix ($\mathbf{H}$) in the Krylov subspace. + The approach follows Equations (17) and (19) in + `"Exact and efficient Lanczos method on a quantum computer" `_. Parameters ---------- - Tk_expectation : dictionary - Mapping k to ⟨T_k(H)⟩. - - D : integer + Tk_expectation : dict + Dictionary of expectations $⟨T_k(H)⟩$ for each Chebyshev polynomial order $k$. + D : int Krylov space dimension. Returns ------- S : ndarray - overlap (Gram) matrix S or Krylov states. - + Overlap (Gram) matrix $\mathbf{S}$ for Krylov states. H_mat : ndarray - Hamiltonian matrix in Krylov subspace. + Hamiltonian matrix $\mathbf{H}$ in Krylov subspace. + """ def Tk(k): k = abs(k) @@ -156,28 +181,32 @@ def Tk(k): ) return S, H_mat + def regularize_S_H(S, H_mat, cutoff=1e-3): r""" - Regularize the overlap matrix S by retaining onlz eigenvectors with sufficientlz large eigenvalues and projects H_mat accordinglz. + Regularize the overlap matrix $\mathbf{S}$ by retaining only eigenvectors with sufficiently large eigenvalues and project the Hamiltonian matrix $\mathbf{H}$ (``H_mat``) accordingly. This function applies a spectral cutoff: only directions in the Krylov subspace with eigenvalues - above ``cutoff * max_eigenvalue`` are kept. Both the overlap matrix (S) and the Hamiltonian matrix (H_mat) + above ``cutoff * max_eigenvalue`` are kept. Both the overlap matrix ($\mathbf{S}$) and the Hamiltonian matrix ($\mathbf{H}$) are projected onto this reduced subspace, ensuring numerical stability for subsequent generalized eigenvalue calculations. Parameters ---------- - S : Overlap matrix - H_mat : Hamiltonian matrix + S : ndarray + Overlap matrix. + H_mat : ndarray + Hamiltonian matrix. cutoff : float - Eigenvalue threshold for regularizing S + Eigenvalue threshold for regularizing $\mathbf{S}$. Returns ------- S_reg : ndarray - regularized overlap matrix S. + Regularized overlap matrix. H_reg : ndarray - regularized Hamiltonian matrix in Krylov subspace. + Regularized Hamiltonian matrix in Krylov subspace. + """ eigvals, eigvecs = np.linalg.eigh(S) max_eigval = eigvals.max() @@ -188,18 +217,19 @@ def regularize_S_H(S, H_mat, cutoff=1e-3): H_reg = U.T @ H_mat @ U return S_reg, H_reg + def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): r""" Exact and efficient Lanczos method on a quantum computer for ground state energy estimation. This function implements the Lanczos method on a quantum computer using block-encodings of Chebyshev polynomials $T_k(H)$, closely following the algorithm proposed in - "Exact and efficient Lanczos method on a quantum computer" (arXiv:2208.00567). + `"Exact and efficient Lanczos method on a quantum computer" `_. The quantum Lanczos algorithm efficiently constructs a Krylov subspace by applying Chebyshev polynomials - of the Hamiltonian to an initial state. The Krylov space dimension D determines the accuracy of ground + of the Hamiltonian to an initial state. 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 + ($|\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, @@ -207,34 +237,34 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): Implements the following steps: 1. Run quantum Lanczos subroutine to obtain Chebyshev expectation values $\langle T_k(H)\rangle$. - 2. Build overlap and Hamiltonian subspace matrices (S, H). - 3. Regularize overlap matrix S and H by projecting onto the subspace with well conditioned eigenvalues. + 2. Build overlap and Hamiltonian subspace matrices $(\mathbf{S}, \mathbf{H})$. + 3. Regularize overlap matrix $\mathbf{S}$ and $\mathbf{H}$ by projecting onto the subspace with well conditioned eigenvalues. 4. Solve generalized eigenvalue problem $\mathbf{H}\vec{v}=\epsilon\mathbf{S}\vec{v}$. 5. Return lowest eigenvalue $\epsilon_{\text{min}}$ as ground state energy estimate. Parameters ---------- - H : Hamiltonian (QuantumOperator) - Hamiltonian operator expressed as the Pauli sum. + H : QubitOperator + Hamiltonian for which to estimate the ground-state energy. D : int Krylov space dimension. state_prep_func : callable - Function preparing the initial system state (operand) $\ket{\psi_0}$ on a QuantumVariable. - mes_kwargs : dictionary + Function preparing the initial system state $\ket{\psi_0}$ on a (operand) QuantumVariable. + mes_kwargs : dict The keyword arguments for the measurement function. cutoff : float - Regularization cutoff threshold for overlap matrix S. + Regularization cutoff threshold for overlap matrix $\mathbf{S}$. Returns ------- energy : float Estimated ground state energy of the Hamiltonian H. - results : dictionary + results : dict Full details including: - - 'Tk_expvals': dict of Chebyshev expectation values - - 'energy': ground state estimate + - 'Tk_expvals': dictionary of Chebyshev expectation values + - 'energy': ground-state energy estimate - 'eigvals': eigenvalues of regularized problem - - 'eigvecs': eigenvectors + - 'eigvecs': eigenvectors of regularized problem - 'S_reg': regularized overlap matrix - 'H_reg': regularized Hamiltonian matrix @@ -243,6 +273,8 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): :: + from qrisp.lanczos import lanczos_alg + from qrisp.operators import X, Y, Z from qrisp.vqe.problems.heisenberg import create_heisenberg_init_function import networkx as nx @@ -253,6 +285,8 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): # Define Hamiltonian e.g. Heisenberg with custom couplings H = (1/4)*sum((X(i)*X(j) + Y(i)*Y(j) + 0.5*Z(i)*Z(j)) for i,j in G.edges()) + print(f"Ground state energy: {H.ground_state_energy()}") + # Prepare initial state function (tensor product of singlets) M = nx.maximal_matching(G) U_singlet = create_heisenberg_init_function(M) From 679fe2a5b3edc735925859f273f8a5b682e26872 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Tue, 26 Aug 2025 14:59:46 +0200 Subject: [PATCH 13/32] Removed jrange import --- src/qrisp/algorithms/lanczos.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 2f0866fa1..cc5331ea3 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -21,7 +21,6 @@ from qrisp.alg_primitives.prepare import prepare from qrisp.alg_primitives.switch_case import qswitch from qrisp.algorithms.grover import diffuser -from qrisp.jasp import jrange import numpy as np import scipy @@ -72,7 +71,7 @@ def UR(case_indicator, operand, unitaries): meas_res = {} - for k in jrange(0, 2*D): + for k in range(0, 2*D): case_indicator = QuantumFloat(n) operand = QuantumVariable(H.find_minimal_qubit_amount()) state_prep_func(operand) # prepare operand $\ket{\psi_0}$ @@ -80,7 +79,7 @@ def UR(case_indicator, operand, unitaries): if k % 2 == 0: # EVEN k: Figure 1 top with conjugate(state_prep)(case_indicator): - for _ in jrange(k//2): + for _ in range(k//2): UR(case_indicator, operand, unitaries) meas = case_indicator.get_measurement(**mes_kwargs) @@ -89,7 +88,7 @@ def UR(case_indicator, operand, unitaries): else: # ODD k: Figure 1 bottom state_prep(case_indicator) - for _ in jrange(k//2): + for _ in range(k//2): UR(case_indicator, operand, unitaries) qv = QuantumVariable(1) h(qv) # Hadamard test for From 607e00230ec02ce997c0f607cc13d8ae8b090965 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Wed, 12 Nov 2025 15:27:53 +0100 Subject: [PATCH 14/32] obtaining block encoding via H.pauli_block_encoding() in inner_lanczos --- src/qrisp/algorithms/lanczos.py | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index cc5331ea3..67d27de24 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -17,10 +17,8 @@ """ #from qrisp.operators import * -from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate -from qrisp.alg_primitives.prepare import prepare +from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate, reflection from qrisp.alg_primitives.switch_case import qswitch -from qrisp.algorithms.grover import diffuser import numpy as np import scipy @@ -57,17 +55,11 @@ def inner_lanczos(H, D, state_prep_func, mes_kwargs): """ # Extract unitaries and size of the case_indicator QuantumFloat. - unitaries, coeffs = H.unitaries() - num_unitaries = len(unitaries) - n = np.int64(np.ceil(np.log2(num_unitaries))) - - # Prepare $\ket{G} = \sum_i \sqrt{\alpha_i}\ket{i}$ on the auxiliary register (case_indicator) (paper Def. 1, Eq. (4)) - def state_prep(case): - prepare(case, np.sqrt(coeffs)) + U, state_prep, n = (H.pauli_block_encoding()) def UR(case_indicator, operand, unitaries): - qswitch(operand, case_indicator, unitaries) #qswitch applies $U = \sum_i\ket{i}\bra{i}\otimes P_i$. - diffuser(case_indicator, state_function=state_prep) # diffuser implements the reflection operator R about $\ket{G}$. + U(case_indicator, operand) # applies $U = \sum_i\ket{i}\bra{i}\otimes P_i$. + reflection(case_indicator, state_function=state_prep) # reflection operator R about $\ket{G}$. meas_res = {} @@ -80,7 +72,7 @@ def UR(case_indicator, operand, unitaries): # EVEN k: Figure 1 top with conjugate(state_prep)(case_indicator): for _ in range(k//2): - UR(case_indicator, operand, unitaries) + UR(case_indicator, operand) meas = case_indicator.get_measurement(**mes_kwargs) meas_res[k] = meas @@ -89,11 +81,11 @@ def UR(case_indicator, operand, unitaries): # ODD k: Figure 1 bottom state_prep(case_indicator) for _ in range(k//2): - UR(case_indicator, operand, unitaries) + UR(case_indicator, operand) qv = QuantumVariable(1) h(qv) # Hadamard test for with control(qv[0]): - qswitch(operand, case_indicator, unitaries) # control-U on the case_indicator QuantumFloat + U(case_indicator, operand) # control-U on the case_indicator QuantumFloat h(qv) # Hadamard test for meas = qv.get_measurement(**mes_kwargs) meas_res[k] = meas From cee3c99c8cca9ce761dd1c1d2834376f750c7811 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Thu, 13 Nov 2025 15:24:35 +0100 Subject: [PATCH 15/32] minor changes, importing qache in qubit_operator --- src/qrisp/algorithms/lanczos.py | 2 +- src/qrisp/operators/qubit/qubit_operator.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 67d27de24..d6f7e0a4f 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -55,7 +55,7 @@ def inner_lanczos(H, D, state_prep_func, mes_kwargs): """ # Extract unitaries and size of the case_indicator QuantumFloat. - U, state_prep, n = (H.pauli_block_encoding()) + U, state_prep, n = H.pauli_block_encoding() def UR(case_indicator, operand, unitaries): U(case_indicator, operand) # applies $U = \sum_i\ket{i}\bra{i}\otimes P_i$. diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py index b082edac6..f91dba236 100644 --- a/src/qrisp/operators/qubit/qubit_operator.py +++ b/src/qrisp/operators/qubit/qubit_operator.py @@ -31,7 +31,7 @@ from qrisp.operators.qubit.commutativity_tools import construct_change_of_basis from qrisp import cx, cz, h, s, sx_dg, IterationEnvironment, conjugate, merge, invert -from qrisp.jasp import check_for_tracing_mode, jrange +from qrisp.jasp import check_for_tracing_mode, jrange, qache threshold = 1e-9 @@ -2054,6 +2054,7 @@ def state_prep(case): return unitaries, np.array(coefficients, dtype=float) + @qache def pauli_block_encoding(self): r""" Returns a block encoding of the operator. @@ -2193,7 +2194,6 @@ def main(): (7.0, 6.0): 0.08333333084980639} """ - from qrisp.jasp import qache from qrisp.alg_primitives import prepare, qswitch unitaries, coeffs = self.unitaries() From f68a2887f373370d2743f23a4a50540636f70191 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Thu, 13 Nov 2025 19:58:06 +0100 Subject: [PATCH 16/32] Restructuring the code: split inner_lanczos into inner_lanczos (preparing circuit) and lanczos_expvals (obtaining the expectation values). Included distinctions for jasp vs non-jasp using check_for_tracing_mode. Tk_expvals are now arrays and not dictionaries anymore. --- src/qrisp/algorithms/lanczos.py | 103 +++++++++++++++++++------------- 1 file changed, 63 insertions(+), 40 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index d6f7e0a4f..4c81835b0 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -16,14 +16,16 @@ ******************************************************************************** """ -#from qrisp.operators import * from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate, reflection from qrisp.alg_primitives.switch_case import qswitch +from qrisp.jasp import jrange, q_cond, check_for_tracing_mode, expectation_value +import jax.numpy as jnp +import jax import numpy as np import scipy -def inner_lanczos(H, D, state_prep_func, mes_kwargs): +def inner_lanczos(H, k, state_prep_func): r""" Perform the quantum subroutine of the exact and efficient Lanczos method to estimate expectation values of Chebyshev polynomials of a Hamiltonian. @@ -57,44 +59,43 @@ def inner_lanczos(H, D, state_prep_func, mes_kwargs): # Extract unitaries and size of the case_indicator QuantumFloat. U, state_prep, n = H.pauli_block_encoding() - def UR(case_indicator, operand, unitaries): + def UR(case_indicator, operand): U(case_indicator, operand) # applies $U = \sum_i\ket{i}\bra{i}\otimes P_i$. reflection(case_indicator, state_function=state_prep) # reflection operator R about $\ket{G}$. - meas_res = {} - - for k in range(0, 2*D): - case_indicator = QuantumFloat(n) - operand = QuantumVariable(H.find_minimal_qubit_amount()) - state_prep_func(operand) # prepare operand $\ket{\psi_0}$ - - if k % 2 == 0: - # EVEN k: Figure 1 top - with conjugate(state_prep)(case_indicator): - for _ in range(k//2): - UR(case_indicator, operand) - - meas = case_indicator.get_measurement(**mes_kwargs) - meas_res[k] = meas + case_indicator = QuantumFloat(n) + operand = QuantumVariable(H.find_minimal_qubit_amount()) + state_prep_func(operand) # prepare operand $\ket{\psi_0}$ - else: - # ODD k: Figure 1 bottom - state_prep(case_indicator) - for _ in range(k//2): + def even(case_indicator, operand, k): + # EVEN k: Figure 1 top + with conjugate(state_prep)(case_indicator): + for _ in jrange(k//2): UR(case_indicator, operand) - qv = QuantumVariable(1) - h(qv) # Hadamard test for - with control(qv[0]): - U(case_indicator, operand) # control-U on the case_indicator QuantumFloat - h(qv) # Hadamard test for - meas = qv.get_measurement(**mes_kwargs) - meas_res[k] = meas + return case_indicator + + def odd(case_indicator, operand, k): + # ODD k: Figure 1 bottom + state_prep(case_indicator) + for _ in jrange(k//2): + UR(case_indicator, operand) + qv = QuantumFloat(1) + h(qv) # Hadamard test for + with control(qv[0]): + U(case_indicator, operand) # control-U on the case_indicator QuantumFloat + h(qv) # Hadamard test for + return qv - case_indicator.delete() - operand.delete() - - return meas_res - + if check_for_tracing_mode(): + x_cond = q_cond + else: + def x_cond(pred, true_fun, false_fun, *operands): + if pred: + return true_fun(*operands) + else: + return false_fun(*operands) + + return x_cond(k%2==0, even, odd, case_indicator, operand, k) def compute_expectation(meas_res): r""" @@ -129,8 +130,33 @@ def compute_expectation(meas_res): expval += prob * (-1) return expval +def lanczos_expvals(H, D, state_prep_func, mes_kwargs={}): + if check_for_tracing_mode(): + + @jax.jit + def post_processor(x): + """ + Returns 1 if the input integer x is 0, and -1 otherwise, using jax.numpy.where. + """ + return jnp.where(x == 0, 1, -1) + + ev_function = expectation_value(inner_lanczos, shots = 50, post_processor = post_processor) + expvals = jnp.zeros(2*D) + + for k in range(0, 2*D): + expval = ev_function(H, k, state_prep_func) + expvals = expvals.at[k].set(expval) + + else: + expvals = np.zeros(2*D) + for k in range(0, 2*D): + qarg = inner_lanczos(H, k, state_prep_func) + meas = qarg.get_measurement(**mes_kwargs) + expvals[k] = compute_expectation(meas) + + return expvals -def build_S_H_from_Tk(Tk_expectation, D): +def build_S_H_from_Tk(Tk_expvals, D): r""" Construct the overlap matrix $\mathbf{S}$ and the Krylov Hamiltonian matrix $\mathbf{H}$ from Chebyshev polynomial expectation values. @@ -156,7 +182,7 @@ def build_S_H_from_Tk(Tk_expectation, D): """ def Tk(k): k = abs(k) - return Tk_expectation.get(k, 0) + return Tk_expvals[k] S = np.zeros((D, D)) H_mat = np.zeros((D, D)) @@ -292,10 +318,7 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): unitaries, coeffs = H.unitaries() # Step 1: Quantum Lanczos: Get expectation values of Chebyshev polynomials - meas_counts = inner_lanczos(H, D, state_prep_func, mes_kwargs) - - # Convert counts to expectation values - Tk_expvals = {k: compute_expectation(counts) for k, counts in meas_counts.items()} + Tk_expvals = lanczos_expvals(H, D, state_prep_func, mes_kwargs) # Step 2: Build matrices S and H S, H_mat = build_S_H_from_Tk(Tk_expvals, D) From 124ab78f3ea72af8c3e4868fdc37a1be16c64145 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Thu, 13 Nov 2025 21:38:37 +0100 Subject: [PATCH 17/32] fixed reflection import --- src/qrisp/algorithms/lanczos.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 4c81835b0..b839727f1 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -16,8 +16,8 @@ ******************************************************************************** """ -from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate, reflection -from qrisp.alg_primitives.switch_case import qswitch +from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate +from qrisp.alg_primitives.switch_case import reflection from qrisp.jasp import jrange, q_cond, check_for_tracing_mode, expectation_value import jax.numpy as jnp import jax @@ -33,7 +33,7 @@ def inner_lanczos(H, k, state_prep_func): of Chebyshev polynomials $T_k(H)$, as described in `"Exact and efficient Lanczos method on a quantum computer" `_. - At each order $k = 0, \dots, 2D-1$ it prepares and measures circuits corresponding + For each polynomial order $k = 0, \dots, 2D-1$, it prepares and measures circuits corresponding either to $\bra{\psi\lfloor k/2\rfloor}R\ket{\psi\lfloor k/2\rfloor}$ for even $k$, or $\bra{\psi\lfloor k/2\rfloor}U\ket{\psi\lfloor k/2\rfloor}$ for odd $k$. The measured statistics encode the expectation values $\langle T_k(H)\rangle$. From c20e3fbb4bf2093bdc613d31949217daf752cc19 Mon Sep 17 00:00:00 2001 From: mat70593 Date: Thu, 13 Nov 2025 21:40:02 +0100 Subject: [PATCH 18/32] fixed reflection import --- src/qrisp/algorithms/lanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index b839727f1..f231852de 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -17,7 +17,7 @@ """ from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate -from qrisp.alg_primitives.switch_case import reflection +from qrisp.alg_primitives import reflection from qrisp.jasp import jrange, q_cond, check_for_tracing_mode, expectation_value import jax.numpy as jnp import jax From 22e4ed86e32aba70c3273b43d1ccac25381282cc Mon Sep 17 00:00:00 2001 From: mat70593 Date: Fri, 14 Nov 2025 15:46:12 +0100 Subject: [PATCH 19/32] fixed the reflection import. for real this time --- src/qrisp/algorithms/lanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index f231852de..2076fd9e6 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -17,7 +17,7 @@ """ from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate -from qrisp.alg_primitives import reflection +from qrisp.alg_primitives.reflection import reflection from qrisp.jasp import jrange, q_cond, check_for_tracing_mode, expectation_value import jax.numpy as jnp import jax From fde14ff8db448edc12c6331c598f5c762373da76 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 14 Nov 2025 18:52:15 +0100 Subject: [PATCH 20/32] Added default amount of shots in Lanczos_expvals function; vectorized and JAX-traceable function for building S and H from Tk --- src/qrisp/algorithms/lanczos.py | 57 +++++++++++++++++++++------------ 1 file changed, 36 insertions(+), 21 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 2076fd9e6..4c8be6da5 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -16,12 +16,13 @@ ******************************************************************************** """ -from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate -from qrisp.alg_primitives.reflection import reflection -from qrisp.jasp import jrange, q_cond, check_for_tracing_mode, expectation_value +from functools import partial import jax.numpy as jnp import jax import numpy as np +from qrisp import QuantumVariable, QuantumFloat, h, control, conjugate +from qrisp.alg_primitives.reflection import reflection +from qrisp.jasp import jrange, q_cond, check_for_tracing_mode, expectation_value import scipy @@ -97,6 +98,7 @@ def x_cond(pred, true_fun, false_fun, *operands): return x_cond(k%2==0, even, odd, case_indicator, operand, k) + def compute_expectation(meas_res): r""" Convert measurement results into an expectation value. @@ -130,7 +132,13 @@ def compute_expectation(meas_res): expval += prob * (-1) return expval + def lanczos_expvals(H, D, state_prep_func, mes_kwargs={}): + + # Set default options + if not "shots" in mes_kwargs: + mes_kwargs["shots"] = 100000 + if check_for_tracing_mode(): @jax.jit @@ -140,7 +148,7 @@ def post_processor(x): """ return jnp.where(x == 0, 1, -1) - ev_function = expectation_value(inner_lanczos, shots = 50, post_processor = post_processor) + ev_function = expectation_value(inner_lanczos, shots = mes_kwargs["shots"], post_processor = post_processor) expvals = jnp.zeros(2*D) for k in range(0, 2*D): @@ -156,6 +164,8 @@ def post_processor(x): return expvals + +@partial(jax.jit, static_argnums=[1]) def build_S_H_from_Tk(Tk_expvals, D): r""" Construct the overlap matrix $\mathbf{S}$ and the Krylov Hamiltonian matrix $\mathbf{H}$ from Chebyshev polynomial expectation values. @@ -167,35 +177,40 @@ def build_S_H_from_Tk(Tk_expvals, D): Parameters ---------- - Tk_expectation : dict + Tk_expectation : numpy.ndarray Dictionary of expectations $⟨T_k(H)⟩$ for each Chebyshev polynomial order $k$. D : int Krylov space dimension. Returns ------- - S : ndarray + S : numpy.ndarray Overlap (Gram) matrix $\mathbf{S}$ for Krylov states. - H_mat : ndarray + H_mat : numpy.ndarray Hamiltonian matrix $\mathbf{H}$ in Krylov subspace. """ - def Tk(k): - k = abs(k) + def Tk_vec(k): + k = jnp.abs(k) return Tk_expvals[k] - S = np.zeros((D, D)) - H_mat = np.zeros((D, D)) - - for i in range(D): - for j in range(D): - S[i, j] = 0.5 * (Tk(i + j) + Tk(abs(i - j))) - H_mat[i, j] = 0.25 * ( - Tk(i + j + 1) - + Tk(abs(i + j - 1)) - + Tk(abs(i - j + 1)) - + Tk(abs(i - j - 1)) - ) + # Create 2D arrays of indices i and j + i_indices = jnp.arange(D, dtype=jnp.int32)[:, None] # Column vector (D, 1) + j_indices = jnp.arange(D, dtype=jnp.int32)[None, :] # Row vector (1, D) + # The combination of these two will broadcast operations across a (D, D) grid + + # Calculate S matrix using vectorized operations + # i+j and abs(i-j) are performed element-wise across the (D, D) grid + S = 0.5 * (Tk_vec(i_indices + j_indices) + Tk_vec(jnp.abs(i_indices - j_indices))) + + # Calculate H_mat matrix using vectorized operations + H_mat = 0.25 * ( + Tk_vec(i_indices + j_indices + 1) + + Tk_vec(jnp.abs(i_indices + j_indices - 1)) + + Tk_vec(jnp.abs(i_indices - j_indices + 1)) + + Tk_vec(jnp.abs(i_indices - j_indices - 1)) + ) + return S, H_mat From 192b5ea71ce17e61b50d5cd034514149d74a5ff5 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 14 Nov 2025 20:43:58 +0100 Subject: [PATCH 21/32] Separate jax postprocessing functions; returning info for lanczos optional; added documentation for lanczos_expvals --- .../source/reference/Algorithms/Lanczos.rst | 6 +- .../source/reference/Algorithms/Shor.rst | 2 +- .../qrisp.lanczos.compute_expectation.rst | 6 - .../qrisp.lanczos.lanczos_expvals.rst | 6 + src/qrisp/algorithms/lanczos.py | 136 +++++++++++++++--- 5 files changed, 129 insertions(+), 27 deletions(-) delete mode 100644 documentation/source/reference/Algorithms/generated/qrisp.lanczos.compute_expectation.rst create mode 100644 documentation/source/reference/Algorithms/generated/qrisp.lanczos.lanczos_expvals.rst diff --git a/documentation/source/reference/Algorithms/Lanczos.rst b/documentation/source/reference/Algorithms/Lanczos.rst index 0b46e7c9a..da39f4675 100644 --- a/documentation/source/reference/Algorithms/Lanczos.rst +++ b/documentation/source/reference/Algorithms/Lanczos.rst @@ -1,8 +1,8 @@ .. _lanczos_alg: -Lanczos method -============== +Lanczos Algorithm +================= .. currentmodule:: qrisp.lanczos @@ -11,7 +11,7 @@ Lanczos method .. autosummary:: :toctree: generated/ + lanczos_expvals inner_lanczos - compute_expectation build_S_H_from_Tk regularize_S_H \ No newline at end of file diff --git a/documentation/source/reference/Algorithms/Shor.rst b/documentation/source/reference/Algorithms/Shor.rst index f4aa8bb8b..377d6fa00 100644 --- a/documentation/source/reference/Algorithms/Shor.rst +++ b/documentation/source/reference/Algorithms/Shor.rst @@ -1,6 +1,6 @@ .. _Shor: -Shor's algorithm +Shor's Algorithm ================ In the realm of quantum computing, where classical limitations are challenged and new horizons are explored, Shor's Algorithm stands as a testament to the transformative potential of quantum mechanics in the field of cryptography. Developed by mathematician Peter Shor in 1994, this groundbreaking algorithm has the power to revolutionize the world of cryptography by efficiently factoring large numbers—once considered an insurmountable task for classical computers. diff --git a/documentation/source/reference/Algorithms/generated/qrisp.lanczos.compute_expectation.rst b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.compute_expectation.rst deleted file mode 100644 index ad4c215a0..000000000 --- a/documentation/source/reference/Algorithms/generated/qrisp.lanczos.compute_expectation.rst +++ /dev/null @@ -1,6 +0,0 @@ -qrisp.lanczos.compute\_expectation -================================== - -.. currentmodule:: qrisp.lanczos - -.. autofunction:: compute_expectation \ No newline at end of file diff --git a/documentation/source/reference/Algorithms/generated/qrisp.lanczos.lanczos_expvals.rst b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.lanczos_expvals.rst new file mode 100644 index 000000000..669434b0b --- /dev/null +++ b/documentation/source/reference/Algorithms/generated/qrisp.lanczos.lanczos_expvals.rst @@ -0,0 +1,6 @@ +qrisp.lanczos.lanczos\_expvals +============================== + +.. currentmodule:: qrisp.lanczos + +.. autofunction:: lanczos_expvals \ No newline at end of file diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 4c8be6da5..3680b07c9 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -34,10 +34,10 @@ def inner_lanczos(H, k, state_prep_func): of Chebyshev polynomials $T_k(H)$, as described in `"Exact and efficient Lanczos method on a quantum computer" `_. - For each polynomial order $k = 0, \dots, 2D-1$, it prepares and measures circuits corresponding + For each polynomial order $k = 0, \dots, 2D-1$, it prepares circuits corresponding either to $\bra{\psi\lfloor k/2\rfloor}R\ket{\psi\lfloor k/2\rfloor}$ for even $k$, or $\bra{\psi\lfloor k/2\rfloor}U\ket{\psi\lfloor k/2\rfloor}$ for odd $k$. - The measured statistics encode the expectation values $\langle T_k(H)\rangle$. + The measured statistics of the prepared quantum states encode the expectation values $\langle T_k(H)\rangle$. Parameters ---------- @@ -47,13 +47,11 @@ def inner_lanczos(H, k, state_prep_func): Krylov space dimension. Determines maximum Chebyshev order $(2D-1)$. state_prep_func : callable Function preparing the initial system state $\ket{\psi_0}$ on a (operand) QuantumVariable. - mes_kwargs : dict - The keyword arguments for the measurement function. Returns ------- - meas_res : dict - Measurement results for each Chebyshev polynomial order $k$, suitable for calculation of expectation values $\langle T_k(H)\rangle$. + QuantumVariable + A QuantumVariable for which the measured statistics encodes the expectation values $\langle T_k(H)\rangle$. """ @@ -134,6 +132,35 @@ def compute_expectation(meas_res): def lanczos_expvals(H, D, state_prep_func, mes_kwargs={}): + r""" + Perform the quantum subroutine of the exact and efficient Lanczos method to estimate expectation values of Chebyshev polynomials of a Hamiltonian. + + This function implements the Krylov space construction via block-encodings + of Chebyshev polynomials $T_k(H)$, as described in + `"Exact and efficient Lanczos method on a quantum computer" `_. + + For each polynomial order $k = 0, \dotsc, 2D-1$, it prepares and measures circuits corresponding + either to $\bra{\psi\lfloor k/2\rfloor}R\ket{\psi\lfloor k/2\rfloor}$ for even $k$, or + $\bra{\psi\lfloor k/2\rfloor}U\ket{\psi\lfloor k/2\rfloor}$ for odd $k$. + The measured statistics encode the expectation values $\langle T_k(H)\rangle$. + + Parameters + ---------- + H : QubitOperator + Hamiltonian for which to estimate the ground-state energy. + D : int + Krylov space dimension. Determines maximum Chebyshev order $(2D-1)$. + state_prep_func : callable + Function preparing the initial system state $\ket{\psi_0}$ on a (operand) QuantumVariable. + mes_kwargs : dict, optional + The keyword arguments for the measurement function. + + Returns + ------- + expvals : ndarray + The expectation values $\langle T_k(H)\rangle$ for $k=0, \dotsc, 2D-1$. + + """ # Set default options if not "shots" in mes_kwargs: @@ -165,7 +192,6 @@ def post_processor(x): return expvals -@partial(jax.jit, static_argnums=[1]) def build_S_H_from_Tk(Tk_expvals, D): r""" Construct the overlap matrix $\mathbf{S}$ and the Krylov Hamiltonian matrix $\mathbf{H}$ from Chebyshev polynomial expectation values. @@ -177,19 +203,46 @@ def build_S_H_from_Tk(Tk_expvals, D): Parameters ---------- - Tk_expectation : numpy.ndarray + Tk_expectation : ndarray Dictionary of expectations $⟨T_k(H)⟩$ for each Chebyshev polynomial order $k$. D : int Krylov space dimension. Returns ------- - S : numpy.ndarray + S : ndarray Overlap (Gram) matrix $\mathbf{S}$ for Krylov states. - H_mat : numpy.ndarray + H_mat : ndarray Hamiltonian matrix $\mathbf{H}$ in Krylov subspace. """ + def Tk_vec(k): + k = np.abs(k) + return Tk_expvals[k] + + # Create 2D arrays of indices i and j + i_indices = np.arange(D, dtype=np.int32)[:, None] # Column vector (D, 1) + j_indices = np.arange(D, dtype=np.int32)[None, :] # Row vector (1, D) + # The combination of these two will broadcast operations across a (D, D) grid + + # Calculate S matrix using vectorized operations + # i+j and abs(i-j) are performed element-wise across the (D, D) grid + S = 0.5 * (Tk_vec(i_indices + j_indices) + Tk_vec(np.abs(i_indices - j_indices))) + + # Calculate H_mat matrix using vectorized operations + H_mat = 0.25 * ( + Tk_vec(i_indices + j_indices + 1) + + Tk_vec(np.abs(i_indices + j_indices - 1)) + + Tk_vec(np.abs(i_indices - j_indices + 1)) + + Tk_vec(np.abs(i_indices - j_indices - 1)) + ) + + return S, H_mat + + +@partial(jax.jit, static_argnums=(1,)) +def build_S_H_from_Tk_jax(Tk_expvals, D): + def Tk_vec(k): k = jnp.abs(k) return Tk_expvals[k] @@ -250,7 +303,30 @@ def regularize_S_H(S, H_mat, cutoff=1e-3): return S_reg, H_reg -def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): +@partial(jax.jit, static_argnums=(2,)) # max_D_out must be static +def regularize_S_H_jax(S, H_mat, max_D_out, cutoff=1e-3): + eigvals, eigvecs = jnp.linalg.eigh(S) + + max_eigval = eigvals.max() + threshold = cutoff * max_eigval + + mask = eigvals > threshold + + # Use nonzero with a fixed size and fill_value to handle static shapes + kept_indices = jnp.nonzero(mask, size=max_D_out, fill_value=0)[0] + + # Use jnp.take to select eigenvectors. The resulting U has a static shape. + # We might need to ensure the indices array has the correct shape for broadcasting + U = jnp.take(eigvecs, kept_indices, axis=1) + + S_reg = U.T @ S @ U + H_reg = U.T @ H_mat @ U + + # Note: S_reg and H_reg will be of shape (max_D_out, max_D_out) + return S_reg, H_reg + + +def lanczos_alg(H, D, state_prep_func, mes_kwargs={}, cutoff=1e-2, show_info=False): r""" Exact and efficient Lanczos method on a quantum computer for ground state energy estimation. @@ -285,13 +361,15 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): mes_kwargs : dict The keyword arguments for the measurement function. cutoff : float - Regularization cutoff threshold for overlap matrix $\mathbf{S}$. + Regularization cutoff threshold for overlap matrix $\mathbf{S}$. The default is 1e-2. + show_info : bool + If True, a dictionary with detailed information is returned. The default is False. Returns ------- energy : float Estimated ground state energy of the Hamiltonian H. - results : dict + info : dict, optional Full details including: - 'Tk_expvals': dictionary of Chebyshev expectation values - 'energy': ground-state energy estimate @@ -330,22 +408,42 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): """ + + if check_for_tracing_mode(): + raise NotImplementedError("Solving the generalized eigenvalue problem not implemented in JAX 0.6") + unitaries, coeffs = H.unitaries() # Step 1: Quantum Lanczos: Get expectation values of Chebyshev polynomials Tk_expvals = lanczos_expvals(H, D, state_prep_func, mes_kwargs) + #if check_for_tracing_mode(): + + # Step 2: Build matrices S and H + #S, H_mat = build_S_H_from_Tk_jax(Tk_expvals, D) + + # Step 3: Regularize matrices via thresholding + #S_reg, H_reg = regularize_S_H_jax(S, H_mat, D, cutoff=cutoff) + + # Step 4: Solve generalized eigenvalue problem $\mathbf{H}\vec{v}=\epsilon\mathbf{S}\vec{v}$ + #evals, evecs = solve_generalized_eigenproblem_jax(H_reg, S_reg) + #evals, evecs = jax.scipy.linalg.eigh(H_reg, S_reg) # Solving the generalized eigenvalue problem not implemented in JAX 0.6 + + #ground_state_energy = jnp.min(evals) * jnp.sum(coeffs) + + # Step 2: Build matrices S and H S, H_mat = build_S_H_from_Tk(Tk_expvals, D) # Step 3: Regularize matrices via thresholding - S_reg, H_reg = regularize_S_H(S, H_mat, cutoff=cutoff) + S_reg, H_reg = regularize_S_H(S, H_mat, cutoff=cutoff) # Step 4: Solve generalized eigenvalue problem $\mathbf{H}\vec{v}=\epsilon\mathbf{S}\vec{v}$ - evals, evecs = scipy.linalg.eigh(H_reg, S_reg) + evals, evecs = scipy.linalg.eigh(H_reg, S_reg) - ground_state_energy = np.min(evals)*sum(coeffs) + ground_state_energy = np.min(evals) * np.sum(coeffs) + results = { 'Tk_expvals': Tk_expvals, 'energy': ground_state_energy, @@ -354,4 +452,8 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs = {}, cutoff=1e-2): 'S_reg': S_reg, 'H_reg': H_reg, } - return ground_state_energy, results \ No newline at end of file + + if show_info: + return ground_state_energy, results + else: + return ground_state_energy \ No newline at end of file From 2586d9389ec952d0661f9a32925c6dae2d4cacb8 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 14 Nov 2025 20:50:27 +0100 Subject: [PATCH 22/32] Fixed minor issue in lanczos example --- src/qrisp/algorithms/lanczos.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 3680b07c9..f4fc4f1f0 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -402,7 +402,7 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs={}, cutoff=1e-2, show_info=Fal U_singlet = create_heisenberg_init_function(M) D = 6 # Krylov dimension - energy, info = lanczos_alg(H, D, U_singlet) + energy, info = lanczos_alg(H, D, U_singlet, show_info=True) print(f"Ground state energy estimate: {energy}") From 1f172e2fca419744cb62b16d6f7be0f7aae26108 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 14 Nov 2025 21:42:18 +0100 Subject: [PATCH 23/32] Changed interface to lanczos algorithm: now receives and init_state function that returns the operand --- src/qrisp/algorithms/lanczos.py | 38 ++++++++++++++++++++------------- 1 file changed, 23 insertions(+), 15 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index f4fc4f1f0..d3ac365cf 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -26,7 +26,7 @@ import scipy -def inner_lanczos(H, k, state_prep_func): +def inner_lanczos(H, k, init_state): r""" Perform the quantum subroutine of the exact and efficient Lanczos method to estimate expectation values of Chebyshev polynomials of a Hamiltonian. @@ -45,8 +45,9 @@ def inner_lanczos(H, k, state_prep_func): Hamiltonian for which to estimate the ground-state energy. D : int Krylov space dimension. Determines maximum Chebyshev order $(2D-1)$. - state_prep_func : callable - Function preparing the initial system state $\ket{\psi_0}$ on a (operand) QuantumVariable. + init_state : callable + Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., + ``operand=init_state()``. Returns ------- @@ -63,8 +64,7 @@ def UR(case_indicator, operand): reflection(case_indicator, state_function=state_prep) # reflection operator R about $\ket{G}$. case_indicator = QuantumFloat(n) - operand = QuantumVariable(H.find_minimal_qubit_amount()) - state_prep_func(operand) # prepare operand $\ket{\psi_0}$ + operand = init_state() def even(case_indicator, operand, k): # EVEN k: Figure 1 top @@ -131,7 +131,7 @@ def compute_expectation(meas_res): return expval -def lanczos_expvals(H, D, state_prep_func, mes_kwargs={}): +def lanczos_expvals(H, D, init_state, mes_kwargs={}): r""" Perform the quantum subroutine of the exact and efficient Lanczos method to estimate expectation values of Chebyshev polynomials of a Hamiltonian. @@ -150,8 +150,9 @@ def lanczos_expvals(H, D, state_prep_func, mes_kwargs={}): Hamiltonian for which to estimate the ground-state energy. D : int Krylov space dimension. Determines maximum Chebyshev order $(2D-1)$. - state_prep_func : callable - Function preparing the initial system state $\ket{\psi_0}$ on a (operand) QuantumVariable. + init_state : callable + Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., + ``operand=init_state()``. mes_kwargs : dict, optional The keyword arguments for the measurement function. @@ -179,13 +180,13 @@ def post_processor(x): expvals = jnp.zeros(2*D) for k in range(0, 2*D): - expval = ev_function(H, k, state_prep_func) + expval = ev_function(H, k, init_state) expvals = expvals.at[k].set(expval) else: expvals = np.zeros(2*D) for k in range(0, 2*D): - qarg = inner_lanczos(H, k, state_prep_func) + qarg = inner_lanczos(H, k, init_state) meas = qarg.get_measurement(**mes_kwargs) expvals[k] = compute_expectation(meas) @@ -326,7 +327,7 @@ def regularize_S_H_jax(S, H_mat, max_D_out, cutoff=1e-3): return S_reg, H_reg -def lanczos_alg(H, D, state_prep_func, mes_kwargs={}, cutoff=1e-2, show_info=False): +def lanczos_alg(H, D, init_state, mes_kwargs={}, cutoff=1e-2, show_info=False): r""" Exact and efficient Lanczos method on a quantum computer for ground state energy estimation. @@ -356,8 +357,9 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs={}, cutoff=1e-2, show_info=Fal Hamiltonian for which to estimate the ground-state energy. D : int Krylov space dimension. - state_prep_func : callable - Function preparing the initial system state $\ket{\psi_0}$ on a (operand) QuantumVariable. + init_state : callable + Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., + ``operand=init_state()``. mes_kwargs : dict The keyword arguments for the measurement function. cutoff : float @@ -383,6 +385,7 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs={}, cutoff=1e-2, show_info=Fal :: + from qrisp import QuantumVariable from qrisp.lanczos import lanczos_alg from qrisp.operators import X, Y, Z from qrisp.vqe.problems.heisenberg import create_heisenberg_init_function @@ -401,8 +404,13 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs={}, cutoff=1e-2, show_info=Fal M = nx.maximal_matching(G) U_singlet = create_heisenberg_init_function(M) + def init_state(): + qv = QuantumVariable(H.find_minimal_qubit_amount()) + U_singlet(qv) + return qv + D = 6 # Krylov dimension - energy, info = lanczos_alg(H, D, U_singlet, show_info=True) + energy, info = lanczos_alg(H, D, init_state, show_info=True) print(f"Ground state energy estimate: {energy}") @@ -415,7 +423,7 @@ def lanczos_alg(H, D, state_prep_func, mes_kwargs={}, cutoff=1e-2, show_info=Fal unitaries, coeffs = H.unitaries() # Step 1: Quantum Lanczos: Get expectation values of Chebyshev polynomials - Tk_expvals = lanczos_expvals(H, D, state_prep_func, mes_kwargs) + Tk_expvals = lanczos_expvals(H, D, init_state, mes_kwargs) #if check_for_tracing_mode(): From a7f1b4cc81b3335c9b26b192045ecb15b42ef2cc Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 14 Nov 2025 21:49:21 +0100 Subject: [PATCH 24/32] Added test case for lanczos algorithm --- .../test_lanczos_algorithm.py | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 tests/algorithms_tests/test_lanczos_algorithm.py diff --git a/tests/algorithms_tests/test_lanczos_algorithm.py b/tests/algorithms_tests/test_lanczos_algorithm.py new file mode 100644 index 000000000..546fd2ac4 --- /dev/null +++ b/tests/algorithms_tests/test_lanczos_algorithm.py @@ -0,0 +1,52 @@ +""" +******************************************************************************** +* Copyright (c) 2025 the Qrisp authors +* +* This program and the accompanying materials are made available under the +* terms of the Eclipse Public License 2.0 which is available at +* http://www.eclipse.org/legal/epl-2.0. +* +* This Source Code may also be made available under the following Secondary +* Licenses when the conditions for such availability set forth in the Eclipse +* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 +* with the GNU Classpath Exception which is +* available at https://www.gnu.org/software/classpath/license.html. +* +* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 +******************************************************************************** +""" + +import networkx as nx +import numpy as np +from qrisp import QuantumVariable +from qrisp.lanczos import lanczos_alg +from qrisp.operators import X, Y, Z +from qrisp.vqe.problems.heisenberg import create_heisenberg_init_function + + +def test_lanczos_algorithm(): + + L = 6 + G = nx.Graph() + G.add_edges_from([(k, (k+1) % L) for k in range(L - 1)]) + + # Define Hamiltonian e.g. Heisenberg with custom couplings + H = (1/4)*sum((X(i)*X(j) + Y(i)*Y(j) + 0.5*Z(i)*Z(j)) for i,j in G.edges()) + + energy_classical = H.ground_state_energy() + print(f"Ground state energy: {energy_classical}") + + # Prepare initial state function (tensor product of singlets) + M = nx.maximal_matching(G) + U_singlet = create_heisenberg_init_function(M) + + def init_state(): + qv = QuantumVariable(H.find_minimal_qubit_amount()) + U_singlet(qv) + return qv + + D = 6 # Krylov dimension + energy, info = lanczos_alg(H, D, init_state, show_info=True) + print(f"Ground state energy estimate: {energy}") + + assert np.abs(energy_classical-energy) < 1e-2 \ No newline at end of file From 9027652401d37fce339a71f1a239c2506daaf02c Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 14 Nov 2025 21:50:41 +0100 Subject: [PATCH 25/32] Minor change in test for lanczos --- tests/algorithms_tests/test_lanczos_algorithm.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/algorithms_tests/test_lanczos_algorithm.py b/tests/algorithms_tests/test_lanczos_algorithm.py index 546fd2ac4..ca788cbe6 100644 --- a/tests/algorithms_tests/test_lanczos_algorithm.py +++ b/tests/algorithms_tests/test_lanczos_algorithm.py @@ -33,8 +33,8 @@ def test_lanczos_algorithm(): # Define Hamiltonian e.g. Heisenberg with custom couplings H = (1/4)*sum((X(i)*X(j) + Y(i)*Y(j) + 0.5*Z(i)*Z(j)) for i,j in G.edges()) - energy_classical = H.ground_state_energy() - print(f"Ground state energy: {energy_classical}") + energy_exact = H.ground_state_energy() + print(f"Ground state energy: {energy_exact}") # Prepare initial state function (tensor product of singlets) M = nx.maximal_matching(G) @@ -49,4 +49,4 @@ def init_state(): energy, info = lanczos_alg(H, D, init_state, show_info=True) print(f"Ground state energy estimate: {energy}") - assert np.abs(energy_classical-energy) < 1e-2 \ No newline at end of file + assert np.abs(energy_exact-energy) < 1e-2 \ No newline at end of file From 20efe561720bdb4b101ee08bc8f47a807dca0326 Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 14 Nov 2025 21:54:25 +0100 Subject: [PATCH 26/32] Completed Jasp version of lanczos; currently deactivated since solving generalized eigenvalue problem is not implemented --- src/qrisp/algorithms/lanczos.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index d3ac365cf..edf146a14 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -425,20 +425,22 @@ def init_state(): # Step 1: Quantum Lanczos: Get expectation values of Chebyshev polynomials Tk_expvals = lanczos_expvals(H, D, init_state, mes_kwargs) - #if check_for_tracing_mode(): + """ + if check_for_tracing_mode(): # Step 2: Build matrices S and H - #S, H_mat = build_S_H_from_Tk_jax(Tk_expvals, D) + S, H_mat = build_S_H_from_Tk_jax(Tk_expvals, D) # Step 3: Regularize matrices via thresholding - #S_reg, H_reg = regularize_S_H_jax(S, H_mat, D, cutoff=cutoff) + S_reg, H_reg = regularize_S_H_jax(S, H_mat, D, cutoff=cutoff) # Step 4: Solve generalized eigenvalue problem $\mathbf{H}\vec{v}=\epsilon\mathbf{S}\vec{v}$ - #evals, evecs = solve_generalized_eigenproblem_jax(H_reg, S_reg) - #evals, evecs = jax.scipy.linalg.eigh(H_reg, S_reg) # Solving the generalized eigenvalue problem not implemented in JAX 0.6 + evals, evecs = jax.scipy.linalg.eigh(H_reg, S_reg) # Solving the generalized eigenvalue problem not implemented in JAX 0.6 - #ground_state_energy = jnp.min(evals) * jnp.sum(coeffs) + ground_state_energy = jnp.min(evals) * jnp.sum(coeffs) + return ground_state_energy + """ # Step 2: Build matrices S and H S, H_mat = build_S_H_from_Tk(Tk_expvals, D) From 6fb36f77db16da5e9cbb042b65d75c50c0173aea Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 14 Nov 2025 23:20:21 +0100 Subject: [PATCH 27/32] Some improvements for lanczos documentation --- src/qrisp/algorithms/lanczos.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index edf146a14..5c76ee996 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -155,6 +155,7 @@ def lanczos_expvals(H, D, init_state, mes_kwargs={}): ``operand=init_state()``. mes_kwargs : dict, optional The keyword arguments for the measurement function. + By default, 100_000 ``shots`` are executed for measuring each expectation value. Returns ------- @@ -361,10 +362,11 @@ def lanczos_alg(H, D, init_state, mes_kwargs={}, cutoff=1e-2, show_info=False): Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., ``operand=init_state()``. mes_kwargs : dict - The keyword arguments for the measurement function. + The keyword arguments for the measurement function. + By default, 100_000 ``shots`` are executed for measuring each expectation value. cutoff : float Regularization cutoff threshold for overlap matrix $\mathbf{S}$. The default is 1e-2. - show_info : bool + show_info : bool, optional If True, a dictionary with detailed information is returned. The default is False. Returns @@ -373,12 +375,18 @@ def lanczos_alg(H, D, init_state, mes_kwargs={}, cutoff=1e-2, show_info=False): Estimated ground state energy of the Hamiltonian H. info : dict, optional Full details including: - - 'Tk_expvals': dictionary of Chebyshev expectation values - - 'energy': ground-state energy estimate - - 'eigvals': eigenvalues of regularized problem - - 'eigvecs': eigenvectors of regularized problem - - 'S_reg': regularized overlap matrix - - 'H_reg': regularized Hamiltonian matrix + - 'Tk_expvals' : ndarray + Chebyshev expectation values + - 'energy' : float + Ground-state energy estimate + - 'eigvals' : ndarray + Eigenvalues of regularized problem + - 'eigvecs' : ndarray + Eigenvectors of regularized problem + - 'S_reg' : ndarray + Regularized overlap matrix + - 'H_reg' : ndarray + Regularized Hamiltonian matrix Examples -------- From 4f0b72310ef1b84e1f3713b312b1a9b95d544bed Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Sat, 15 Nov 2025 11:56:55 +0100 Subject: [PATCH 28/32] Minor fix for VQE test --- tests/operators_tests/test_VQE_electronic_structure.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/operators_tests/test_VQE_electronic_structure.py b/tests/operators_tests/test_VQE_electronic_structure.py index d237c2880..ad066c744 100644 --- a/tests/operators_tests/test_VQE_electronic_structure.py +++ b/tests/operators_tests/test_VQE_electronic_structure.py @@ -38,7 +38,7 @@ def test_vqe_electronic_structure_H2(): basis = 'sto-3g') H = create_electronic_hamiltonian(mol).to_qubit_operator() - assert np.abs(H.ground_state_energy()-(-1.85238817356958)) + assert np.abs(H.ground_state_energy()-(-1.85238817356958)) < 1e-9 vqe = electronic_structure_problem(mol) @@ -90,7 +90,7 @@ def test_vqe_electronic_structure_BeH2(): basis = 'sto-3g') H = create_electronic_hamiltonian(mol,active_orb=6,active_elec=4).to_qubit_operator() - assert np.abs(H.ground_state_energy()-(-16.73195995959339)) + assert np.abs(H.ground_state_energy()-(-16.73195995959339)) < 1e-9 # runs for >1 minute vqe = electronic_structure_problem(mol,active_orb=6,active_elec=4) From 8c79a3f7ae043554201fd5da012dbc0725e6e19f Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 21 Nov 2025 17:47:25 +0100 Subject: [PATCH 29/32] Rename init_state to operand_prep in Lanczos algorithm --- .../source/reference/Algorithms/index.rst | 2 +- src/qrisp/algorithms/lanczos.py | 30 +++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/documentation/source/reference/Algorithms/index.rst b/documentation/source/reference/Algorithms/index.rst index 5438bf075..0e32c4593 100644 --- a/documentation/source/reference/Algorithms/index.rst +++ b/documentation/source/reference/Algorithms/index.rst @@ -27,7 +27,7 @@ This algorithms submodule of Qrisp provides a collection of commonly used quantu - estimating the amount of solutions for a given Grover oracle * - :ref:`Quantum Monte Carlo Integration ` - numerical integration - * - :ref:`Lanczos Method ` + * - :ref:`Lanczos Algorithm ` - finding the ground state energy of a Hamiltonian * - :ref:`CKS Algorithm ` - solving the quantum linear system problem diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 5c76ee996..bbc8d6e9b 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -26,7 +26,7 @@ import scipy -def inner_lanczos(H, k, init_state): +def inner_lanczos(H, k, operand_prep): r""" Perform the quantum subroutine of the exact and efficient Lanczos method to estimate expectation values of Chebyshev polynomials of a Hamiltonian. @@ -45,9 +45,9 @@ def inner_lanczos(H, k, init_state): Hamiltonian for which to estimate the ground-state energy. D : int Krylov space dimension. Determines maximum Chebyshev order $(2D-1)$. - init_state : callable + operand_prep : callable Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., - ``operand=init_state()``. + ``operand=operand_prep()``. Returns ------- @@ -64,7 +64,7 @@ def UR(case_indicator, operand): reflection(case_indicator, state_function=state_prep) # reflection operator R about $\ket{G}$. case_indicator = QuantumFloat(n) - operand = init_state() + operand = operand_prep() def even(case_indicator, operand, k): # EVEN k: Figure 1 top @@ -131,7 +131,7 @@ def compute_expectation(meas_res): return expval -def lanczos_expvals(H, D, init_state, mes_kwargs={}): +def lanczos_expvals(H, D, operand_prep, mes_kwargs={}): r""" Perform the quantum subroutine of the exact and efficient Lanczos method to estimate expectation values of Chebyshev polynomials of a Hamiltonian. @@ -150,9 +150,9 @@ def lanczos_expvals(H, D, init_state, mes_kwargs={}): Hamiltonian for which to estimate the ground-state energy. D : int Krylov space dimension. Determines maximum Chebyshev order $(2D-1)$. - init_state : callable + operand_prep : callable Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., - ``operand=init_state()``. + ``operand=operand_prep()``. mes_kwargs : dict, optional The keyword arguments for the measurement function. By default, 100_000 ``shots`` are executed for measuring each expectation value. @@ -181,13 +181,13 @@ def post_processor(x): expvals = jnp.zeros(2*D) for k in range(0, 2*D): - expval = ev_function(H, k, init_state) + expval = ev_function(H, k, operand_prep) expvals = expvals.at[k].set(expval) else: expvals = np.zeros(2*D) for k in range(0, 2*D): - qarg = inner_lanczos(H, k, init_state) + qarg = inner_lanczos(H, k, operand_prep) meas = qarg.get_measurement(**mes_kwargs) expvals[k] = compute_expectation(meas) @@ -328,7 +328,7 @@ def regularize_S_H_jax(S, H_mat, max_D_out, cutoff=1e-3): return S_reg, H_reg -def lanczos_alg(H, D, init_state, mes_kwargs={}, cutoff=1e-2, show_info=False): +def lanczos_alg(H, D, operand_prep, mes_kwargs={}, cutoff=1e-2, show_info=False): r""" Exact and efficient Lanczos method on a quantum computer for ground state energy estimation. @@ -358,9 +358,9 @@ def lanczos_alg(H, D, init_state, mes_kwargs={}, cutoff=1e-2, show_info=False): Hamiltonian for which to estimate the ground-state energy. D : int Krylov space dimension. - init_state : callable + operand_prep : callable Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., - ``operand=init_state()``. + ``operand=operand_prep()``. mes_kwargs : dict The keyword arguments for the measurement function. By default, 100_000 ``shots`` are executed for measuring each expectation value. @@ -412,13 +412,13 @@ def lanczos_alg(H, D, init_state, mes_kwargs={}, cutoff=1e-2, show_info=False): M = nx.maximal_matching(G) U_singlet = create_heisenberg_init_function(M) - def init_state(): + def operand_prep(): qv = QuantumVariable(H.find_minimal_qubit_amount()) U_singlet(qv) return qv D = 6 # Krylov dimension - energy, info = lanczos_alg(H, D, init_state, show_info=True) + energy, info = lanczos_alg(H, D, operand_prep, show_info=True) print(f"Ground state energy estimate: {energy}") @@ -431,7 +431,7 @@ def init_state(): unitaries, coeffs = H.unitaries() # Step 1: Quantum Lanczos: Get expectation values of Chebyshev polynomials - Tk_expvals = lanczos_expvals(H, D, init_state, mes_kwargs) + Tk_expvals = lanczos_expvals(H, D, operand_prep, mes_kwargs) """ if check_for_tracing_mode(): From 08fabdcf9ab159c9c0d4d80d4ae37b982cbfce3a Mon Sep 17 00:00:00 2001 From: Rene Zander Date: Fri, 21 Nov 2025 17:50:26 +0100 Subject: [PATCH 30/32] Remove workaround for q_cond --- src/qrisp/algorithms/lanczos.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index bbc8d6e9b..9a59dd19e 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -84,17 +84,8 @@ def odd(case_indicator, operand, k): U(case_indicator, operand) # control-U on the case_indicator QuantumFloat h(qv) # Hadamard test for return qv - - if check_for_tracing_mode(): - x_cond = q_cond - else: - def x_cond(pred, true_fun, false_fun, *operands): - if pred: - return true_fun(*operands) - else: - return false_fun(*operands) - return x_cond(k%2==0, even, odd, case_indicator, operand, k) + return q_cond(k%2==0, even, odd, case_indicator, operand, k) def compute_expectation(meas_res): From a3024d1341fced745477c86d83e58e190f2f6b55 Mon Sep 17 00:00:00 2001 From: renezander90 Date: Sun, 30 Nov 2025 23:49:04 +0000 Subject: [PATCH 31/32] JAX implementation for solving generalized eigenvalue problem; improved JAX implementation for regularization of H, S matrices; Lanczos algorithm is now Jasp compatible --- src/qrisp/algorithms/lanczos.py | 128 +++++++++++++++++++++----------- 1 file changed, 85 insertions(+), 43 deletions(-) diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py index 9a59dd19e..3e242d454 100644 --- a/src/qrisp/algorithms/lanczos.py +++ b/src/qrisp/algorithms/lanczos.py @@ -184,6 +184,7 @@ def post_processor(x): return expvals +# Postprocessing def build_S_H_from_Tk(Tk_expvals, D): r""" @@ -233,34 +234,7 @@ def Tk_vec(k): return S, H_mat -@partial(jax.jit, static_argnums=(1,)) -def build_S_H_from_Tk_jax(Tk_expvals, D): - - def Tk_vec(k): - k = jnp.abs(k) - return Tk_expvals[k] - - # Create 2D arrays of indices i and j - i_indices = jnp.arange(D, dtype=jnp.int32)[:, None] # Column vector (D, 1) - j_indices = jnp.arange(D, dtype=jnp.int32)[None, :] # Row vector (1, D) - # The combination of these two will broadcast operations across a (D, D) grid - - # Calculate S matrix using vectorized operations - # i+j and abs(i-j) are performed element-wise across the (D, D) grid - S = 0.5 * (Tk_vec(i_indices + j_indices) + Tk_vec(jnp.abs(i_indices - j_indices))) - - # Calculate H_mat matrix using vectorized operations - H_mat = 0.25 * ( - Tk_vec(i_indices + j_indices + 1) - + Tk_vec(jnp.abs(i_indices + j_indices - 1)) - + Tk_vec(jnp.abs(i_indices - j_indices + 1)) - + Tk_vec(jnp.abs(i_indices - j_indices - 1)) - ) - - return S, H_mat - - -def regularize_S_H(S, H_mat, cutoff=1e-3): +def regularize_S_H(S, H_mat, cutoff=1e-2): r""" Regularize the overlap matrix $\mathbf{S}$ by retaining only eigenvectors with sufficiently large eigenvalues and project the Hamiltonian matrix $\mathbf{H}$ (``H_mat``) accordingly. @@ -295,9 +269,40 @@ def regularize_S_H(S, H_mat, cutoff=1e-3): H_reg = U.T @ H_mat @ U return S_reg, H_reg +# JAX postprocessing + +@partial(jax.jit, static_argnums=(1,)) +def build_S_H_from_Tk_jitted(Tk_expvals, D): + + def Tk_vec(k): + k = jnp.abs(k) + return Tk_expvals[k] + + # Create 2D arrays of indices i and j + i_indices = jnp.arange(D, dtype=jnp.int32)[:, None] # Column vector (D, 1) + j_indices = jnp.arange(D, dtype=jnp.int32)[None, :] # Row vector (1, D) + # The combination of these two will broadcast operations across a (D, D) grid + + # Calculate S matrix using vectorized operations + # i+j and abs(i-j) are performed element-wise across the (D, D) grid + S = 0.5 * (Tk_vec(i_indices + j_indices) + Tk_vec(jnp.abs(i_indices - j_indices))) + + # Calculate H_mat matrix using vectorized operations + H_mat = 0.25 * ( + Tk_vec(i_indices + j_indices + 1) + + Tk_vec(jnp.abs(i_indices + j_indices - 1)) + + Tk_vec(jnp.abs(i_indices - j_indices + 1)) + + Tk_vec(jnp.abs(i_indices - j_indices - 1)) + ) + + return S, H_mat + + +@jax.jit +def regularize_S_H_jitted(S, H_mat, cutoff=1e-2): + + D = S.shape[0] -@partial(jax.jit, static_argnums=(2,)) # max_D_out must be static -def regularize_S_H_jax(S, H_mat, max_D_out, cutoff=1e-3): eigvals, eigvecs = jnp.linalg.eigh(S) max_eigval = eigvals.max() @@ -306,17 +311,56 @@ def regularize_S_H_jax(S, H_mat, max_D_out, cutoff=1e-3): mask = eigvals > threshold # Use nonzero with a fixed size and fill_value to handle static shapes - kept_indices = jnp.nonzero(mask, size=max_D_out, fill_value=0)[0] + kept_indices = jnp.nonzero(mask, size=D, fill_value=0)[0] # Use jnp.take to select eigenvectors. The resulting U has a static shape. # We might need to ensure the indices array has the correct shape for broadcasting U = jnp.take(eigvecs, kept_indices, axis=1) - S_reg = U.T @ S @ U - H_reg = U.T @ H_mat @ U + # Replace the vectors for the remaining indices by all zeros vectors + new_mask = jnp.sort(mask)[::-1] + zeros = jnp.zeros(D) + V = jnp.zeros((D, D), dtype=U.dtype) + V = V.at[jnp.arange(D), :].set(jnp.where(new_mask, U, zeros)) + + S_reg_ext = V.T @ S @ V + H_reg_ext = V.T @ H_mat @ V + + I = jnp.eye(D) + S_reg_ext = jnp.where(new_mask, S_reg_ext, I) + H_reg_ext =jnp.where(new_mask, H_reg_ext, I) - # Note: S_reg and H_reg will be of shape (max_D_out, max_D_out) - return S_reg, H_reg + # S_reg and H_reg will be a block matrix of shape (D, D): + # upper left: k x k block of regularized matrices S_reg, H_reg, respectively; + # lower right: (D-k) x (D-k) identity; + # upper right and lower left are zeros + return S_reg_ext, H_reg_ext + + +@jax.jit +def generalized_eigh_jitted(A, B): + # Solves the generalized eigenvalue problem A*v = lambda*B*v + # for symmetric/hermitian matrices A and symmetric positive-definite B. + + # Compute Cholesky decomposition of B (L*L.T) + # The 'lower=True' parameter is typically the default but good to be explicit. + L = jax.scipy.linalg.cholesky(B, lower=True) + + # Compute L_inv * A * L_inv_T. + # Use solve_triangular for efficiency and numerical stability instead of explicit inverse. + A_prime = jax.scipy.linalg.solve_triangular(L, A, lower=True) + A_prime = jax.scipy.linalg.solve_triangular(L.T, A_prime.T, lower=False).T + + # Solve the standard eigenvalue problem for A' + eigenvalues, eigenvectors_prime = jax.scipy.linalg.eigh(A_prime) + + # Transform eigenvectors back to the original basis (v = L_inv_T * v_prime) + eigenvectors = jax.scipy.linalg.solve_triangular(L.T, eigenvectors_prime, lower=False) + + # Normalize eigenvectors if needed (optional) + eigenvectors = eigenvectors / jnp.linalg.norm(eigenvectors, axis=0) + + return eigenvalues, eigenvectors def lanczos_alg(H, D, operand_prep, mes_kwargs={}, cutoff=1e-2, show_info=False): @@ -415,31 +459,29 @@ def operand_prep(): """ - - if check_for_tracing_mode(): - raise NotImplementedError("Solving the generalized eigenvalue problem not implemented in JAX 0.6") unitaries, coeffs = H.unitaries() # Step 1: Quantum Lanczos: Get expectation values of Chebyshev polynomials Tk_expvals = lanczos_expvals(H, D, operand_prep, mes_kwargs) - """ + if check_for_tracing_mode(): # Step 2: Build matrices S and H - S, H_mat = build_S_H_from_Tk_jax(Tk_expvals, D) + S, H_mat = build_S_H_from_Tk_jitted(Tk_expvals, D) # Step 3: Regularize matrices via thresholding - S_reg, H_reg = regularize_S_H_jax(S, H_mat, D, cutoff=cutoff) + S_reg, H_reg = regularize_S_H_jitted(S, H_mat, cutoff=cutoff) # Step 4: Solve generalized eigenvalue problem $\mathbf{H}\vec{v}=\epsilon\mathbf{S}\vec{v}$ - evals, evecs = jax.scipy.linalg.eigh(H_reg, S_reg) # Solving the generalized eigenvalue problem not implemented in JAX 0.6 + #evals, evecs = jax.scipy.linalg.eigh(H_reg, S_reg) # Solving the generalized eigenvalue problem not implemented in JAX 0.6 + evals, evecs = generalized_eigh_jitted(H_reg, S_reg) ground_state_energy = jnp.min(evals) * jnp.sum(coeffs) return ground_state_energy - """ + # Step 2: Build matrices S and H S, H_mat = build_S_H_from_Tk(Tk_expvals, D) From a1fa6575929123924b3e5038411cc6410af3ca35 Mon Sep 17 00:00:00 2001 From: renezander90 Date: Sun, 30 Nov 2025 23:55:12 +0000 Subject: [PATCH 32/32] Added test for jasp Lanczos --- .../test_lanczos_algorithm.py | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/tests/algorithms_tests/test_lanczos_algorithm.py b/tests/algorithms_tests/test_lanczos_algorithm.py index ca788cbe6..21a18beff 100644 --- a/tests/algorithms_tests/test_lanczos_algorithm.py +++ b/tests/algorithms_tests/test_lanczos_algorithm.py @@ -19,6 +19,7 @@ import networkx as nx import numpy as np from qrisp import QuantumVariable +from qrisp.jasp import jaspify from qrisp.lanczos import lanczos_alg from qrisp.operators import X, Y, Z from qrisp.vqe.problems.heisenberg import create_heisenberg_init_function @@ -49,4 +50,36 @@ def init_state(): energy, info = lanczos_alg(H, D, init_state, show_info=True) print(f"Ground state energy estimate: {energy}") + assert np.abs(energy_exact-energy) < 1e-2 + + +def test_jasp_lanczos_algorithm(): + + @jaspify(terminal_sampling=True) + def main(): + + L = 6 + G = nx.Graph() + G.add_edges_from([(k, (k+1) % L) for k in range(L - 1)]) + + # Define Hamiltonian e.g. Heisenberg with custom couplings + H = (1/4)*sum((X(i)*X(j) + Y(i)*Y(j) + 0.5*Z(i)*Z(j)) for i,j in G.edges()) + + energy_exact = H.ground_state_energy() + + # Prepare initial state function (tensor product of singlets) + M = nx.maximal_matching(G) + U_singlet = create_heisenberg_init_function(M) + + def init_state(): + qv = QuantumVariable(H.find_minimal_qubit_amount()) + U_singlet(qv) + return qv + + D = 6 # Krylov dimension + energy = lanczos_alg(H, D, init_state) + return energy_exact, energy + + energy_exact, energy = main() + assert np.abs(energy_exact-energy) < 1e-2 \ No newline at end of file