LRQ-Solver: Low-Rank Query PDE Solver

• LRQ-Solver is a framework that replaces globally expensive dense operations with efficient low‐rank surrogates through precomputation and hierarchical updates.
• It encompasses diverse formulations—boundary-integral updates, HODLR/BDLR multifrontal solvers, and transformer-based neural operators—to handle local geometry perturbations and large point clouds.
• The approach achieves significant speedups and memory savings by leveraging techniques like the Woodbury formula and covariance-based low-rank attention while preserving accuracy.

Low-Rank Query-based PDE Solver (LRQ-Solver) denotes, in the supplied literature, a class of low-rank, query-oriented PDE solution frameworks rather than a single invariant algorithm. In one formulation, it is a two-stage direct-solver workflow for boundary value problems on locally perturbed geometries: one precomputes a direct factorization of an original boundary-integral system and answers each new query by a low-rank update. In another, it is an HODLR- and BDLR-based solver design for 3D elliptic finite-element matrices within a multifrontal elimination process. In a third, it is the name of a transformer-based neural operator for large-scale 3D PDEs that combines Parameter Conditioned Lagrangian Modeling (PCLM) with Low-Rank Query Attention (LR-QA) (Zhang et al., 2017, Aminfar et al., 2014, Zeng et al., 13 Oct 2025).

1. Scope and unifying idea

Across the supplied sources, the label is attached to three distinct constructions. Their common motif is the replacement of globally expensive dense operations by compressed, low-rank surrogates targeted at repeated solves, multifrontal subproblems, or large point-cloud inference. This suggests a shared organizing principle: expensive global interactions are amortized by a one-time factorization, hierarchical compression, or covariance-based attention.

Formulation	Query object	Low-rank mechanism
Locally perturbed boundary solver	New, locally perturbed geometry	Low rank update to the original system plus Sherman–Morrison–Woodbury
3D elliptic multifrontal solver	Dense frontal matrices in multifrontal elimination	HODLR representation plus BDLR pseudo-skeleton compression
Transformer neural operator	Global design vector and large point cloud	PCLM plus LR-QA with covariance decomposition

A common misconception is that “low-rank” has a single technical meaning here. In the supplied formulations it refers, respectively, to numerically low-rank cut-and-paste updates in boundary integral equations, off-diagonal low-rank structure in hierarchical dense fronts, and low-rank structure in attention induced by second-order statistics of physical fields (Zhang et al., 2017, Aminfar et al., 2014, Zeng et al., 13 Oct 2025).

2. Boundary-integral LRQ-Solver for locally perturbed geometries

In the boundary-integral formulation, the continuous problem is the Laplace Dirichlet boundary value problem on a simply-connected domain $\Omega \subset \mathbb{R}^2$ with smooth boundary $\Gamma$:

$$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$

The solution is represented by a double-layer potential

$$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$

and the boundary limit yields the second-kind integral equation

$$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$

abbreviated as $K[\sigma]=f$ on $\Gamma$ (Zhang et al., 2017).

The geometric query is a local perturbation of an original boundary $\Gamma_0$. One removes a small patch $\Gamma_c \subset \Gamma_0$ with $N_c$ points, glues in a new patch $\Gamma$0 with $\Gamma$1 points, and denotes the remaining part by $\Gamma$2. The perturbed boundary is

$\Gamma$3

After Nyström or Galerkin discretization on $\Gamma$4, one obtains

$\Gamma$5

A fast direct solver is then built for $\Gamma$6, using a data-sparse factorization such as HBS, HODLR, or an $\Gamma$7-matrix approximation $\Gamma$8 satisfying $\Gamma$9, together with an approximate inverse $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$0 (Zhang et al., 2017).

The new perturbed system is rewritten as an extended block diagonal matrix plus a small update $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$1, acting on the unknown densities $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$2. The update collects the couplings caused by removing $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$3 and adding $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$4. By proxy-surface arguments or HBS leaf-box extraction, the subblocks

$$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$5

as well as the diagonal-with-zero-diagonal block $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$6, are numerically low rank. Hence

$$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$7

and, after discarding zero rows and columns, the perturbed system is represented in the compact form

$$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$8

