Quantum algorithm for solving generalized eigenvalue problems with application to the Schrödinger equation

TL;DR

This work tackles the challenge of efficiently computing excited-state spectra for quantum Hamiltonians, where classical methods suffer from the curse of dimensionality. It introduces a quantum landscape-scanning algorithm that encodes a grid of matrix parameters into a quantum superposition and uses Quantum Phase Estimation with Amplitude Amplification to identify parameter values where a parameterized matrix has near-zero singular values; energy information is embedded in the quantum state, not in the energy register. The method is applied to the Schrödinger equation via a pseudospectral (collocation) framework, replacing explicit matrix inversion with a residue-based scan and achieving favorable scaling in problem size $N$ (up to $\widetilde{O}(\sqrt{N})$ versus $\widetilde{O}(N)$ classically) for certain well-behaved potentials, with extensions to rectangular generalized eigenvalue problems. The paper also analyzes classical versus quantum resource requirements, discusses the regimes where quantum advantage is plausible for high-dimensional, dense spectra, and outlines practical considerations for implementing the approach on fault-tolerant quantum hardware.

Abstract

Accurate computation of multiple eigenvalues of quantum Hamiltonians is essential in quantum chemistry, materials science, and molecular spectroscopy. Estimating excited-state energies is challenging for classical algorithms due to exponential scaling with system size, posing an even harder problem than ground-state calculations. We present a quantum algorithm for estimating eigenvalues and singular values of parameterized matrix families, including solving generalized eigenvalue problems that frequently arise in quantum simulations. Our method uses quantum amplitude amplification and phase estimation to identify matrix eigenvalues by locating minima in the singular value spectrum. We demonstrate our algorithm by proposing a quantum-computing formulation of the pseudospectral collocation method for the Schrödinger equation. We estimate fault-tolerant quantum resource requirements for the quantum collocation method, showing favorable scaling in the size of the problem $N$ (up to $\widetilde{\mathcal{O}}(\sqrt{N})$) compared to classical implementations with $\widetilde{\mathcal{O}}(N)$, for certain well-behaved potentials. Additionally, unlike the standard collocation method, which results in a generalized eigenvalue problem requiring matrix inversion, our algorithm circumvents the associated numerical instability by scanning a parameterized matrix family and detecting eigenvalues through singular value minimization. This approach is particularly effective when multiple eigenvalues are needed or when the generalized eigenvalue problem involves a high condition number. In the fault-tolerant era, our method may thus be useful for simulating high-dimensional molecular systems with dense spectra involving highly excited states, such as those encountered in molecular photodynamics or quasi-continuum regimes in many-body and solid-state systems.

Figures (7)

• Figure 1: A sketch of the generic landscape of the eigenvalues for a one-parameter family of Hermitian matrices $\mathbf M(\alpha)$. Here, for all equidistant grid values of $\alpha \in A = A_1 \cup A_2 \cup A_3$, the matrices $\mathbf M(\alpha)$ have eigenvalue $\lambda_0$ with accuracy $\varepsilon$.
• Figure 2: A quantum circuit realizing the search for eigenvalues of a matrix that are close to the target value $\lambda_0$, effectively constructing the algorithm described in Theorem \ref{['thm:quantum_algorithm_square']} for $N = 2^n$ and $K = 2^k$. The costs of the block-encoding are assumed to depend on the technique, and for example, for $d$-sparse they will scale as $\mathcal{O}(\log NK+1+\beta)$, where $\beta$ is the precision for the elements of the matrix $\mathbf M$ in the chosen block-encoding. The entire circuit (a) is a modification of the standard Grover algorithm. First, via $c$ steps of Quantum Phase Estimation, we prepare a state in an equal superposition of the eigenvectors and the corresponding eigenvalues on $c+1$ registers of $b$ qubits each. The "median" register stores the best approximation to the true eigenvalues Kerzner_2024. Additionally, we apply the region oracle on the uppermost qubit, which gives $\ket{1}$ if the corresponding eigenvalue is in the region of interest and $\ket{0}$ otherwise. This part of the procedure is called $\hat{\mathcal{F}}$, depicted on part (b). Observe that its action is equivalent to the tensor product of Hadamard gates in the standard Grover algorithm. Having prepared the state, in the first round, the oracle $Z$ acting on the qubit with the characteristic function $\chi$ is applied. What follows is the reflection over the equal superposition state, realized via $\hat{\mathcal{F}}^\dagger\hat{\mathcal{R}}_0\hat{\mathcal{F}}$, where $\hat{\mathcal{R}}_0 = 2\ket{0}\!\!\bra{0}^{\otimes 2(n+k)+b(c+1)+1} - \mathbb{I}$. Finally, after $\xi$ such rounds, the qubits are measured. The measurement of the $k$ bottom qubits, corresponding to the eigenvector registers, extracts the information about the relevant eigenvector $\ket{\psi_i} = \ket{\alpha_i}_k\otimes\ket{\phi_i}_n$ of matrix $\widetilde{\mathbf M}$, as defined in Eq. \ref{['eq:firstMtylda']}. Therefore, this determines the value $\alpha_i$ that yields matrices $\mathbf{M}(\alpha_i)$ with eigenvalues in the target region.
• Figure 3: Comparison of the matrix inverse method described in Sec. \ref{['subsec:collocation_method']} with the classical landscape scanning method for solving the collocation Schrödinger equation for the 1D harmonic oscillator potential, using 26 basis functions and 80 grid points. The condition number of $\mathbf B^\dagger \mathbf B$ equals to $\kappa = 1.53\times 10^{9}$. The exact energies $E_n = 2n+1$ are denoted by solid, vertical green lines. The blue curve corresponds to the minimal singular value produced by the landscape method, while the dashed vertical lines correspond to the eigenvalue estimates obtained by the inversion method.
• Figure 4: The same setup as in Fig. \ref{['fig:26functions']} with 35 (blue) and 36 (orange) basis functions, with 80 grid points. The condition numbers respectively $1.21\times 10^{16}$ and $5.47\times 10^{16}$. Two different sets of basis functions are combined to account for the parity of the solutions to the Schrödinger equation.
• Figure 5: Comparison of the matrix inverse method described in Sec. \ref{['subsec:collocation_method']} with the classical landscape scanning method for solving the collocation Schrödinger equation for the 1D Morse potential, given by Eq. \ref{['eq:Morse_potential']}, using 35 basis functions and 80 grid points, with the condition number of $\mathbf B^\dagger \mathbf B$ equal to $1.23 \times 10^{16}$. The curves have the same meaning as in Fig. \ref{['fig:26functions']}. Note that the matrix inversion method for solving the collocation equations (red dashed line) fails to correctly predict even a single energy level (green solid lines).

Theorems & Definitions (5)

• Definition 1: Hogben_2013
• Theorem 1
• proof
• Corollary 2
• proof