⚛️ quantum physics

Quantum algorithms for the exponentiation of Toeplitz matrices and applications in partial differential equations

This paper presents quantum algorithms that circumvent the large-norm limitations of banded Toeplitz matrices by leveraging their relationship to circulant and skew-circulant generators to efficiently construct block encodings for matrix exponentiation, which are then applied to solve discretized heat equations with various boundary conditions.

Original authors: Xabier Gutiérrez, Nicola Mariella, Javier González-Conde, Sergiy Zhuk, Mikel Sanz

Published 2026-09-28

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Xabier Gutiérrez, Nicola Mariella, Javier González-Conde, Sergiy Zhuk, Mikel Sanz

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Science often deals with equations that describe how things change over time, from the flow of heat through a metal rod to the movement of fluids in the atmosphere. These are known as partial differential equations, and they are the language of physics and engineering. To solve them on a computer, scientists break the continuous world into a grid of tiny points, turning the smooth equations into massive lists of numbers. The solution to these problems usually involves a mathematical operation called exponentiation, which tells us how the system evolves from a starting point to a future moment. For decades, the hope has been that quantum computers could solve these problems much faster than classical machines, offering a speedup that grows exponentially with the size of the problem. However, a significant roadblock has stood in the way: the standard way to prepare these calculations on a quantum computer requires a "normalization" step that becomes impossibly expensive as the grid gets finer. The numbers involved in the equations grow so large that the quantum computer struggles to handle them, effectively canceling out the potential speed advantage.

A team of researchers has now developed a new method to bypass this roadblock, specifically for a common type of matrix that appears in these grid-based calculations. These matrices, known as Toeplitz matrices, have a special repeating pattern where the numbers along any diagonal are identical. While these patterns are crucial for modeling physical systems, they are notoriously difficult to work with on quantum computers because they cannot be easily broken down into simpler parts. The researchers found a way to rewrite these complex matrices as combinations of two simpler, rotating structures that are much easier for a quantum computer to handle. By doing this, they created a direct path to calculate the time evolution of the system without needing the expensive normalization step that usually slows things down.

The core of their discovery lies in how they treat the mathematical building blocks of these matrices. Instead of trying to force the quantum computer to handle the difficult, non-repeating parts directly, the team showed that these difficult parts can be expressed as a sum of two types of shifting patterns. One type shifts information in a circle, like beads on a necklace, while the other shifts them with a slight twist. Both of these patterns have a special property: they can be perfectly understood by a quantum computer using a tool called the Quantum Fourier Transform, which acts like a prism that separates light into its individual colors, but here it separates the complex numbers into their fundamental frequencies. Because these patterns are so well-behaved, the researchers could approximate their behavior using a series of simple, controlled rotations on individual quantum bits.

To make this practical, the team introduced a method to cut off the parts of the calculation that contribute very little to the final answer. In many physical systems, such as the diffusion of heat, the most important information is concentrated in the low-frequency parts of the signal, while the high-frequency parts fade away quickly. By focusing only on the significant low-frequency components and ignoring the rest, the researchers could drastically reduce the size of the calculation while keeping the error under strict control. This allowed them to construct a simplified version of the time-evolution operator that is small enough to be handled efficiently, yet accurate enough to be useful. They then combined these simplified pieces using a step-by-step approach, similar to taking small steps to walk a long distance, to reconstruct the full solution.

The researchers tested this framework on the classic problem of the heat equation, which describes how heat spreads through a material. They showed that their method works for different types of boundaries, including cases where the material is a loop, where the ends are held at a fixed temperature, or where the ends are insulated. In each case, they demonstrated that the new approach avoids the massive scaling costs that plague previous methods. Instead of the computational cost exploding as the grid becomes finer, their method keeps the cost manageable. This is a significant step forward because it removes the normalization bottleneck that has prevented quantum computers from solving these specific types of physics problems efficiently.

While the method is powerful, the authors are careful to note its limits. The approach works best when the repeating pattern in the matrix is narrow compared to the total size of the system, a condition that is common in many physical simulations but not universal. They also point out that while the error bounds are well-defined, the exact number of steps needed to reach a certain level of precision depends on the specific coefficients of the problem. Furthermore, the selection of which parts of the calculation to keep is currently based on observed patterns rather than a strict mathematical proof for every possible case. Despite these open questions, the work provides a clear and concrete pathway for quantum computers to tackle a class of problems that were previously out of reach, turning a theoretical possibility into a practical algorithm for simulating the physical world.

Technical Summary: Quantum Algorithms for the Exponentiation of Toeplitz Matrices

Problem Statement
Partial Differential Equations (PDEs) are fundamental to modeling physical phenomena, often reduced via spatial discretization to the time evolution of a linear system $\partial_t u = Au$ , with the solution $u(t) = e^{At}u(0)$ . The matrix $A$ arising from such discretizations is typically a banded Toeplitz or quasi-Toeplitz matrix. While quantum computers offer potential exponential advantages for solving such problems, a significant bottleneck exists: standard quantum algorithms rely on unitary transformations, whereas the matrix exponential $e^{At}$ for PDE operators is generally non-unitary.

Existing approaches utilize the block-encoding framework to embed a subnormalized matrix $A/\alpha$ into a larger unitary. However, for differential operators like the Laplacian, the norm $\|A\|$ diverges as the grid spacing $\delta_x$ decreases (e.g., $\|\Delta\| \sim 4/\delta_x^2$ ). This necessitates a large subnormalization factor $\alpha$ , which scales the complexity of algorithms based on Quantum Signal Processing (QSP) or Quantum Singular Value Transformation (QSVT) prohibitively. Furthermore, directly block-encoding the exponential $e^{At}$ has not been previously addressed, as arbitrary Toeplitz matrices are not generally diagonalizable.