The significance of this formulation is operational rather than merely algebraic: the local perturbation is treated as a query against a precomputed inverse, so the original factorization is reused instead of rebuilt.

3. Precomputation, Woodbury acceleration, and asymptotics

For the low-rank updated system $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$9, the core identity is the Sherman–Morrison–Woodbury formula

$$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$0

The efficient application sequence is: compute $$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$1; form $$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$2; factor $$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$3 once by LU or Cholesky; compute $$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$4; solve $$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$5; and return $$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$6 (Zhang et al., 2017).

The precomputation stage is performed once for $$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$7: discretize $$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$8 on $$u(x)=\int_\Gamma D(x,y)\,\sigma(y)\,ds(y),\qquad D(x,y)=\partial_{\nu_y}G(x,y),\qquad G=-(1/2\pi)\log|x-y|,$$9, build the fast direct solver, and form $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$0. The online stage is executed for each placement of the patch $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$1: discretize the new patch and identify removal indices; assemble $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$2 so that $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$3; compute $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$4 at cost $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$5; form $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$6 at cost $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$7; factor or invert $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$8 in $$-\frac12 \sigma(x)+\int_\Gamma D(x,y)\sigma(y)\,ds(y)=g(x),\qquad x\in\Gamma,$$9; and, for a given right-hand side, compute the solution by the Woodbury sequence in $K[\sigma]=f$0. The total per-query cost is therefore $K[\sigma]=f$1 (Zhang et al., 2017).

For the original factorization, the stated construction and memory costs are $K[\sigma]=f$2 or $K[\sigma]=f$3 for HBS/HODLR, and applying $K[\sigma]=f$4 to a vector also costs $K[\sigma]=f$5 or $K[\sigma]=f$6. The update becomes asymptotically faster when $K[\sigma]=f$7. For fixed-size local patches, where $K[\sigma]=f$8 is roughly constant, the reported precompute-stage speedup is approximately $K[\sigma]=f$9 relative to building a new fast direct solver from scratch, while the solve stage is comparable or slightly slower (Zhang et al., 2017).

The limitations are explicit. If the patch grows proportionally with $\Gamma$0, then $\Gamma$1 and the speedup disappears. In practice the method is most advantageous as long as the patch contains fewer than $\Gamma$2 points, often less than $\Gamma$3 of $\Gamma$4. Supported perturbations include translated patches of fixed shape, which are described as ideal, and locally refined meshes, which remain efficient as long as $\Gamma$5 stays a small fraction of $\Gamma$6. The same methodology is stated to extend to other boundary-integral operators, including Helmholtz and elasticity, and to 3D provided a fast inverse $\Gamma$7 can be built (Zhang et al., 2017).

4. HODLR/BDLR realization for 3D elliptic multifrontal systems

A second LRQ-Solver formulation addresses 3D elliptic PDEs through the dense frontal matrices generated by nested-dissection and multifrontal sparse elimination. At each node of the elimination tree, a dense frontal matrix $\Gamma$8 is partitioned at the top level into

$\Gamma$9

with $\Gamma_0$0 and $\Gamma_0$1 treated recursively and the off-diagonal blocks approximated by rank-$\Gamma_0$2 factors,

$\Gamma_0$3

Recursive subdivision yields a $\Gamma_0$4-level HODLR tree (Aminfar et al., 2014).

The low-rank compression is constructed by the boundary-distance low-rank approximation (BDLR), a pseudo-skeleton scheme based on graph distance in the sparse matrix graph. For an off-diagonal block $\Gamma_0$5, one forms a separator graph on the DOFs in $\Gamma_0$6, identifies boundary sets $\Gamma_0$7 and $\Gamma_0$8, computes graph distances $\Gamma_0$9 and $\Gamma_c \subset \Gamma_0$0, and selects rows and columns closest to the separator. With the selected index sets $\Gamma_c \subset \Gamma_0$1 and $\Gamma_c \subset \Gamma_0$2, one forms

$\Gamma_c \subset \Gamma_0$3

then computes an LU factorization with full pivoting of $\Gamma_c \subset \Gamma_0$4, truncates to numerical rank $\Gamma_c \subset \Gamma_0$5, and returns factors $\Gamma_c \subset \Gamma_0$6 and $\Gamma_c \subset \Gamma_0$7 so that $\Gamma_c \subset \Gamma_0$8 (Aminfar et al., 2014).

