quantum linear system, quantum algorithm, domain decomposition, partial differential equation, finite element method

1 Introduction↩︎

The numerical solution of partial differential equations, or PDEs, plays a central role in a wide range of scientific and engineering applications, including fluid dynamics, electromagnetism, and solid mechanics. Standard approaches rely on spatial discretization, which typically leads to large-scale sparse linear systems whose size grows rapidly as higher accuracy is sought. For fine-grid discretizations, the resulting linear algebra problems can become prohibitively expensive to solve using classical (direct or iterative) methods, due to both computational complexity and memory constraints. This scalability bottleneck has motivated the exploration of alternative computational paradigms, including quantum algorithms for linear systems and PDEs. In [1], Harrow, Hassidim, and Lloyd introduced the first quantum algorithm for solving linear systems of equations (the HHL algorithm). The complexity of quantum linear system solvers (QLSS) has since been improved to a dependency that is linear in the matrix condition number and logarithmic in the inverse of the precision see e.g., [2]–[8] and [9] for a survey of the different propositions. Variational approaches to quantum linear system solving have also been proposed for the Noisy Intermediate-Scale Quantum regime, known as NISQ, in which quantum resources are limited [10].

Quantum linear solvers have already been applied to the linear systems arising from the finite element discretization of partial differential equations. It has been shown in [11] that, for smooth solutions of the Poisson equation, a quantum advantage can be achieved for large problem dimensions. As in the classical setting, the performance of QLSSs depends on the condition number of the system matrix. In the finite element context, the condition number typically grows with the problem size, which may ultimately negate the expected quantum speedup.

In classical numerical linear algebra, this issue is typically fixed by preconditioning: a preconditioner \(H \approx A^{-1}\) is applied to improve the conditioning of the original system \(Ax = b\), leading to the preconditioned system \(HA x = Hb\). If \(H\) is well chosen, the preconditioned system is solved more efficiently than the original one. Quantum preconditioning strategies have been proposed by [12]–[17]. In particular, [13] introduces a QLSS using a sparse approximate inverse (SPAI) preconditioner, while [16] studies a circulant preconditioning approach.

Recently, [18] proposed a quantum linear system solver based on the multilevel Bramble–Pasciak–Xu (BPX) preconditioner [19] via a block-encoding of the preconditioned system in its symmetrized form. The BPX preconditioning strategy has also been incorporated into the Schrödingerization framework [20], where singular value amplification is used to extract the solution. In a related direction, [21] develops quantum algorithms for the neutron diffusion \(k\)-eigenvalue problem using BPX-based techniques and fast inversion methods. Quantum linear system solvers have also been integrated into multiscale numerical homogenization frameworks to accelerate the computation of local fine-scale problems [22].

The present work is inspired by the, previously mentioned, article [18]. Whereas [18] considers a multigrid-type preconditioner [19], [23], we focus on preconditioning by domain decomposition. We borrow several ideas from [18] and also apply their quantum BPX algorithm within our own preconditioner. In classical numerical linear algebra, domain decomposition methods constitute a well-established framework for constructing scalable preconditioners for elliptic problems [24]. They rely on a decomposition of the computational domain into subdomains, where local problems are solved and combined through a coarse space correction to obtain global convergence. Classical theory shows that, with appropriately designed coarse spaces, such methods yield condition number bounds that are robust with respect to mesh size and coefficient heterogeneities [24]–[30]. The contribution of this work is to show that domain decomposition preconditioners can be incorporated into a QLSS. This establishes the feasibility of implementing domain decomposition preconditioners in a quantum linear system framework, thereby complementing existing quantum preconditioning strategies based on multilevel, sparse approximate inverse and circulant techniques. We note for completeness that [31] utilizes domain decomposition concepts in variational quantum algorithms, albeit in a distinct setting and for a different objective to ours. Indeed, [31] introduces a domain decomposition based strategy to mitigate Barren plateau–like phenomena in variational quantum algorithms. The approach relies on localizing the cost function by partitioning the spatial domain into overlapping subdomains, each associated with a localized parametrized quantum circuit and measurement operator. This improves trainability compared to global variational formulations.

The article is organized as follows. Section 2 presents the finite element setting for the scalar elliptic problem. Section 3 introduces the concepts of domain decomposition and discusses the potential advantages of domain decomposition in the quantum framework (Section 3.4). A particular domain decomposition preconditioner, the two-level Additive Schwarz preconditioner is defined. Spectral bounds for the resulting preconditioned operator are stated in Theorem 1. Following ideas from [18], the preconditioned system is written in a symmetrized form. This way, the spectral bounds imply a bound on the condition number which is later needed for the analysis of the QLSS. In Section 4, we construct a block-encoding of the symmetric preconditioned system, describe its implementation and analyze the cost of inverting it by a quantum singular value transformation. Emphasis is put on the case where BPX is used to solve each local problem within the domain decomposition preconditioner. In Section 5 we present a proof-of-concept numerical simulation of the quantum circuits.

1.0.0.1 Notation

For any integer \(m\), we denote by \(I_m\) the identity matrix in \(\mathbb{R}^{m \times m}\). For any matrix \(A\), the spectral condition number is defined by \[\kappa(A)=\|A\|_2\,\|A^+\|_2 =\frac{\sigma_{\max}(A)}{\sigma_{\min}^+(A)},\] where \(A^+\) denotes the pseudoinverse and \(\sigma_{\min}^+(A)\) is the smallest nonzero singular value of \(A\). In the special case where \(A\) is a non-singular positive normal matrix (e.g., Hermitian positive definite), the singular values are also the eigenvalues so \[\kappa(A)=\frac{\lambda_{\max}(A)}{\lambda_{\min}(A)}.\] The conjugate transpose of \(A\) is denoted by \(A^\dagger\), its pseudoinverse is denoted by \(A^+\). For real matrices, the transpose of \(A\) is denoted by \(A^\top\).

2 Finite element setting↩︎

In the present article we focus on a particular linear system that arises from discretizing a simple PDE on a simple geometry. We are confident that the methods introduced could be extended to more difficult problems in the future.

2.1 Poisson equation↩︎

