Efficiently decomposing finite element stiffness matrices into the Pauli basis is challenging due to the exponential growth of Pauli strings with problem size. A naive Pauli expansion requires $Θ(8^{\lceil \log_2 N \rceil})$ operations, where $N$ denotes the number of degrees of freedom, rendering direct decomposition infeasible for large systems. Existing approaches exploit algebraic sparsity or operator structure but do not incorporate the geometric organization intrinsic to finite element discretizations, and consequently exhibit poor scaling for stiffness matrices. To address this problem, we introduce PHASE, a hierarchical, geometry-aware Pauli decomposition algorithm that leverages recursive mesh partitioning to organize element contributions across multiple spatial scales. PHASE employs a hybrid strategy that combines full- and reduced-space Tensorized Pauli Decomposition with Fast Walsh-Hadamard Transform-based aggregation to assemble global Pauli coefficients efficiently. We show that this approach yields a dimension-dependent reduction in the exponential scaling exponent of Pauli assembly asymptotic complexity relative to existing methods, reducing the cost from $2^{2{\lceil \log_2 N \rceil}}$ to $2^{γ_d{\lceil \log_2 N \rceil}}$ with $γ_d < 2$ under standard mesh regularity and balanced partition assumptions. These results substantially improve the feasibility of quantum-compatible operator synthesis for large-scale finite element models.
Efficiently decomposing finite element stiffness matrices into the Pauli basis is challenging due to the exponential growth of Pauli strings with problem size. A naive Pauli expansion requires Θ(8⌈log2N⌉) operations, where N denotes the number of degrees of freedom, rendering direct decomposition infeasible for large systems. Existing approaches exploit algebraic sparsity or operator structure but do not incorporate the geometric organization intrinsic to finite element discretizations, and consequently exhibit poor scaling for stiffness matrices. To address this problem, we introduce PHASE, a hierarchical, geometry-aware Pauli decomposition algorithm that leverages recursive mesh partitioning to organize element contributions across multiple spatial scales. PHASE employs a hybrid strategy that combines full- and reduced-space Tensorized Pauli Decomposition with Fast Walsh-Hadamard Transform-based aggregation to assemble global Pauli coefficients efficiently. We show that this approach yields a dimension-dependent reduction in the exponential scaling exponent of Pauli assembly asymptotic complexity relative to existing methods, reducing the cost from 22⌈log2N⌉ to 2γd⌈log2N⌉ with γd<2 under standard mesh regularity and balanced partition assumptions. These results substantially improve the feasibility of quantum-compatible operator synthesis for large-scale finite element models.
Decomposing classical matrices into linear combinations of Pauli strings is a major bottleneck for end-to-end implementations of near-term quantum algorithms. In this work, we consider a promise version of this Pauli decomposition problem in which the matrix is guaranteed to have support on only k=poly(n) Pauli strings and is given through classical sparse query access. Existing Pauli decomposition algorithms are designed for the generic, dense problem and do not inherently take advantage of this promised sparsity, so these approaches take time that is exponential in n. We present a randomized classical algorithm that does take advantage of this sparsity and recovers the exact Pauli decomposition with success probability at least 1−δ, for any δ. Under the stated access model, the algorithm executes with query and runtime complexity that is polynomial in n, k, and log(1/δ). These results show that, even though finding the Pauli decomposition is exponentially hard for general matrices, it becomes efficiently solvable for matrices that are known to be sparse in the Pauli basis, a regime that is relevant to near-term quantum algorithms operating on structured classical input.
Decomposing a matrix into a weighted sum of Pauli strings is a common chore of the quantum computer scientist, whom is not easily discouraged by exponential scaling. But beware, a naive decomposition can be cubically more expensive than necessary! In this manuscript, we derive a fixed-memory, branchless algorithm to compute the inner product between a 2^N-by-2^N complex matrix and an N-term Pauli tensor in O(2^N) time, by leveraging the Gray code. Our scheme permits the embarrassingly parallel decomposition of a matrix into a weighted sum of Pauli strings in O(8^N) time. We implement our algorithm in Python, hosted open-source on Github, and benchmark against a recent state-of-the-art method called the "PauliComposer" which has an exponentially growing memory overhead, achieving speedups in the range of 1.5x to 5x for N < 8. Note that our scheme does not leverage sparsity, diagonality, Hermitivity or other properties of the input matrix which might otherwise enable optimised treatment in other methods. As such, our algorithm is well-suited to decomposition of dense, arbitrary, complex matrices which are expected dense in the Pauli basis, or for which the decomposed Pauli tensors are a priori unknown.
Abstract The decomposition of a square matrix into a sum of Pauli strings is a classical pre-processing step required to realize many quantum algorithms. Such a decomposition requires significant computational resources for large matrices. We present an exact and explicit formula for the Pauli string coefficients which inspires an efficient algorithm to compute them. More specifically, we show that up to a permutation of the matrix elements, the decomposition coefficients are related to the original matrix by a multiplication of a generalised Hadamard matrix. This allows one to use the Fast Walsh-Hadamard transform and calculate all Pauli decomposition coefficients in O ( N 2 log N ) time and using O ( 1 ) additional memory, for an N × N matrix. A numerical implementation of our equation outperforms currently available solutions.
Pauli decomposition of an n-qubit operator is expensive if each of the 4n Pauli strings is treated as an unstructured object. Labeling those strings by the configuration space (Z2)n makes them displacement operators X⊗xZ⊗z on the computational space: shifts by~x and multiplications by characters Z⊗z in the computational basis. That labeling organizes the operator into 2n fibers of length 2n, one per shift~x, each collecting the computational-basis amplitudes that hop by~x. Expanding each fiber function in the character basis of (Z2)n --- a length-2n Walsh--Hadamard transform --- recovers the Pauli coefficients at total cost O(n⋅4n). The same fiber / shift--character picture extends to (Zp)n (Heisenberg--Weyl qudits). \texttt{paulikit} is a memory-bounded realization for qubits: for matrix-sparse inputs, peak memory need not densify to a full 2n×2n operator. Performance is checked on a stated low-end GNU/Linux platform for a coupled-oscillator sweep beyond a billion Pauli terms.
We present a novel framework for Linear Combination of Unitaries (LCU)-style decomposition tailored to structured sparse matrices, which frequently arise in the numerical solution of partial differential equations (PDEs). While LCU is a foundational primitive in both variational and fault-tolerant quantum algorithms, conventional approaches based on the Pauli basis can require a number of terms that scales quadratically with matrix size. We introduce the Sigma basis, a compact set of simple, non-unitary operators that can better capture sparsity and structure, enabling decompositions with only polylogarithmic scaling in the number of terms. We develop both numerical and semi-analytical methods for computing Sigma basis decompositions of arbitrary matrices. Given this new basis is comprised of non-unitary operators, we leverage the concept of unitary completion to design efficient quantum circuits for evaluating observables in variational quantum algorithms and for constructing block encodings in fault-tolerant quantum algorithms. We compare our method to related techniques like unitary dilation, and demonstrate its effectiveness through several PDE examples, showing exponential improvements in decomposition size while retaining circuit efficiency.
Classical numerical techniques for solving partial differential equations (PDEs) become computationally expensive as the dimension of the discretized differential operator increases. For PDEs giving rise to Sturm--Liouville problems, tensor network (TN) methods can be highly productive: an operator of dimension N×N can be represented as a matrix product operator (MPO) using only n=log2(N) qubits, enabling computation of eigenvalues and eigenvectors via imaginary time evolution (ITE). However, this remains computationally challenging. First, most methods for generating MPOs of large operators without explicit tensor-product structure require prohibitively large memory. Second, the number of Trotterization steps for convergence in conventional ITE increases rapidly with n. We present techniques to mitigate both challenges for certain sparse, structured differential operators. To address the first, we construct the MPO by expanding the operator in the Pauli-string basis, enabled by an analytical expression for the Pauli basis coefficients that reduces the memory requirement from O(2n+1) to O(2n). To address the second, we propose a multistage state-refinement heuristic that accelerates ITE convergence, reducing convergence time by up to two orders of magnitude. Using this TN framework, we compute the first 32 eigenstates of a Laplacian of dimension exceeding 106 with fidelity above 0.95 using a 20-qubit MPO. We further validate the method on the 2D anharmonic oscillator and investigate disordered systems, where increasing random potential strength degrades accuracy and limits the approach.