The direct solve is itself Woodbury-style. For

$\Gamma_c \subset \Gamma_0$9

one introduces auxiliary variables $N_c$0 and $N_c$1, factors the diagonal blocks, and solves a small Schur complement

$N_c$2

The same procedure recurses to leaf size $N_c$3. The stated complexity is $N_c$4 if the off-diagonal ranks remain approximately constant, or $N_c$5 in practice; when rank growth follows the 3D pattern $N_c$6, the total can rise to $N_c$7, though with a small prefactor (Aminfar et al., 2014).

This formulation can also be used as a GMRES preconditioner. The stated strategy is to use a low-tolerance HODLR factorization, with $N_c$8, so that GMRES converges in $N_c$9–$\Gamma$00 iterations to machine precision $\Gamma$01. Practical settings include leaf size $\Gamma$02–$\Gamma$03, direct-solve tolerance $\Gamma$04, preconditioner tolerance $\Gamma$05–$\Gamma$06, BDLR depth $\Gamma$07–$\Gamma$08 for $\Gamma$09 and up to $\Gamma$10–$\Gamma$11 for $\Gamma$12 on unstructured meshes, and a rank cap that switches to dense storage when $\Gamma$13 (Aminfar et al., 2014).

5. Transformer-based neural-operator LRQ-Solver

The 2025 LRQ-Solver is an end-to-end differentiable neural operator built on a transformer-style backbone for large-scale 3D PDEs. Its inputs are a global design parameter vector $\Gamma$14 and a large point cloud $\Gamma$15 representing the geometry. Internally it combines two submodules: Parameter Conditioned Lagrangian Modeling (PCLM), which encodes $\Gamma$16 into a latent control field $\Gamma$17, and Low-Rank Query Attention (LR-QA), which replaces standard $\Gamma$18 self-attention by a covariance-based approximation with complexity $\Gamma$19, where $\Gamma$20 (Zeng et al., 13 Oct 2025).

The architecture proceeds in five stages. First, the Parameter-Conditioned Encoder (PCE) maps $\Gamma$21 to $\Gamma$22 by cross-attention. In the stated implementation, the PCE uses $\Gamma$23-head cross-attention, $\Gamma$24 learnable queries, and $\Gamma$25. Second, $\Gamma$26 is concatenated with each spatial coordinate $\Gamma$27 to form an extended input $\Gamma$28. Third, a stack of $\Gamma$29 transformer-style layers processes these features; each layer uses LR-QA with $\Gamma$30 heads and $\Gamma$31 per head together with a residual MLP of width $\Gamma$32. Fourth, a $\Gamma$33-layer MLP decodes the final features into pseudo-physics fields such as $\Gamma$34. Fifth, system-level quantities such as drag coefficient, total force, and heat flux are obtained by differentiable control-volume integrals (Zeng et al., 13 Oct 2025).

PCLM is defined by treating the local state as

$\Gamma$35

with $\Gamma$36 a global latent control vector summarizing $\Gamma$37. The conservation laws are then enforced over the extended input space $\Gamma$38 through a physics residual

$\Gamma$39

For aerodynamic tasks, a pair-wise ranking loss

$\Gamma$40

is added to preserve monotonicity in drag predictions (Zeng et al., 13 Oct 2025).

LR-QA is motivated by the observation that, in many PDE fields, the key and value matrices admit a low-rank structure. At layer $\Gamma$41, one forms

$\Gamma$42

each in $\Gamma$43. Instead of forming the $\Gamma$44 similarity matrix, LR-QA computes the $\Gamma$45 covariances

$\Gamma$46

and thereby reduces attention complexity from $\Gamma$47 to $\Gamma$48. The stated approximation guarantee is

$\Gamma$49

with the interpretation that decaying singular values imply discarding only low-energy modes (Zeng et al., 13 Oct 2025).

The training setup uses AdamW with learning rate $\Gamma$50 and weight decay $\Gamma$51, a learning-rate decay by $\Gamma$52 at epoch $\Gamma$53, batch size $\Gamma$54 per A100 GPU (64 GB), and DDP over $\Gamma$55 GPUs. The total loss is