In detail, we solve the Poisson equation in a domain \(\Omega \in \mathbb{R}^d\) with homogeneous Dirichlet boundary condition. We make the assumption that \(\Omega\) is a rectangle (if \(d=2\)), cuboid (if \(d=3\)) or, more generally, a hyperrectangle (any \(d\)). The source term \(f\) is taken in the Lebesgue space \(L^2(\Omega)\) of functions that are square integrable over \(\Omega\). As is standard, we denote by \(H_0^1(\Omega)\) the Sobolev space of functions in \(L^2(\Omega)\) that vanish on the boundary of \(\Omega\) and have square-integrable first derivatives.

Then, the formulation of our problem is: Find \(u \in H_0^1(\Omega)\) such that \[\label{eq:varf} \int_\Omega \rho(x) \, \nabla u(x) \cdot \nabla v(x) \mathrm d x = \int_\Omega f(x) v(x) \mathrm d x \quad \forall v \in H_0^1(\Omega),\tag{1}\] where, \(\rho\) is into the space of \({d \times d}\) symmetric positive definite matrices with measurable entries. We assume that there exist \(\rho_{\min},\rho_{\max}>0\) such that \[\label{eq:assumption-rho} \rho_{\min} \eta^\top \eta \leq \eta^\top \rho(x) \eta \leq \rho_{\max} \eta^\top \eta ~ \text{for all } \eta \in \mathbb{R}^d \text{ and almost all } x \in \Omega.\tag{2}\] As is well-known, this is the weak formulation of the partial differential equation (PDE) \[\label{eq:pde} \left\{ \begin{align} - \nabla \cdot \rho \nabla u &= f \quad \text{in } \Omega, \\ u &= 0 \quad \text{on } \partial \Omega, \end{align} \right.\tag{3}\] in which \(\rho\) is the, possibly anisotropic, diffusion parameter. If \(\rho = I_{d}\), the first line in the equation rewrites as \(- \Delta u = f\).

2.2 Finite element setting↩︎

In order to discretize the PDE, we introduce a Cartesian mesh \(\mathcal{T}_L\) of \(\Omega\). Throughout the article, the subscript \(\cdot_L\) refers to quantities defined on the mesh \(\mathcal{T}_L\). We apply the \(\mathbb{Q}_1\) finite element method, meaning that we consider the restriction of the variational formulation 1 to a particular finite dimensional subspace \(V_L\) of \(H_0^1(\Omega)\). The \(\mathbb{Q}_1\) finite element space is the space of functions that are globally continuous over \(\Omega\), \(d\)-linear on each mesh element and that vanish on \(\partial \Omega\). Formally, we set \[V_L := \left\{u_L \in C^0(\bar{\Omega}) \cap H_0^1(\Omega);\, {u_L}_{|K} \in \mathbb{Q}_1(K), \, \forall K \in \mathcal{T}_L \right\},\] where \[\mathbb{Q}_1(K) = \operatorname{span}\left\{ x_1^{\alpha_1} \cdots x_d^{\alpha_d} :\alpha_i \in \{0,1\} \right\} \text{ with } x= (x_1,\dots,x_d) \text{ for } x \in K.\]

We let \(\{ \phi_j \}_{1 \leq j \leq \mathcal{N}}\) be the standard nodal basis of \(V_L\), i.e., the set of functions that are \(1\) at an interior vertex of the mesh and \(0\) at all others. We then obtain the linear system of equations \[\label{eq:Axequalb} A x = b,\tag{4}\] where \(A \in \mathbb{R}^{\mathcal{N} \times \mathcal{N}}\) and \(b \in \mathbb{R}^{\mathcal{N}}\) are defined by \[A_{ij} := \int_\Omega \rho \, \nabla \phi_i \cdot \nabla \phi_j \text{ for all }1 \leq i,j \leq \mathcal{N} \text{ and } b_i := \int_\Omega f \phi_i \text{ for all } 1 \leq i \leq \mathcal{N}.\] Note that the matrix \(A\) is symmetric positive definite and that its order, \(\mathcal{N}\), is the number of vertices of \(\mathcal{T}_L\) that are in the interior of \(\Omega\).

2.3 Factorized form of the linear system↩︎

We follow the idea from [18] to factorize \(A\) into discrete gradient operators weighted by the (matrix-valued) coefficient \(\rho\). We assume that \(\rho\) is constant over each mesh element \(K \in \mathcal{T}_L\).

We introduce the discontinuous piecewise \(\mathbb{Q}_1\) space \[\label{eq:QL} Q_L := \left\{u_L \in L^2(\Omega);\, {u_L}_{|K} \in \mathbb{Q}_1(K), \;\forall K \in \mathcal{T}_L \right\},\tag{5}\] where each element \(K\) carries \(2^d\) local degrees of freedom. We equip \(Q_L\) with an \(L^2\)-orthonormal basis consisting of \(2^d\) functions \(\psi_{K,1}, \dots, \psi_{K,2^d}\) on each element \(K \in \mathcal{T}_L\).

We next define \(Q_L^d := Q_L \otimes \mathbb{R}^d\), and equip it with the basis functions \(\psi_{K,k} \otimes e_s\) where \(e_s\) (\(s \in \{1, \dots, d\}\)) is the \(s\)-th canonical basis vector in \(\mathbb{R}^{d}\). Let also \(C_L\) denote the matrix representation of the (elementwise) gradient operator \[\nabla : V_L \to Q_L^d.\]

Then, as in [18], the stiffness matrix admits the decomposition \[A = C_L^\top (D_\rho \otimes I_{2^d}) C_L, \label{eq:decomposition95stiffness95matrix}\tag{6}\] where \(D_\rho\) is a block-diagonal matrix with one block per element \(K \in \mathcal{T}_L\), each block being the diffusion tensor \(\rho_{|K} \in \mathbb{R}^{d \times d}\). The tensor product structure reflects the assumption that \(\rho\) is constant over each element.

The factorized form of \(A\) plays a crucial role in block-encoding the linear system. Before moving to the quantum framework we introduce domain decomposition preconditioners.

3 Domain decomposition↩︎

This section focuses on the construction of a particular domain decomposition method: the two-level Additive Schwarz preconditioner with inexact local solves. The presentation is for a simple geometry but most definitions carry over to more complex geometries. The spectrum of the Poisson problem preconditioned by additive Schwarz is bounded by Theorem 1 (which is well-known in the domain decomposition literature [24]). A split, or factorized, form of the preconditioner is proposed that is well-suited to its quantum implementation. In Section 3.4, the advantages of domain decomposition are discussed as well as its potential for accelerating quantum linear solvers.

3.1 Geometry of the domain and subdomains↩︎

We introduce a partition of \(\Omega\) into \(N\) overlapping subdomains denoted \(\Omega^{(i)}\). The subdomains \(\Omega^{(i)}\) are open, mesh conforming, subsets of \(\Omega\) which satisfy \[\bigcup\limits_{i=1}^N \Omega^{(i)} = \Omega.\] Remark that, in order to satisfy the definition, it is clear that the subdomains must overlap, i.e., \(\Omega^{(i)} \cap \Omega^{(j)} \neq 0\) for some indices \(j\). It is natural to define local meshes \(\mathcal{T}_L^{(i)}\) which consist of the restriction of \(\mathcal{T}_L\) to the elements in \(\Omega^{(i)}\) and to the vertices in \(\overline{\Omega^{(i)}}\). For convenience in introducing the quantum framework, we next introduce a simple geometry (but we stress that this is not necessary for domain decomposition). We first introduce the subdomains \(\Omega^{(i)}\) and then set \(\Omega\) to be their union.

We assume that each subdomain \(\Omega^{(i)}\) is a unit hypercube discretized by a Cartesian mesh \(\mathcal{T}_L^{(i)}\) of size \(2^{-L}\). Consequently, each subdomain contains \(2^{dL}\) elements (\(2^L\) in each coordinate direction). We further assume that there are \(N_s\) subdomains along each direction \(s= 1, \dots, d\), and that they overlap by a fixed number of mesh elements given by \(2 \delta / 2^{-L}\). (This ensures that \(\delta\) matches the standard overlap parameter used in the domain decomposition literature.) Figure 1 shows an example of this geometry for \(d=2\), \(N_1 = 4\), \(N_2 = 2\), \(L = 3\), and \(\delta = 2^{-L}\).

The global domain is then \[\Omega = \prod_{s=1}^d \left[ 0, N_s - 2\delta (N_s - 1) \right].\] The global mesh \(\mathcal{T}_L\) is the Cartesian mesh of \(\Omega\) with discretization step \(2^{-L}\). The finite element space \(V_L\) (already defined in Section 2.2) has \(\mathcal{N}\) degrees of freedom (i.e., unknowns) at the interior mesh vertices with \(\mathcal{N} = \prod_{s=1}^d \left[2^L (N_s - 2 \delta (N_s-1)) -1\right]\) . This is also the order of the matrix \(A\) in the linear system 4 .

We next introduce in each subdomain a local finite element space \[V_L^{(i)} := \left\{u_L \in C^0(\overline{\Omega^{(i)}}); \, u_L = 0 \text{ on } \partial \Omega \text{ and } {u_L}_{|K} \in \mathbb{Q}_1(K), \, \forall K \in \mathcal{T}_L^{(i)} \right\}.\] The dimension of each \(V_L^{(i)}\) is \((2^L-1)^{d}\) (the number of inner vertices of \(\Omega^{(i)}\)). There is a natural embedding of \(V_L^{(i)}\) into \(V_L\) since the functions in \(V_L^{(i)}\) can be seen as restrictions of certain functions in \(V_L\). To make things concrete, we also introduce this interpolation operator in the algebraic formalism. The degrees of freedom (interior vertices) of \(\Omega\) are numbered from \(1\) to \(\mathcal{N}\). The degrees of freedom (interior vertices) of \(\Omega^{(i)}\) are a subset of vertices that we can denote by \(\omega_i\). The restriction matrix \(R_L^{(i)} \in \mathbb{R}^{(2^L - 1)^d \times \mathcal{N}}\) is \[R_L^{(i)} := I_\mathcal{N}(\omega_i, :),\] by which it is meant that we select all columns and only rows that are in \(\omega_i\) . Matrix \(R_L^{(i)}\) is the representation of the restriction operator \(V_L \rightarrow V_L^{(i)}\). Its transpose \({R_L^{(i)}}^\top\) represents the prolongation by zero operator \(V_L^{(i)} \rightarrow V_L\), which we also call the natural embedding \(V_L^{(i)} \hookrightarrow V_L\).

Example: In the one-dimensional case with 2 subdomains, \(L = 2\), and overlap \(\delta=1/4\), the extension matrices are \[{R_L^{(1)}}^\top = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix} \text{ and } {R_L^{(2)}}^\top = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{pmatrix}.\]

