A Quantum Algorithm for the Finite Element Method

TL;DR

This work presents Qu-FEM, a fault-tolerant quantum algorithm for the finite element method that preserves geometric flexibility by introducing the unit of interaction $\hat{A}^{\parallel}_{jk}$ and the local-to-global node number indicator $\hat{A}^{\parallel}_{j}$. For constant elemental contributions, global FE arrays are assembled via compact block-encodings with gate costs scaling as $\tilde{\mathcal{O}}(d^2 p^2 n)$ in $d$ dimensions and element order $p$, while in the variable-coefficient case numerical integration on quantum hardware yields $\tilde{\mathcal{O}}(dn(D^d G^d + p^{2d}m))$ gate costs using Gauss–Legendre quadrature and MQET/QSP. Dirichlet boundary conditions are handled without modifying existing encodings by a boundary Lagrange-multiplier system, enabling a projector-based BE framework. The method achieves explicit BE circuits for Cartesian geometries and demonstrates Poisson and unidirectional Navier–Stokes benchmarks, illustrating consistency with classical FE results. Overall, Qu-FEM offers a path to quantum advantage in time for high-dimensional problems and space for low-dimensional, non-constant-coefficient cases, while outlining practical challenges such as preconditioning and observable retrieval.

Abstract

The finite element method (FEM) is a cornerstone numerical technique for solving partial differential equations (PDEs). Here, we present $\textbf{Qu-FEM}$, a fault-tolerant era quantum algorithm for the finite element method. In contrast to other quantum PDE solvers, Qu-FEM preserves the geometric flexibility of FEM by introducing two new primitives, the unit of interaction and the local-to-global indicator matrix, which enable the assembly of global finite element arrays with a constant-size linear combination of unitaries. We study the modified Poisson equation as an elliptic problem of interest, and provide explicit circuits for Qu-FEM in Cartesian domains. For problems with constant coefficients, our algorithm admits block-encodings of global arrays that require only $\tilde{\mathcal{O}}\left(d^2 p^2 n\right)$ Clifford+$T$ gates for $d$-dimensional, order-$p$ tensor product elements on grids with $2^n$ degrees of freedom in each dimension, where $n$ is the number of qubits representing the $N=2^n$ discrete grid points. For problems with spatially varying coefficients, we perform numerical integration directly on the quantum computer to assemble global arrays and force vectors. Dirichlet boundary conditions are enforced via the method of Lagrange multipliers, eliminating the need to modify the block-encodings that emerge from the assembly procedure. This work presents a framework for extending the geometric flexibility of quantum PDE solvers while preserving the possibility of a quantum advantage.

Figures (12)

• Figure 1: Circuit symbols for the Clifford $+$$T$ gate set. This gate set is used as a universal gate set for fault-tolerant quantum computing.
• Figure 2: The Toffoli gate expressed in terms of Clifford $+$$T$ gates, as given by nielsen2010quantum.
• Figure 3: Two-dimensional tensor product element.
• Figure 4: A two-dimensional mesh of first-order quadrilateral Lagrange elements with $\textnormal{numel}{} = 9$ elements, $\textnormal{nen}{} = 4$ elemental nodes, and $\textnormal{numnp}{} = 16$ total nodal points. A representative element is highlighted in orange. The local and global node numbers are written in binary, and the numbering convention is induced from the tensor product of one-dimensional elements.
• Figure 5: Comparison between classical and quantum strategies for the parallel assembly of finite element arrays. On top (a), we depict how a typical classically assembly process takes place in parallel. An array of classical processors (CPU's) computes the elemental contributions from different parts of the mesh, and inserts the contributions directly into the global array in memory. On the bottom (b), we depict our proposed method for the quantum assembly in parallel. Instead of partitioning the mesh into groups of elements (i.e., domain decomposition), we group local $ij$ interactions and compute the contributions arising from all such interactions simultaneously using a quantum computer (depicted abstractly as a 'QPU'). All pairs of interactions are then summed together in an LCU to obtain the assembled global array. As an example, we highlight the interactions between local node numbers $0$ and $3$ in (blue and red, respectively) in a two-dimensional mesh. The sparsity pattern for the corresponding contributions to the global array are shown in the right-hand matrix.

Theorems & Definitions (53)

• Definition 2.1: $(\alpha,m,\epsilon)$-Block-Encoding gilyenQuantumSingularValue2019alin2022lecture
• Lemma 2.2: LCU of diagonal position matrix, adapted from mukhopadhyay2024PauliFierz
• Definition 2.3: Quantum Linear Systems Problem (QLSP)
• Definition 2.4: Standard One-Dimensional Lagrange Finite Element of Order $p$
• Lemma 2.5: Partition of unity
• proof
• Definition 2.6: Tensor Product Element
• Proposition 2.7
• proof
• Definition 3.1: Unit of Interaction