Generalized Rough Polyharmonic Splines (GRPS)

• GRPS are operator-adapted basis functions for numerical homogenization that extend classic polyharmonic splines to handle rough coefficients and divergence-form elliptic operators.
• Their formulation incorporates Bayesian inference and constrained energy minimization to derive coarse spaces that avoid assumptions like periodicity while achieving O(H) energy-norm accuracy.
• Localization techniques yield exponentially decaying basis functions, enabling efficient repeated solves in multiscale PDEs and nonlinear problems such as the Landau–Lifshitz equation.

Generalized Rough Polyharmonic Splines (GRPS) are operator-adapted basis functions for numerical homogenization of PDEs with rough coefficients. They extend polyharmonic splines from Laplacian-based settings to divergence-form elliptic operators $L=\operatorname{div}(a\nabla\cdot)$ with merely $L^\infty$ coefficients, and in later formulations to general integro-differential operators and measurement-constrained coarse spaces derived from Bayesian conditional expectation or equivalent constrained energy minimization. The framework is designed to avoid assumptions of periodicity, scale separation, or ergodicity, while retaining $O(H)$ energy-norm accuracy, exponential localization on patches of size $O(H\log(1/H))$, and reuse of a precomputed coarse basis across repeated solves. Recent work extends GRPS to the fully nonlinear Landau–Lifshitz equation with rough, non-periodic, and nonseparable coefficients (Owhadi et al., 2012, Liu et al., 2021, Chen et al., 2019, Ma et al., 4 Aug 2025).

1. Foundational operator-theoretic definition

In its original rough-coefficient form, GRPS are defined for the elliptic problem

$$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$

with $a(x)$ symmetric, uniformly elliptic, and in $L^\infty(\Omega)$. For $d\leq 3$, the solution space is

$$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$

where $L:=\operatorname{div}(a\nabla\cdot)$. Given scattered interpolation points $L^\infty$0 with mesh norm

$L^\infty$1

the $L^\infty$2-th basis function is the unique minimizer

$L^\infty$3

Away from the interpolation nodes, the Euler–Lagrange equation is

$L^\infty$4

with boundary conditions $L^\infty$5 and $L^\infty$6 on $L^\infty$7. For $L^\infty$8, the resulting functions are therefore biharmonic with respect to the rough operator $L^\infty$9; for $O(H)$0, one chooses $O(H)$1 with $O(H)$2 and minimizes $O(H)$3, leading to $O(H)$4-harmonicity with respect to $O(H)$5 (Owhadi et al., 2012).

The interpolation and Galerkin properties are central. If $O(H)$6 solves $O(H)$7 and $O(H)$8, then

$O(H)$9

The corresponding Galerkin solution $O(H\log(1/H))$0 in $O(H\log(1/H))$1 satisfies the same $O(H\log(1/H))$2 energy-norm estimate. The analysis uses a higher-order Poincaré inequality for functions vanishing at the interpolation nodes,

$O(H\log(1/H))$3

A notable structural feature is that the error depends on the fill distance $O(H\log(1/H))$4, not on coarse-element aspect ratios or regularity of node placement (Owhadi et al., 2012).

2. Bayesian reformulation and constrained energy minimization

A generalized formulation casts numerical homogenization as a Bayesian inference problem for

$O(H\log(1/H))$5

with $O(H\log(1/H))$6 mapping a Hilbert space $O(H\log(1/H))$7 into $O(H\log(1/H))$8. The prototype treated in detail is the rough elliptic operator

$O(H\log(1/H))$9

where $$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$0 is symmetric positive definite and uniformly elliptic (Liu et al., 2021).

The Bayesian construction randomizes the right-hand side: $$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$1 with $$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$2 a centered Gaussian field. Two prior choices are emphasized: white noise, $$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$3, and “$$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$4-noise,” for which the covariance operator of $$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$5 equals $$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$6. For a finite family of linearly independent measurement functions $$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$7, the measurements are

$$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$8

and the posterior mean is

$$-\operatorname{div}(a(x)\nabla u)=g \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\partial\Omega,$$9

The span of $a(x)$0 is the GRPS coarse space (Liu et al., 2021).

The same basis arises variationally. With

$a(x)$1

each basis function is the unique minimizer

$a(x)$2

This formulation yields a strict convexity statement, a saddle-point system, and the orthogonality relation

$a(x)$3

Within this perspective, the prior determines the energy being minimized, the measurements determine the coarse information being preserved, and the posterior mean is the minimum-variance estimator as well as the minimal-energy field consistent with those measurements (Liu et al., 2021).

3. Measurement families and the structure of coarse information

The generalized theory replaces nodal constraints by coarse observables on a triangulation $a(x)$4. Three families of measurement functions are treated explicitly, together with a combined form that is preferred in implementation because of linear independence (Liu et al., 2021).