Methodology
The authors propose a quantum framework that circumvents the normalization bottleneck by directly constructing a block encoding of the exponential of banded Toeplitz matrices. The methodology relies on exploiting the algebraic structure of these matrices:

Decomposition into Diagonalizable Generators:
Although a general Toeplitz matrix $T$ is not diagonalizable, it can be expressed as a polynomial of the lower ( $L_n$ ) and upper ( $L_n^\top$ ) shift operators. These shift operators are nilpotent and non-diagonalizable. However, the authors relate them to the circulant generator $C_n$ and the skew-circulant generator $N_n$ . Both $C_n$ and $N_n$ are unitary and diagonalizable by the Quantum Fourier Transform (QFT). Specifically, the shift operators can be written as linear combinations of powers of $C_n$ and $N_n$ .
Approximated Block Encoding via Truncated Pauli Strings:
The eigenvalues of the circulant generator $C_n$ form a diagonal matrix $\Gamma_n$ . The authors introduce an efficient approximation of $\Gamma_n$ by truncating the smaller phase components of its eigenvalues. This approximation, $\tilde{\Gamma}_{d,n}$ , is decomposed into a sum of commuting Pauli- $Z$ strings with a closed-form Frobenius-norm error bound. Because the resulting Pauli strings commute, their exponentiation reduces to a product of single-qubit rotations, which can be implemented efficiently.
QFT-Trotter Product Formula:
To handle the general non-diagonalizable Toeplitz matrix, the authors apply the Lie-Trotter product formula. By decomposing the Toeplitz matrix into a sum of circulant and skew-circulant components (which are diagonalized by QFT), the exponential $e^{Tt}$ is approximated by the product of the exponentials of these diagonalizable components. This yields a specialized QFT-Trotter product formula.
Application to the Heat Equation with Frequency Cutoff:
As a concrete application, the framework is applied to the discretized heat equation ( $\partial_t u = \Delta u$ ). The discretized Laplacian $\Delta$ under periodic boundary conditions is circulant. The authors observe that the eigenvalues of the propagator $e^{\Delta t}$ decay rapidly away from low-frequency modes. By imposing a frequency cutoff $k_c$ , they truncate the operator to a narrow band of low frequencies. This truncation removes the $\delta_x^{-2}$ scaling of the norm.
- For the resulting non-unitary truncated propagator, the authors employ Linear Combination of Hamiltonian Simulation (LCHS) to express the evolution as a linear combination of skew-Hermitian unitaries.
- Each unitary in this combination is implemented using the truncated Pauli-string decomposition described earlier.
- This approach is extended to Dirichlet and Neumann boundary conditions by utilizing the Discrete Sine Transform (DST) and Discrete Cosine Transform (DCT), respectively, which are embedded as block encodings in an $(n+1)$ -qubit system.

Key Contributions and Results

Direct Exponentiation: The paper presents the first direct construction of a block encoding for the exponential of banded Toeplitz matrices, avoiding the intermediate step of encoding the matrix $A$ and then applying QSP/QSVT.
Efficient Decomposition: A novel, controllable truncation of the circulant eigenphases into Pauli strings is derived, providing a closed-form error bound that depends on the truncation depth $\delta = n-d$ rather than the system size $n$ .
QFT-Trotter Formula: A specialized product formula is established for exponentiating Toeplitz matrices by leveraging their decomposition into QFT-diagonalizable unitaries.
Elimination of Normalization Bottleneck: In the application to the heat equation, the frequency cutoff strategy allows the construction of a propagator with an $O(1)$ subnormalization factor. This contrasts sharply with direct block encodings of the Laplacian, where the subnormalization scales as $O(1/\delta_x^2)$ , making QSP/QSVT infeasible for fine grids.
Error Bounds: The total error is bounded by terms related to the Trotter steps and the truncation error, specifically $O(t^2/c \|\alpha\|_1^2) + O(\pi t m / (2\delta) \|\alpha\|_1)$ , where $c$ is the number of Trotter steps and $m$ is the bandwidth.

Significance
The paper claims that this work provides a general methodology for exploiting the algebraic structure of banded Toeplitz operators to overcome the prohibitive costs associated with subnormalization factors in PDE solvers. By directly addressing the exponentiation of these matrices, the method achieves $O(1)$ subnormalization, making it suitable for regimes where the operator norm is large (typical in differential operator discretizations) but the bandwidth is small relative to the system size ( $m \ll 2^n$ ).

The authors note that while prior works have addressed block encoding of Toeplitz operators or the Laplacian itself, they did not directly construct the exponential form $e^{Tt}$ required for time-evolution problems without incurring high subnormalization costs. This work fills that gap, offering a specialized, efficient alternative to general PDE solvers for structured matrices.

Limitations
The authors explicitly state two limitations:

The Trotter error bound in Theorem V.1 includes constants that may grow exponentially with $t\|\alpha\|_1$ for non-contractive coefficient sets; a tighter bound for general banded Toeplitz operators remains an open question.
The selection of the dominant subset of Pauli strings ( $\Omega$ ) used in the approximation is currently heuristic, based on observed patterns. A rigorous a-priori bound on the size of $\Omega$ versus the target error is needed to fully strengthen the result.

Technical Summary: Quantum Algorithms for the Exponentiation of Toeplitz Matrices

More like this