Quantum Rational Transformation Using Linear Combinations of Hamiltonian Simulations

Rational functions are exceptionally powerful tools in scientific computing, yet their abilities to advance quantum algorithms remain largely untapped. In this paper, we introduce effective implementations of rational transformations of a target operator on quantum hardware. By leveraging suitable integral representations of the operator resolvent, we show that rational transformations can be performed efficiently with Hamiltonian simulations using a linear-combination-of-unitaries (LCU). We formulate two complementary LCU approaches, discrete-time and continuous-time LCU, each providing unique strategies to decomposing the exact integral representations of a resolvent. We consider quantum rational transformation for the ubiquitous task of approximating functions of a Hermitian operator, with particular emphasis on the elementary signum function. For illustration, we discuss its application to the ground and excited state problems. Combining rational transformations with observable dynamic mode decomposition (ODMD), our recently developed noise-resilient quantum eigensolver, we design a fully real-time approach for resolving many-body spectra. Our numerical demonstration on spin systems indicates that our real-time framework is compact and achieves accurate estimation of the low-lying energies.

💡 Research Summary

The paper introduces a novel framework for implementing rational transformations of Hermitian operators on quantum hardware by leveraging linear combinations of Hamiltonian simulations (LCHS). Traditional quantum algorithms for function approximation, such as Quantum Signal Processing (QSP) and Quantum Singular Value Transformation (QSVT), rely on polynomial approximations. While powerful for smooth functions, these methods become inefficient for non‑smooth or discontinuous functions like the signum function, requiring high‑degree polynomials that translate into deep circuits and significant error accumulation.

Rational functions, expressed as ratios of polynomials, can be decomposed into a sum of partial fractions, each term being a resolvent ((z_k - H)^{-1}) where (z_k) are poles. The core technical contribution is to represent each resolvent using unitary time‑evolution operators (e^{-iHt}) through suitable integral kernels. For complex poles with non‑zero imaginary parts, a Laplace‑type representation reduces the double integral to a single integral with an exponentially decaying factor, enabling short‑time simulations. For real poles, a Gaussian kernel introduces rapid decay, allowing the integral to be approximated efficiently despite the lack of intrinsic exponential damping.

Two complementary approximation strategies are developed:

1. Discrete‑time LCHS – An optimal quadrature rule (e.g., Gauss‑Legendre) selects a set of time points ({t_j}) and corresponding weights ({w_j}). The resolvent is approximated by a linear combination (\sum_j w_j e^{-iHt_j}). By carefully choosing the nodes according to the weight function of the integral kernel, the method achieves high accuracy with a number of samples proportional to the square root of the desired precision, matching the optimal query complexity for linear‑combination‑of‑unitaries constructions.

2. Continuous‑time LCHS – The weight function itself is expanded in a finite basis of Gaussian functions. Each basis integral can be evaluated analytically or via Monte‑Carlo sampling, yielding a weighted sum of unitary evolutions. This approach is especially advantageous when the kernel is naturally Gaussian, as it avoids explicit discretization of the time domain and can exploit variance‑reduction techniques.

Both approaches maintain favorable scaling of the maximum simulation time (the longest individual Hamiltonian evolution) and the total simulation time (the sum over all evolutions). The authors also discuss an ancillary bosonic mode (e.g., a harmonic oscillator) that can encode the evolution time, allowing the physical runtime to remain constant at the expense of more demanding ancilla preparation.

The framework is applied to construct a high‑quality approximation of the signum function, a bounded discontinuous function that acts as a spectral projector: (\operatorname{sgn}(H) = \sum_k c_k (z_k - H)^{-1}). By combining the signum approximation with a spectral filter (\frac{1}{2}(I + \operatorname{sgn}(H - \mu))), the algorithm isolates eigenvalues below a chosen threshold (\mu). This enables efficient extraction of ground‑state and low‑lying excited‑state energies.

Numerical experiments on several spin models (1‑D Heisenberg chains, 2‑D transverse‑field Ising models) demonstrate that the rational‑transformation approach achieves comparable or better accuracy than polynomial‑based methods while using significantly fewer Hamiltonian simulation steps. The authors further integrate their rational filter with Observable Dynamic Mode Decomposition (ODMD), a noise‑resilient eigensolver that processes time‑series data of observable expectations. The combined “Rational‑LC‑ODMD” pipeline operates entirely in real time, requiring only short‑depth circuits and modest numbers of measurement shots, and exhibits robustness against realistic noise models.

The paper concludes that rational transformations, when expressed via LCHS, provide a powerful and hardware‑friendly toolbox for quantum algorithm designers. It opens avenues for tackling interior‑eigenvalue problems, matrix functions with branch cuts, and other tasks where polynomial approximations are inefficient. Future work may explore extensions to degenerate poles, adaptive kernel selection, and experimental implementation on both digital and analog quantum platforms, as well as error‑mitigation strategies tailored to the linear‑combination structure.