$\Gamma$56

with $\Gamma$57, $\Gamma$58, $\Gamma$59, and $\Gamma$60 (Zeng et al., 13 Oct 2025).

6. Quantitative behavior, invariance, and stated limitations

The three formulations report different performance regimes. For the locally perturbed boundary solver, fixed-size local patches yield an approximately $\Gamma$61 speedup in the precompute stage, while the solve stage is comparable or slightly slower; the method is most advantageous when the perturbation remains a small fraction of the full geometry (Zhang et al., 2017). For the HODLR/BDLR multifrontal solver, the reported benchmarks include a 3D unstructured beam front of size approximately $\Gamma$62 K, where HODLR-preconditioned GMRES with $\Gamma$63 takes approximately $\Gamma$64 s versus LU at $\Gamma$65 s, reducing iterations from $\Gamma$66 to $\Gamma$67; a $\Gamma$68 K cylinder-head front with approximately $\Gamma$69 speedup at $\Gamma$70; and a FETI local cube problem with a $\Gamma$71 K front and approximately $\Gamma$72–$\Gamma$73 speedup (Aminfar et al., 2014).

For the neural-operator LRQ-Solver, the DrivAer++ benchmark uses $\Gamma$74 car models, $\Gamma$75 deformable parameters, and $\Gamma$76 k surface points per sample, with a $\Gamma$77 train/validation/test split and metrics MSE, MAE, MaxAE, and MRE. The reported result is MSE $\Gamma$78, MAE $\Gamma$79, MaxAE $\Gamma$80, MRE $\Gamma$81, training time $\Gamma$82 h, and inference time $\Gamma$83 s, corresponding to a $\Gamma$84 relative MSE reduction and a $\Gamma$85 inference speedup over Transolver++ at $\Gamma$86 s. On the 3D Beam benchmark, using $\Gamma$87 beam configurations with point clouds ranging from $\Gamma$88 to $\Gamma$89 k points and a $\Gamma$90 split, the reported result is MAE_sub $\Gamma$91 MPa, MAE_all $\Gamma$92 MPa, training time $\Gamma$93 h, and inference time $\Gamma$94 ms, corresponding to a $\Gamma$95 MAE reduction over Geom-DeepONet and $\Gamma$96 faster training. The same source states that RegDGCNN runs OOM on the full geometry, whereas LRQ-Solver scales to $\Gamma$97 million points on one A100 (Zeng et al., 13 Oct 2025).

The ablations identify complementary contributions from LR-QA and PCLM. On DrivAer++, the baseline MLP has MSE $\Gamma$98 and MRE $\Gamma$99; adding LR-QA reduces this to $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$00 and $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$01; adding PCLM gives $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$02 and $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$03; and the full model gives $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$04 and $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$05. On 3D Beam, the baseline MAE_all is $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$06 MPa, with $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$07 for baseline plus LR-QA, $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$08 for baseline plus PCLM, and $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$09 for the full model (Zeng et al., 13 Oct 2025).

Discretization invariance is reported by varying the number of points from $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$10 to $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$11 k, while maintaining stable error, exemplified by MAE $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$12 MPa on 3D Beam, and approximately constant inference latency of about $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$13 ms (Zeng et al., 13 Oct 2025). The stated limitations are likewise formulation-specific. In the classical update solver, large patches drive $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$14 toward $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$15 and remove the speedup (Zhang et al., 2017). In the HODLR setting, 3D rank growth can increase total cost to $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$16 (Aminfar et al., 2014). In the neural-operator setting, the covariance approximation relies on a rapidly decaying singular spectrum; highly turbulent or multi-scale flows with full-rank interactions may require larger $$-\Delta u(x)=0,\quad x\in\Omega,\qquad u(x)=g(x),\quad x\in\Gamma.$$17 or hierarchical rank-adaptation, the current model addresses steady or quasi-steady PDEs, and topology changes or additional multiphysics may require domain-decomposition or hybrid graph/transformer architectures (Zeng et al., 13 Oct 2025).

Taken together, these formulations position LRQ-Solver as a name for PDE solvers that use low-rank structure to answer queries efficiently, whether the query is a local geometric perturbation, a multifrontal dense front, or a parameterized point-cloud configuration.