Family	Measurement function(s)	Comment
V	$a(x)$5, $a(x)$6	Volume averages on coarse elements
E	$a(x)$7, $a(x)$8	Edge averages on coarse edges
D'	$a(x)$9	First-order derivative measurements
D	$L^\infty(\Omega)$0	Same span as $L^\infty(\Omega)$1 for Dirichlet problems; linearly independent

For derivative measurements,

$L^\infty(\Omega)$2

and the edge–derivative equivalence is made explicit through

$L^\infty(\Omega)$3

when $L^\infty(\Omega)$4. This provides a concrete realization of derivative information through edge combinations (Liu et al., 2021).

The choice of measurements has direct approximation consequences. In the multiscale PDE analysis with rough coefficients, GRPS-E uses edge averages, GRPS-V uses volume averages, and GRPS-D incorporates higher-order information. The same source reports that edge averages are natural when elongated channels or cracks imply that edge integrals better capture the physics, whereas derivative measurements encode higher-order information through combinations of edge data. In numerical experiments on the multiscale trigonometric coefficient, GRPS-E and GRPS-D produce bases that are more localized and “spiky” than GRPS-V and classical RPS, indicating faster spatial decay (Liu et al., 2021).

4. Localization, exponential decay, and approximation theory

Global GRPS basis functions are supported on all of $L^\infty(\Omega)$5, so practical computation relies on localization. If $L^\infty(\Omega)$6 is the smallest union of coarse elements containing $L^\infty(\Omega)$7, the $L^\infty(\Omega)$8-layer patch is defined recursively by

$L^\infty(\Omega)$9

The localized basis is then

$d\leq 3$0

The localized space is $d\leq 3$1 (Liu et al., 2021).

The localization theory is based on additive Schwarz decomposition. Under uniform ellipticity and shape regularity, one obtains exponential decay of the corrector and hence of the global basis away from its measurement support. A representative localization theorem states that there exist constants $d\leq 3$2, depending only on $d\leq 3$3, such that

$d\leq 3$4

The corresponding localized coarse solve

$d\leq 3$5

retains the global $d\leq 3$6 rate when $d\leq 3$7. The analysis combines Galerkin orthogonality with measurement-dependent inequalities: standard Poincaré estimates for volume measurements, zero mean boundary trace inequalities of Nazarov–Repin type for edge measurements, and mixed arguments for combined measurements (Liu et al., 2021).

The earlier rough polyharmonic spline theory gives a complementary localization statement for pointwise constraints. Localized basis functions $d\leq 3$8 are computed by minimizing $d\leq 3$9 over a local admissible set with point constraints inside $$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$0. The a posteriori estimate

$$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$1

contains explicit boundary-layer terms $$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$2, and numerical evidence shows exponential decay of these tails. Choosing $$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$3 with thickness $$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$4 yields $$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$5, so the optimal $$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$6 energy-norm accuracy is preserved. Because each localized basis is supported only on its patch, the resulting linear systems are sparse and banded (Owhadi et al., 2012).

The same localization mechanism underlies repeated-solve settings. In coarse optimal-control discretizations, localized basis functions on overlapping patches of diameter $$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$7 are precomputed once and reused in all state and adjoint solves. This decouples the fine-scale heterogeneity from the subsequent iterative solves (Chen et al., 2019).

5. Relation to adjacent multiscale methods and established applications

GRPS generalize several earlier operator-adapted constructions. In the generalized Bayesian formulation, RPS correspond to white noise and pointwise measurements, while Gamblets correspond to $$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$8-noise and volume averages. GRPS retain the same posterior-mean and constrained-minimization structure but enlarge the admissible measurement families to include edge averages and first-order derivative information. Compared with GMsFEM, GRPS do not rely on local spectral problems; compared with LOD, they can be placed in the same constrained coarse-space framework by choosing suitable measurement functions or quasi-interpolation; compared with MsFEM, their analysis is not tied to cell problems, scale separation, or ergodicity (Liu et al., 2021, Chen et al., 2019).

Several application domains are explicitly documented. For elliptic optimal control with rough $$V:=\{u\in H_0^1(\Omega):Lu\in L^2(\Omega)\},\qquad \|u\|_V:=\|Lu\|_{L^2(\Omega)},$$9 coefficients, the state and co-state both solve multiscale elliptic equations of the form

$L:=\operatorname{div}(a\nabla\cdot)$0

and GRPS provide a coarse space $L:=\operatorname{div}(a\nabla\cdot)$1 of dimension $L:=\operatorname{div}(a\nabla\cdot)$2 such that the localized coarse solution achieves

$L:=\operatorname{div}(a\nabla\cdot)$3

In the optimality system, the discrete state, adjoint, and control errors satisfy

$L:=\operatorname{div}(a\nabla\cdot)$4

The computational strategy is a one-time, parallelizable precomputation of localized basis functions followed by repeated low-dimensional coarse solves; the projected-gradient iteration on the coarse system converges linearly with rate $L:=\operatorname{div}(a\nabla\cdot)$5 under the stated assumptions (Chen et al., 2019).