The second last ingredient for defining domain decomposition preconditioners is to define, for each \(i\), a local solver that acts in the local subspaces \(V_L^{(i)}\). In a first step, we choose the exact local solver \({A^{(i)}}^{-1}\) with \[A^{(i)} := R_L^{(i)} A {R_L^{(i)}}^\top = A(\omega_i, \omega_i),\] by which it is meant that we select the rows and columns of \(A\) that are in \(\omega_i\). Later on, these \(A^{(i)}\) will be replaced by approximations.

It is to be noted that \(A^{(i)}\) is also the matrix associated with the finite element discretization, in \(V_L^{(i)}\), of the PDE \[\left\{ \begin{align} - \nabla \cdot \rho \nabla u &= f \quad \text{in } \Omega^{(i)}, \\ u &= 0 \quad \text{on } \partial \Omega^{(i)}. \end{align} \right.\]

3.2 Domain Decomposition preconditioners↩︎

Having introduced the interpolation operators (\({R_L^{(i)}}\), \({R_L^{(i)}}^\top\)) between local and global functions, as well as the local solvers \({A^{(i)}}^{-1}\) we can next define the one level Additive Schwarz preconditioner as \[H_1:= \sum_{i=1}^N {R_L^{(i)}}^\top {A^{(i)}}^{-1} R_L^{(i)}.\] To say it with words, the inverse of \(A\) is approximated by a sum of local inverses \({A^{(i)}}^{-1}\) in the subdomains. This approximation is, in general, not sufficient to guarantee that the eigenvalues of the preconditioned operator \(H_1A\) are bounded independently from the number of subdomains (what is called the scalability property). For this reason, an additional correction is performed called the coarse correction. The coarse correction solves one more problem in a space called the coarse space (which is a subset of \(V_L\)). We choose to parametrize the coarse space by a matrix \(Z\), the columns of which form a basis of the coarse space. The coarse problem is then \(Z^\top A Z\), and the two-level Additive Schwarz preconditioner is defined by \[\label{eq:HASM2} H_2:= \sum_{i=1}^N {R_L^{(i)}}^\top {A^{(i)}}^{-1} R_L^{(i)} + Z (Z^\top A Z)^{-1} Z^\top.\tag{7}\]

The preconditioned matrix \(H_2A\) that results from preconditioning by domain decomposition has the very appealing property that its eigenvalues can be bounded independently of the mesh size (parameterized by \(L\) in our notation), and, if \(Z\) is chosen correctly, of the number of subdomains (\(N\) in our notation). This is fundamental in classical computing. The solver of choice would indeed be the preconditioned conjugate gradient (PCG) method, and it converges fast if the eigenvalues of the preconditioned matrix are clustered away from zero (see [32] and [24]). We state below a fundamental result from the domain decomposition literature.

We have injected in brackets the formula that is particular to our choice of geometry.

The local solves \({A^{(i)}}^{-1}\) in the definition of the domain decomposition preconditioner can be replaced by approximations while retaining spectral bounds.

We note that we could also consider an inexact coarse solve \((Z^\top A Z)^{-1}\) with no extra difficulty.

All above results are for the spectrum of the preconditioned operator. Since the preconditioned matrix is non-symmetric, these do not imply bounds for the condition number of the preconditioned matrix³. In order to obtain condition number bounds, we use again an idea from [18] and split the preconditioner to symmetrize the preconditioned operator.

3.3 Split form of the preconditioned operator↩︎

Let’s assume that, for any \(i=1,\dots,N\), the inexact local solvers \(H^{(i)}\) can be factorized as \(H^{(i)} =: F^{(i)} {F^{(i)}}^\top\) and that the coarse solve can be factorized as \((Z^\top A Z)^{-1} =: F^{(0)} {F^{(0)}}^\top\). We do not require the matrices \(F^{(i)}\) and \(F^{(0)}\) to be square. Then \(\widetilde{H_2}\) rewrites \[\widetilde{H_2}= \sum_{i=1}^N {R_L^{(i)}}^\top F^{(i)} {F^{(i)}}^\top R_L^{(i)} + Z F^{(0)} {F^{(0)}}^\top Z^\top = \tilde{F} \tilde{F}^\top, \label{eq:decomposition95preconditioner}\tag{8}\] with \[\tilde{F} := \left({R_L^{(1)}}^\top F^{(1)} | ~\dots~| {R_L^{(i)}}^\top F^{(i)} | \dots~ | {R_L^{(N)}}^\top F^{(N)} | Z F^{(0)} \right).\]

By this notation, we mean that we have concatenated all matrices of the form \({R_L^{(i)}}^\top F^{(i)}\) horizontally with the matrix \(Z F^{(0)}\). The matrix \(\tilde{F}\) is a rectangular matrix. It has \(\mathcal{N}\) rows.

Based on the decomposition of the stiffness matrix 6 and the decomposition of the preconditioner 8 , we arrive at a symmetric preconditioned system \[\tilde{F}^\top C_L^\top (D_{\rho} \otimes Id_{2^d}) C_L \tilde{F} y = \tilde{F}^\top b , \text{ with } x=\tilde{F} y. \label{eq:prec95system}\tag{9}\] This is the system that is solved in fine by our methodology. The symmetric preconditioned operator is \(\tilde{F}^\top C_L^\top (D_{\rho} \otimes Id_{2^d}) C_L \tilde{F}\). Its order is larger than the order \(\mathcal{N}\) of the original preconditioned matrix \(\widetilde{H_2}A\). Its rank is \(\mathcal{N}\). The non-zero eigenvalues of \(\tilde{F}^\top C_L^\top (D_{\rho} \otimes Id_{2^d}) C_L \tilde{F}\) are the eigenvalues of \(\widetilde{H_2}A\) so that \[\kappa(\tilde{F}^\top C_L^\top (D_{\rho} \otimes Id_{2^d}) C_L \tilde{F}) = \frac{\lambda_{\max}(\widetilde{H_2}A)}{\lambda_{\min}(\widetilde{H_2}A)},\] where \(\lambda_{\max}(\widetilde{H_2}A)\) and \(\lambda_{\min}(\widetilde{H_2}A)\) are the extreme eigenvalues of \(\widetilde{H_2}A\). These can be bounded by Theorems 1 and 2. The solution \(y\) to 9 is defined only up to an element in the kernel of \(\tilde{F}\) but this does not modify the corresponding solution \(x = \tilde{F} y\) of the original system.

In Section 4, a block-encoding of the matrix \(\tilde{F}^\top C_L^\top (D_{\rho} \otimes Id_{2^d}) C_L \tilde{F}\) is proposed. Because the matrices \(C_L\) and \(\tilde{F}\) are not well-conditioned it is chosen to avoid implementing it as the product of block-encodings. Instead, the product \[C_L \tilde{F} = \left(C_L {R_L^{(1)}}^\top F^{(1)} | ~\dots~| C_L{R_L^{(i)}}^\top F^{(i)} | \dots~ | C_L{R_L^{(N)}}^\top F^{(N)} | C_L Z F^{(0)} \right)\] is reformulated so that each of its blocks is block-encoded at once.

For each \(i = 1, \dots, N\) the local discontinuous Galerkin space is introduced as \[Q_L^{(i)} := \left\{u_L^{(i)} \in L^2(\Omega^{(i)});\, {u_L^{(i)}}_{|K} \in \mathbb{Q}_1(K), \;\forall K \in \mathcal{T}_L^{(i)} \right\},\] i.e., the local counterpart of \(Q_L\) from equation 5 . The vector-valued space \({Q_L^{(i)}}^d := Q_L^{(i)} \otimes \mathbb{R}^d\) is also introduced. These spaces are equipped, respectively, with with an \(L^2\)-orthonormal basis consisting of \(2^d\) functions \(\psi_{K,1}, \dots, \psi_{K,2^d}\) on each element \(K \in \mathcal{T}_L^{(i)}\), and with the basis functions \(\psi_{K,k} \otimes e_s\) where \(e_s\) is the \(s\)-th canonical basis vector in \(\mathbb{R}^{d}\). We also let \(C_L^{(i)}\) denote the matrix representation of the (elementwise) gradient operator \[\nabla : V_L^{(i)} \to {Q_L^{(i)}}^d.\] There is a natural embedding of \({Q_L^{(i)}}^d\) into \(Q_L^{d}\). It induces an extension operator from \({Q_L^{(i)}}^d\) into \(Q_L^{d}\) that is represented by a matrix which we denote by \({\mathcal{R}_L^{(i)}}^\top\). Each term \(C_L {R_L^{(i)}}^\top\) can then be rewritten as \[C_L {R_L^{(i)}}^\top = {\mathcal{R}_L^{(i)}}^\top C_L^{(i)}, \label{eq:switch95gradient}\tag{10}\] and both are equal to \(C_L(:,\omega_i)\), i.e., the selection of columns in \(C_L\) that corresponds to degrees of freedom of \(V_L\) that are in the interior of \(\Omega^{(i)}\). This can be verified by considering the operators represented by these matrices: for a function in \(V_L^{(i)}\), it is equivalent to first prolongate it by zero to \(\Omega\) and then compute the element-wise (discontinuous) gradient, or to first compute the element-wise (discontinuous) gradient and then prolongate it by zero to \(\Omega\) as summarized in the diagram below \[\begin{tikzcd}[row sep=normal, column sep=huge] V_L^{(i)} \arrow[r,hook, "{ {R_L^{(i)}}^\top }"] \arrow[d, "C_L^{(i)}"'] & V_L \arrow[d, "C_L"] \\ {Q_L^{(i)}}^d \arrow[r,hook, "{ {\mathcal{R}_L^{(i)}}^\top }"'] & Q_L^d \end{tikzcd}\]

We arrive at \[C_L \tilde{F} = \left( {\mathcal{R}_L^{(1)}}^\top C_L^{(1)} F^{(1)} |~\dots~| {\mathcal{R}_L^{(i)}}^\top C_L^{(i)} F^{(i)} | \dots~| {\mathcal{R}_L^{(N)}}^\top C_L^{(N)} F^{(N)} |C_L Z F^{(0)} \right), \label{eq:precond95gradient}\tag{11}\] and propose to block-encode this matrix in the next section, as a step towards block-encoding \(\tilde{F}^\top C_L^\top (D_{\rho} \otimes Id_{2^d}) C_L \tilde{F}\) and proposing a QLSS for 9 . We have not transformed the term \(C_L Z F^{(0)}\).

3.4 Advantages of domain decomposition↩︎

Now that we have introduced domain decomposition, also known as DD, we present in more detail some of its strengths that could carry over to the quantum formalism.

• As shown in Theorem 1, preconditioning a linear system by domain decomposition allows to improve significantly the eigenvalue distribution and to make it independent of some parameters such as the mesh size. Theorem 1 is only one example of such a result and we want to stress that domain decomposition preconditioners have been developed for a very wide range of PDEs. For this reason, the feasibility of quantum domain decomposition that is achieved in this article paves the way to solving many more problems than the Poisson problem on a hyperrectangle.

• Another advantage of domain decomposition is that it can be viewed as an abstract framework into which various local and coarse solvers can be plugged. This is by now quite standard in classical computing (in particular via, the PETSc library [34]) and such versatility could be carried over to quantum implementations.

• One particular strength of domain decomposition is to address the solution of PDEs with heterogeneous coefficients. The level of difficulty depends on whether the heterogeneities are across the subdomains boundaries, inside the subdomains or whether they cross the boundaries. A good choice of coarse space is instrumental to ensure convergence with respect to the coefficient distribution [24], [26]–[30].

• The coarse space formalism mentioned above may turn out not to be the best suited for quantum implementation. However, it remains from the corresponding analysis that only a small number of eigenvalues are isolated and responsible for an unfavourable condition number. These could be filtered away by other techniques.

• Domain decomposition is also a suitable framework for addressing multiphysics formulations.

• Although it is not the point of view adopted in this article, domain decomposition algorithms could be useful to parallelize computations over multiple quantum processing units.

• A recent contribution by [31] is to alleviate Barren Plateau-like effects in variational quantum algorithms by replacing a global cost function by localized ones by domain decomposition techniques.

• A final group of ideas concerns the reformulation of \(Ax=b\) through substructuring [24]. After partitioning the degrees of freedom into interior and overlapped variables \(x = (x_I, x_\Gamma)\), the interior unknowns can be eliminated via the local subdomain operators \(A^{(i)}\). The global system can then be reduced to an equivalent Schur complement problem posed only on the interface variables: \[S x_\Gamma = b_\Gamma.\] In this formulation, the unknown that must be represented or accessed by a quantum algorithm is no longer the full vector \(x\), but only its interface component \(x_\Gamma\). The interior degrees of freedom can be reconstructed classically through independent local subdomain solves, so that the quantum state preparation is restricted to a reduced interface space.

• Elimination of the interior degrees of freedom also changes the structure of the data that must be encoded in quantum routines for linear system solving. Indeed, after static condensation, the effective right-hand side for the reduced system is supported only on the overlapping degrees of freedom. As a consequence, both the state preparation for \(b\) and the representation of residuals in iterative schemes can be restricted to a smaller overlap-supported set of components. This suggests a reduction in the number of amplitudes required to encode problem data for quantum linear algebra subroutines.

• Finally, this structure could be exploited within the hybrid QLSS proposed in [35]. The algorithm is an iterative refinement strategy inspired by mixed-precision numerical linear algebra techniques. A first quantum solution is computed using QSVT in low precision, and then refined in higher precision until a satisfactory accuracy is achieved. Since the residual driving the refinement inherits the overlap-supported structure induced by domain decomposition, the QSVT correction step can be reformulated to act on reduced interface residuals rather than full-system vectors. This suggests that each quantum linear system call within the iterative refinement loop may be applied to lower-dimensional data, potentially reducing quantum resources per iteration while preserving the convergence properties of the original scheme.

3.5 Numerical illustration of domain decomposition↩︎

The purpose of this section is to illustrate the domain decomposition knowledge that we have introduced, and to test the proposed preconditioner before carrying it over to the quantum formalism. We consider the two-dimensional Poisson problem (\(d=2\)) with \(f=1\) and \(\rho=I_d\) over the whole domain \(\Omega\). The reported results are for the ratio between the largest and smallest eigenvalue of the (un)preconditioned operator. This is also the condition number of the split preconditioned operator so when we speak of condition number in this section we refer to \({\lambda_{\max}}/{\lambda_{\min}}\).

In 1, the number of subdomains varies in both directions. The mesh size is \(h=2^{-4}\) (i.e. L = 4), and the overlap is \(\delta = 2^{-4}\). Two-level Additive Schwarz preconditioners \(H_2\) and \(\widetilde{H_2}\) with exact and inexact local solves are compared to two-level hybrid Schwarz preconditioners \(H_{\text{hyb}}\) and \(\widetilde{H_{\text{hyb}}}\) with exact and inexact local solves. For the inexact local solves, we replace \({A^{(i)}}^{-1}\) by the standard BPX preconditioner [19]. Its definition is deferred to 4.4. We first observe that all preconditioners drastically improve the condition number: from \(1499\) to at most \(76\) when there are \(6 \times 6\) subdomains. In agreement with Theorem 1, the condition numbers of the preconditioned operators do not depend significantly on the number of subdomains. This is particularly true for the hybrid two-level preconditioner \(H_{\text{hyb}}\) with an almost constant condition number close to \(25\). For \(H_2\), the increase in condition number does not contradict the theory, it indicates that the asymptotic regime has not been reached. The preconditioners with inexact local solves perform only slightly worse than their exact counterparts which justifies their use.

We consider a banded domain for 2 and 3: subdomains are only added along direction \(1\) (i.e. \(N_2=1\)). In this case, a classical result in domain decomposition is that the coarse correction is not necessary because all subdomains are connected through the Dirichlet boundary condition. 2 illustrates the influence of the mesh size for the one-level Additive Schwarz preconditioners \(H_1\) (exact) and \(\tilde{H}_1\) (inexact). The overlap is kept constant at \(\delta = 2^{-3}\) and the number of subdomains is \(N_1=8\). Using exact local solves, the eigenvalue ratio does not depend on the mesh size. This is in agreement with Theorem 1. In the case of inexact local solves, the eigenvalue ratio slightly increases from \(4.2\) to \(8.5\) when the mesh size decreases from \(2^{-3}\) to \(2^{-5}\).

In 3, we study the influence of the overlap \(\delta\) for the one-level Additive Schwarz preconditioners. The number of subdomains and the mesh size are kept constant (\(N_1=8\) and \(h=2^{-5}\) corresponding to \(L=5\)). For both the exact and inexact local solves, the eigenvalue ratio decreases when the overlap increases. Once more, this corresponds to the prediction of Theorem 1.

4 Quantum algorithm↩︎

We first introduce the block-encoding framework and state the complexity of a particular quantum linear system solver. Then, we discuss the block-encoding of the symmetric preconditioned matrix \(\tilde{F}^\top A \tilde{F}\), which was introduced in the previous section by symmetrizing \(\widetilde{H_2}A\). We also discuss the implementation of the block-encoding and the cost of inverting the linear system (in the quantum sense).

4.1 Preliminaries↩︎

The block-encoding framework allows to embed a general matrix inside a larger unitary matrix [6]. We use a general definition [6], [18] also known as projected unitary encoding.

Considering a general, well-posed, system of the form \(\mathcal{A}x = b\), quantum linear system solvers output a quantum state close to \(| x \rangle \propto \mathcal{A}^{+} | b \rangle\). We recall the complexity of a QSVT-based QLSS from [17] and [36].

The complexity can be further improved to almost linear in \(\kappa(\mathcal{A})\) using variable time amplitude amplification [2]. In the cited references, \(\mathcal{A}\) is non-singular and the lemma is written with \(\mathcal{A}^{-1}\) instead of the pseudo-inverse \(\mathcal{A}^+\). However pseudo-inversion is not a difficulty for QSVT and was already considered in the original article [6].

From this result, we observe that QSVT is efficient if \(\mathcal{A}\) is well conditioned, and if there exists a block-encoding of \(\mathcal{A}\) with a small subnormalization factor. To improve the condition number, we precondition by domain decomposition as proposed in Section 3. It is then the condition number of the split preconditioned operator \[\tilde{F}^\top A \tilde{F} = \tilde{F}^\top C_L^\top (D_{\rho} \otimes Id_{2^d}) C_L \tilde{F}.\] which plays a role and we have already analyzed it. We are left with the analysis of the subnormalization factor \(\tilde{\alpha}(U_{\tilde{F}^\top A \tilde{F}})\). We will propose a particular choice of local solvers \({F^{(i)}}^\top F^{(i)}\) in the domain decomposition preconditioner for which \(\tilde{\alpha}(U_{\tilde{F}^\top A \tilde{F}})\) remains bounded independently of the problem parameters.

4.2 Block-encoding of the preconditioned system↩︎

The following theorem presents a block-encoding of the preconditioned system provided access to a block-encoding of the coefficient \(D_\rho\) (following e.g, [18]), as well as to block-encodings of the local matrices \(C_L^{(i)} F^{(i)}\) and to a block encoding of the coarse correction \(C_L Z F^{(0)}\). We recall that these matrices form the blocks of \(C_L \tilde{F}\) in 11 , the definition of which we state again for clarity: \[C_L \tilde{F} = \left( {\mathcal{R}_L^{(1)}}^\top C_L^{(1)} F^{(1)} |~\dots~| {\mathcal{R}_L^{(i)}}^\top C_L^{(i)} F^{(i)} | \dots~| {\mathcal{R}_L^{(N)}}^\top C_L^{(N)} F^{(N)} | C_L Z F^{(0)} \right).\]

In order to propose a full quantum domain decomposition result, we discuss implementation choices in Section 4.3, particularly for the block encoding of the matrices \(\mathcal{R}_L^{(i)}\). Then, in Section 4.4, we set the local solvers \(F^{(i)}{F^{(i)}}^\top\) to be the BPX approximation of \({A^{(i)}}^{-1}\), and give our final result in 5.

4.3 Implementation↩︎

Up until now, we have not considered the concrete implementation of the block-encodings. An efficient circuit that block-encodes the preconditioned system \(U_{\tilde{F}^\top A \tilde{F}}\) is necessary in order to ensure the efficiency of the QLSS. The same remark applies to the state-preparation oracles.

We first detail the different registers that are involved in the block-encodings. There are five:

• Registers \(| j \rangle\), \(| k \rangle\) and \(| s \rangle\) are the same as in [18]. Register \(| j \rangle\) encodes the mesh elements or the degrees of freedom of the local Lagrange finite element spaces. Based on the tensor product structure of the Cartesian mesh and of the Lagrange spaces, it consists of \(d\) sub-registers \(| j_1 \rangle \cdots | j_d \rangle\), each containing \(L\) qubits. Register \(| k \rangle\) is associated with the \(2^d\) degrees of freedom per mesh element of the local Discontinuous Galerkin spaces, so it contains \(d\) qubits \(| k_1 \rangle \cdots | k_d \rangle\). Register \(| s \rangle\), which consists of \(\log_2(d)\) qubits, encodes the \(d\) components of the gradient.

• We introduce a register \(| i \rangle\) associated with the \(N\) subdomains. It is composed of \(d\) registers \(| i_1 \rangle \cdots | i_d \rangle\) (see Remark [rem:i-multiindex]). For \(1 \leq s \leq d\), register \(| i_s \rangle\) contains \(n_s= log_2(N_s)\) qubits (where \(N_s\) is the number of subdomains in direction \(s\)).

• Finally, we introduce a new register \(| \tilde{i} \rangle\) which has the same structure as \(| i \rangle\) in order to perform quantum domain decomposition.

4.3.0.1 Local to global Prolongation operators

We now propose a block-encoding of the prolongation operators \({\mathcal{R}^{(i)}_L}^\top : {Q_L^{(i)}}^d \hookrightarrow Q_L^{d}\) defined in Section 3.3. Each matrix \({\mathcal{R}^{(i)}_L}^\top\) is of order \(d 2^{(L+1)d} \prod_{s=1}^d (N_s - 2 \delta (N_s-1)) \times d 2^{(L+1)d}\). Having split the domain \(\Omega\) into an overlapping Cartesian grid of subdomains, the prolongation operators display a tensor-product structure, \[{\mathcal{R}^{(i)}_L}^\top = I_d \otimes {\mathcal{R}^{(i_1)}_{L,1D}}^\top \otimes \cdots \otimes {\mathcal{R}^{(i_d)}_{L,1D}}^\top,\] where each \(\mathcal{R}^{(\text{ind}_s)}_{L,1D}\) performs the restriction to local elements for the \(s\)-th coordinate. Therefore, we can focus on the implementation of one of the \({\mathcal{R}^{(\text{ind}_s)}_{L,1D}}^\top\) with \(1 \leq \text{ind}_s \leq N_s\). The index of a Discontinuous Galerkin degree of freedom encoded by \(| k_s \rangle | j_s \rangle | i_s \rangle\) is \(k_s + 2 j_s +2^{L+1} i_s - 2^{L+2} \delta i_s\). The term \(-2^{L+2} \delta i_s\) accounts for the overlap. Quantum arithmetic based on the Fourier transform [37] enables the implementation of the operator \(U_{{\mathcal{R}^{(\text{ind}_s)}_{L,1D}}^\top}\) \[| k_s,j_s,i_s \rangle | \tilde{i}_s \rangle \rightarrow | k_s + 2j_s +2^{L+1}i_s - 2^{L+2}\tilde{i}_s \delta \rangle | \tilde{i}_s \rangle.\]

The row index encoded by \(| k_s,j_s,i_s,\tilde{i}_s \rangle\) is the integer \(p:=k_s + 2j_s +2^{L+1}i_s + 2^{n_s+L+1}\tilde{i}_s\). The \(C_{\Pi_{row}} \text{NOT}\) gate flips the ancilla qubit when the index of the row satisfies \(p \leq M_s := N_s 2^{L+1} - 2^{L+2} \delta(N_s -1) + N_s 2^{L+1} (\text{ind}_s - 1)\) and \(p \geq N_s 2^{L+1} (\text{ind}_s - 1)\). The number of rows selected by \(C_{\Pi_{row}} \text{NOT}\) corresponds to the total number of degrees of freedom for the Discontinuous Galerkin space in direction \(s\), i.e. \(N_s 2^{L+1} - 2^{L+2} \delta(N_s -1)\). Each of the two inequality tests can also be performed by QFT-based quantum arithmetic with one ancilla qubit [38], [39]. For the inequality \(p \leq M_s\), the following operations are carried out in the Fourier spaces \[\begin{align} | k_s,j_s,i_s,\tilde{i}_s \rangle | 0 \rangle &\rightarrow | p - M_s - 1 \rangle_{p \geq M_s+1} | 0 \rangle + | p - M_s -1 + N_s^2 2^{L+1} \rangle_{p < M_s+1} | 1 \rangle \\ & \rightarrow | k_s,j_s,i_s,\tilde{i}_s \rangle | p < M_s+1 \rangle. \end{align}\] Similar operations are performed for the inequality \(p \geq N_s 2^{L+1} (\text{ind}_s - 1)\). The \(C_{\Pi_{col}}\text{NOT}\) gate flips the ancilla qubit when \(| \tilde{i}_s \rangle = | i_s \rangle = \text{ind}_s-1\). This equality test can be performed using quantum arithmetic based on multi-controlled gates. This gate selects the \(2^{L+1}\) columns.

4.3.0.2 Preparation of the right hand side

The efficiency of the quantum linear system solver also relies on the efficiency of the state preparation. We address the preparation of the preconditioned right hand side \[{\tilde{F}}^\top b = \begin{pmatrix} {F^{(1)}}^\top R_L^{(1)} b \\ \vdots \\ {F^{(N)}}^\top R_L^{(N)} b \\ {F^{(0)}}^\top Z^\top b \end{pmatrix}.\]

We assume that we already know how to efficiently encode local components \(\displaystyle | {F^{(i)}}^\top R_L^{(i)} b \rangle := \frac{{F^{(i)}}^\top R_L^{(i)} b }{\|{F^{(i)}}^\top R_L^{(i)} b \|}\) and the coarse component \(\displaystyle | {F^{(0)}}^\top Z^\top b \rangle := \frac{{F^{(0)}}^\top Z^\top b }{\|{F^{(0)}}^\top Z^\top b \|}\) using \(m\) qubits. Supposing access to all of the norms, we can then prepare the \((n+1)\)-qubit state \(| 0 \rangle_{n+1} \rightarrow \frac{\|{F^{(0)}}^\top Z^\top b\|}{\|{\tilde{F}}^\top b\|} | 0 \rangle_{n+1} + \sum_{i=1}^N \frac{\|{F^{(i)}}^\top R_L^{(i)} b\|}{\|{\tilde{F}}^\top b\|} | i \rangle_{n+1}\) with a cost at most linear in \(N=2^n\). Controlling the state preparation of \(| {F^{(i)}}^\top R_L^{(i)} b \rangle\) and \(| {F^{(0)}}^\top Z^\top b \rangle\) on the \(n+1\)-qubit register, we get \(| \tilde{F}^\top b \rangle\). More specifically, these two steps lead to

\[\begin{align} | 0 \rangle_{n+1} | 0 \rangle_m \rightarrow & \left( \frac{\|{F^{(0)}}^\top Z^\top b\|}{\|{\tilde{F}}^\top b\|} | 0 \rangle_{n+1} + \sum_{i=1}^N \frac{\|{F^{(i)}}^\top R_L^{(i)} b\|}{\|{\tilde{F}}^\top b\|} | i \rangle_{{n+1}} \right) | 0 \rangle_{m}, \\ \rightarrow & \frac{\|{F^{(0)}}^\top Z^\top b\|}{\|{\tilde{F}}^\top b\|} | 0 \rangle | {F^{(0)}}^\top Z^\top b \rangle + \sum_{i=1}^N \frac{\|{F^{(i)}}^\top R_L^{(i)} b\|}{\|{\tilde{F}}^\top b\|} | i \rangle | {F^{(i)}}^\top R_L^{(i)} b \rangle = | \tilde{F}^\top b \rangle. \end{align}\]

4.4 Domain decomposition with local BPX solvers↩︎

One possible choice of local solver is the Bramble–Pasciak–Xu (BPX) preconditioner [19] which is spectrally equivalent to \(A^{(i)}\), and for which [18] proposes a quantum implementation. We first define the BPX preconditioner in subdomain \(\Omega^{(i)}\).

Remark that \(F^{(i)} \in\mathbb{R}^{\text{dim } V_L^{(i)} \times \sum_{l=1}^L \text{dim } V_l^{(i)} }\). We next state the spectral equivalence result.

Another advantage of choosing BPX as our local solver is that, the block-encodings \(U_{C_L^{(i)} F^{(i)}}\) in the Assumptions of 3 are exactly the block-encodings proposed in [18] so we can write the following result.

We are interested in solving the linear system 9 and computing some quantities of interest.

We remark that the complexity scales linearly with the number of subdomains. This is due to the particular construction of the block-encoding of the preconditioned system. In domain decomposition, it is standard to eliminate \(N\) and instead have a dependency in the number of colors \(N_C\) needed to color the subdomains so that two overlapping subdomains never share the same color (see e.g, the spectral result in Theorem 1). We leave this for future work in the quantum formalism.

5 Numerical simulations↩︎

We now turn to numerical simulations for two-level Additive Schwarz with BPX local solves. We have not included a coarse space in these simulations. The quantum circuits are emulated in Qiskit [41] for a one-dimensional domain (\(d=1\)). The code associated with [18] is used for BPX local solvers. The inner product is split as \(\left[(\tilde{F}^\top C_L^\top)^+ \tilde{w} \right]^\top (\tilde{F}^\top C_L^\top)^+ \tilde{b}\) [18]. The linear systems are inverted using QSVT with the polynomial from [4] and [6]. This polynomial is normalized to satisfy the requirements of QSVT. Its degree, denoted by \(D\), is computed using the quantity \(\kappa(\tilde{F}^\top C_L^\top) \tilde{\alpha}(\tilde{F}^\top C_L^\top)\) and the target QSVT accuracy \(\epsilon=10^{-6}\). The angles for the polynomials are computed using the library pysqp [42]–[44]. The domain \(\Omega\) depends on the number of subdomains. We choose the right-hand side as \(f=1\) so that \(b=M (1 \cdots 1)^\top\) where \(M\) is the mass matrix. We choose \(w=b\) for the computation of the inner product. In that case, only the norm of the quantum solution is estimated from the probability of success of QSVT. The norm is multiplied by the normalization constant of the polynomial and by the norm of the right-hand side and divided by the normalization constant. The mesh size is \(2^{-L}\) with \(L=4\). Two values of numbers of subdomains are considered \(N_1 = 2\) and \(N_1=4\). The simulations are run \(100\) times, each with \(10^6\) samples. The results are reported in 4. We check that the computed mean value \(\mathcal{Q}_{\text{mean}}\) is a good approximation of \(\mathcal{Q}_{\text{ref}}=w^\top b\) computed classically. We also report the effective condition number \(\kappa(\tilde{F}^\top C_L^\top) \tilde{\alpha}(U_{\tilde{F}^\top C_L^\top})\) and the degree of the inversion polynomial. Both depend on the number of subdomains \(N_1\). Adding the coarse correction would eliminate this dependency.

6 Conclusions↩︎

In this work, we have introduced the first quantum domain decomposition preconditioner. We have demonstrated that it is feasible to block-encode the preconditioned problem that results from preconditioning Poisson by the Additive Schwarz preconditioner. We have studied the overall complexity of measuring a quantity of interest. In future work we plan to investigate other choices for the local solves: sparse approximate inverses, circulant preconditioners, Neumann solves…Particular attention should also be paid to the quantum coarse correction. Finally, we are optimistic that the linear dependency of the complexity in the number of subdomains can be replaced by the coloring constant.

References↩︎

---

1. Thales CortAIx-Labs, Palaiseau, France ( ).↩︎

2. CNRS, CMAP, École polytechnique, Institut Polytechnique de Paris, Palaiseau, France ().↩︎

3. Note that in the DD literature, it is often said that the condition number is bounded. However this refers to the condition number for the norm induced by \(A\), not the spectral condition number.↩︎

4. In this context a projection is a matrix with orthonormal rows.↩︎