diff --git a/documentation/source/reference/Algorithms/Lanczos.rst b/documentation/source/reference/Algorithms/Lanczos.rst new file mode 100644 index 000000000..da39f4675 --- /dev/null +++ b/documentation/source/reference/Algorithms/Lanczos.rst @@ -0,0 +1,17 @@ +.. _lanczos_alg: + + +Lanczos Algorithm +================= + +.. currentmodule:: qrisp.lanczos + +.. autofunction:: lanczos_alg + +.. autosummary:: + :toctree: generated/ + + lanczos_expvals + inner_lanczos + 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.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.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.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/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/documentation/source/reference/Algorithms/index.rst b/documentation/source/reference/Algorithms/index.rst index 759a916ac..0e32c4593 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 Algorithm ` + - finding the ground state energy of a Hamiltonian * - :ref:`CKS Algorithm ` - solving the quantum linear system problem @@ -45,6 +47,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 10a5b0c71..085c934d2 100644 --- a/src/qrisp/__init__.py +++ b/src/qrisp/__init__.py @@ -46,6 +46,7 @@ "qite", "qmci", "cks", + "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 923a01302..cd6901374 100644 --- a/src/qrisp/algorithms/__init__.py +++ b/src/qrisp/algorithms/__init__.py @@ -25,4 +25,5 @@ import qrisp.algorithms.vqe as vqe import qrisp.algorithms.qite as qite import qrisp.algorithms.qmci as qmci +import qrisp.algorithms.lanczos as lanczos import qrisp.algorithms.cks as cks diff --git a/src/qrisp/algorithms/lanczos.py b/src/qrisp/algorithms/lanczos.py new file mode 100644 index 000000000..3e242d454 --- /dev/null +++ b/src/qrisp/algorithms/lanczos.py @@ -0,0 +1,510 @@ +""" +******************************************************************************** +* 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 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 + + +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. + + 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, \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 of the prepared quantum states 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)$. + operand_prep : callable + Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., + ``operand=operand_prep()``. + + Returns + ------- + QuantumVariable + A QuantumVariable for which the measured statistics encodes the expectation values $\langle T_k(H)\rangle$. + + """ + + # Extract unitaries and size of the case_indicator QuantumFloat. + U, state_prep, n = H.pauli_block_encoding() + + 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}$. + + case_indicator = QuantumFloat(n) + operand = operand_prep() + + 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) + 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 + + return q_cond(k%2==0, even, odd, case_indicator, operand, k) + + +def compute_expectation(meas_res): + r""" + Convert measurement results into an expectation value. + + 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-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 block-encoding unitary $U$ and then measure the + auxilary QuantumVariable in the $Z$ basis. + + Parameters + ---------- + 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 meas_res.items(): + + if int(outcome) == 0: + expval += prob * 1 + else: + expval += prob * (-1) + return expval + + +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. + + 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)$. + operand_prep : callable + Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., + ``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. + + 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: + mes_kwargs["shots"] = 100000 + + 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 = mes_kwargs["shots"], post_processor = post_processor) + expvals = jnp.zeros(2*D) + + for k in range(0, 2*D): + 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, operand_prep) + meas = qarg.get_measurement(**mes_kwargs) + expvals[k] = compute_expectation(meas) + + return expvals + +# Postprocessing + +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. + + Using Chebyshev recurrence relations, this function generates the matrix elements for + 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 : ndarray + Dictionary of expectations $⟨T_k(H)⟩$ for each Chebyshev polynomial order $k$. + D : int + Krylov space dimension. + + Returns + ------- + S : ndarray + Overlap (Gram) matrix $\mathbf{S}$ for Krylov states. + 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 + + +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. + + This function applies a spectral cutoff: only directions in the Krylov subspace with eigenvalues + 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 : ndarray + Overlap matrix. + H_mat : ndarray + Hamiltonian matrix. + cutoff : float + Eigenvalue threshold for regularizing $\mathbf{S}$. + + Returns + ------- + S_reg : ndarray + Regularized overlap matrix. + 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_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] + + 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=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) + + # 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) + + # 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): + 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" `_. + + 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 $\langle T_k(H)\rangle$. + 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 : QubitOperator + Hamiltonian for which to estimate the ground-state energy. + D : int + Krylov space dimension. + operand_prep : callable + Function returning the (operand) QuantumVariable in the initial system state $\ket{\psi_0}$, i.e., + ``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. + cutoff : float + Regularization cutoff threshold for overlap matrix $\mathbf{S}$. The default is 1e-2. + show_info : bool, optional + If True, a dictionary with detailed information is returned. The default is False. + + Returns + ------- + energy : float + Estimated ground state energy of the Hamiltonian H. + info : dict, optional + Full details including: + - '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 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 + 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()) + + 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) + + 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, operand_prep, show_info=True) + + print(f"Ground state energy estimate: {energy}") + + + """ + + 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_jitted(Tk_expvals, D) + + # Step 3: Regularize matrices via thresholding + 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 = 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) + + # Step 3: Regularize matrices via thresholding + 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) + + ground_state_energy = np.min(evals) * np.sum(coeffs) + + + results = { + 'Tk_expvals': Tk_expvals, + 'energy': ground_state_energy, + 'eigvals': evals, + 'eigvecs': evecs, + 'S_reg': S_reg, + 'H_reg': H_reg, + } + + if show_info: + return ground_state_energy, results + else: + return ground_state_energy \ No newline at end of file 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() diff --git a/tests/algorithms_tests/test_lanczos_algorithm.py b/tests/algorithms_tests/test_lanczos_algorithm.py new file mode 100644 index 000000000..21a18beff --- /dev/null +++ b/tests/algorithms_tests/test_lanczos_algorithm.py @@ -0,0 +1,85 @@ +""" +******************************************************************************** +* 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.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 + + +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_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) + 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_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 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)