Time-dependent and inverse problems are also part of the original scope. The same spatial GRPS basis is used for parabolic equations $L:=\operatorname{div}(a\nabla\cdot)$6 and hyperbolic equations $L:=\operatorname{div}(a\nabla\cdot)$7, with coarse approximations of the form $L:=\operatorname{div}(a\nabla\cdot)$8. For recovery from point measurements, if $L:=\operatorname{div}(a\nabla\cdot)$9 and $L^\infty$00, then

$L^\infty$01

This makes GRPS a reconstruction device as well as a homogenization basis (Owhadi et al., 2012).

The numerical record emphasizes robustness in highly heterogeneous media. On the multiscale trigonometric coefficient with contrast $L^\infty$02, fixed-localization studies show saturation at the optimal $L^\infty$03 rate, while for $L^\infty$04 the GRPS-D space exhibits approximately second-order convergence and requires fewer layers to reach saturation than RPS, GRPS-V, or GRPS-E. On the SPE10 benchmark, GRPS-D achieves the best accuracy for the same degrees of freedom and reaches stable convergence with fewer layers, including a reported layer-63 test with contrast $L^\infty$05. In a wave equation test with heterogeneous stationary medium, errors measured in $L^\infty$06 show nearly linear convergence versus coarse degrees of freedom when $L^\infty$07 is large enough (Liu et al., 2021).

6. Adaptation to the Landau–Lifshitz equation with rough coefficients

A recent extension applies GRPS to the fully nonlinear Landau–Lifshitz equation

$L^\infty$08

with

$L^\infty$09

The coefficient $L^\infty$10 is rough, non-periodic, and nonseparable, so direct fine-mesh discretization is prohibitively expensive, while classical HMM approaches that depend on cell problems and scale separation are not robust in this setting (Ma et al., 4 Aug 2025).

The GRPS basis is derived from the Landau–Lifshitz energy

$L^\infty$11

Given measurement functions $L^\infty$12, the global and localized basis problems are

$L^\infty$13

and

$L^\infty$14

Because anisotropy is aligned with $L^\infty$15, the localized energy penalizes $L^\infty$16 and $L^\infty$17 but not $L^\infty$18 (Ma et al., 4 Aug 2025).

The time-discrete variational problem is written as

$L^\infty$19

with coercive part

$L^\infty$20

while $L^\infty$21 collects the scheme-dependent terms for the Cimrák, Gao, or An time discretizations. The precession term is skew-symmetric and does not define an energy. Accordingly, basis construction uses the coercive part of the operator, and the full $L^\infty$22 is then solved in the reduced space. To reflect the anisotropic vector structure, the construction splits by component: $L^\infty$23 for the first component, and

$L^\infty$24

for the second and third components (Ma et al., 4 Aug 2025).

The coarse space is $L^\infty$25 for each component, assembled into $L^\infty$26. The reduced solve is

$L^\infty$27

with

$L^\infty$28

component-wise. Initial data are projected by

$L^\infty$29

An $L^\infty$30-projection $L^\infty$31 defined by

$L^\infty$32

is used to turn 4-valence tensors into 3-valence ones for efficient nonlinear assembly (Ma et al., 4 Aug 2025).

The reported approximation and performance properties combine the general GRPS theory with Landau–Lifshitz-specific experiments. The basis functions satisfy exponential decay,

$L^\infty$33

and the coarse approximation obeys

$L^\infty$34

In the numerical study, with $L^\infty$35, both GRPS-E and GRPS-V show roughly first-order convergence in $L^\infty$36; with $L^\infty$37, GRPS-V achieves roughly second-order $L^\infty$38 accuracy, and GRPS-E often shows improved rates close to second order in the reported tests. The basis does not need recomputation at each time step: its construction is offline and parallelizable. On an example with $L^\infty$39 and fine reference $L^\infty$40, the fine solve time is about $L^\infty$41s, GRPS-E about $L^\infty$42s, and GRPS-V about $L^\infty$43s, corresponding to reductions of approximately $L^\infty$44 and $L^\infty$45, respectively, for comparable accuracy. A larger timing table reports wall-clock reductions above $L^\infty$46 when GRPS-V is combined with accelerated assembly, replacing 4-valence tensor costs measured in days by coarse-space assembly measured in seconds (Ma et al., 4 Aug 2025).

The Landau–Lifshitz extension clarifies a frequent misunderstanding about GRPS in nonlinear systems. The method does not require the full operator to be symmetric or energy generating. In this case, the coarse basis is built from the coercive multiscale part—exchange plus anisotropy—while the skew-symmetric precession term is retained only in the reduced solve. This preserves the defining GRPS principle: basis functions are tailored to the multiscale structure that controls the dominant energy, and localization then converts that structure into a computationally tractable coarse model (Ma et al., 4 Aug 2025).