PaperlayerView exact revision on arXiv

Source-linked AI summary

Multigrid with rough coefficients and Multiresolution operator decomposition from Hierarchical Information Games

Houman Owhadi

arXiv:1503.03467v5math.NAcs.AImath.STstat.ML

TL;DR

The paper tackles the open problem of robust multigrid methods for PDEs with rough L∞ coefficients. It uses hierarchical information games to construct localized gamblet bases and an energy-orthogonal multiresolution decomposition, obtaining near-linear complexity and rigorous accuracy estimates.

  • Problem

    Provably robust multigrid methods for elliptic PDEs with rough L∞ coefficients remained an open practical problem because coefficient irregularity can severely impair convergence.

  • Method

    The method reformulates interpolation, recovery, and PDE-solution approximation as hierarchical information games, yielding localized gamblets and nested variational multiresolution systems.

  • Results

    The proposed method achieves near-linear complexity with H1 accuracy, including ϵ ∼N−1/d in O(N ln3d N) operations and subsequent-solve accuracy ϵ ≈N−1/d in O(N lnd+1 N) operations.

  • Takeaways & Limitations

    Gamblets provide energy-orthogonal sparse compression and decompose the PDE into independent linear systems with uniformly bounded condition numbers.

  • Takeaways & Limitations

    Sharper constants use convex measurement subdomains, while periodic or ergodic coefficient structure can limit gamblet computation to periodicity cells.

Abstract

from arXiv · show

We introduce a near-linear complexity (geometric and meshless/algebraic) multigrid/multiresolution method for PDEs with rough ($L^\infty$) coefficients with rigorous a-priori accuracy and performance estimates. The method is discovered through a decision/game theory formulation of the problems of (1) identifying restriction and interpolation operators (2) recovering a signal from incomplete measurements based on norm constraints on its image under a linear operator (3) gambling on the value of the solution of the PDE based on a hierarchy of nested measurements of its solution or source term. The resulting elementary gambles form a hierarchy of (deterministic) basis functions of $H^1_0(Ω)$ (gamblets) that (1) are orthogonal across subscales/subbands with respect to the scalar product induced by the energy norm of the PDE (2) enable sparse compression of the solution space in $H^1_0(Ω)$ (3) induce an orthogonal multiresolution operator decomposition. The operating diagram of the multigrid method is that of an inverted pyramid in which gamblets are computed locally (by virtue of their exponential decay), hierarchically (from fine to coarse scales) and the PDE is decomposed into a hierarchy of independent linear systems with uniformly bounded condition numbers. The resulting algorithm is parallelizable both in space (via localization) and in bandwith/subscale (subscales can be computed independently from each other). Although the method is deterministic it has a natural Bayesian interpretation under the measure of probability emerging (as a mixed strategy) from the information game formulation and multiresolution approximations form a martingale with respect to the filtration induced by the hierarchy of nested measurements.

1 Introduction

The paper addresses the open problem of designing multigrid methods provably robust to rough L∞ coefficients by reformulating solver design through information and decision games. It constructs a multiresolution method with H1 accuracy, near-linear complexity, and energy-orthogonal decompositions.

  • 1.1 Scientific discovery as a decision theory problem: The paper reformulates fast PDE-solver design as hierarchical information games in which incomplete measurements determine interpolation and recovery strategies.The formulation treats computation as operating on finite-dimensional features and uses nested measurements of the solution or source term.
  • 1.1 Scientific discovery as a decision theory problem: Multigrid convergence can be severely affected by nonsmooth coefficients, while provably robust methods for rough L∞ coefficients remained an open practical problem.The difficulty includes identifying interpolation and restriction operators adapted to coefficient microstructure.
  • 1.1 Scientific discovery as a decision theory problem: For ϵ ∼N−1/d, the proposed multiresolution method achieves ϵ-accuracy in H1 norm in O(N ln3d N) operations.For subsequent solves, it achieves accuracy ϵ ≈N−1/d in H1-norm in O(N lnd+1 N) operations.
  • 1.2 Outline of the paper: The decomposition splits the PDE into linear systems with uniformly bounded condition numbers and supports near-linear computation through orthogonality, nesting, and localization.The core decomposition can also support fast simulation or diagonalization of associated wave and parabolic equations.
  • 1.2 Outline of the paper: Gamblets form nested basis functions that orthogonally decompose H1_0(Ω) in the PDE energy scalar product and sparsely compress its solution space.The resulting multiresolution approximations form a martingale under the hierarchy of nested measurements.
  • 1.3 On the correspondence between statistical inference and numerical approximation: The paper connects statistical inference and numerical approximation through probabilistic conditioning, including spline interpolation and quadrature examples.The cubic spline interpolant appears within this correspondence between probabilistic models and numerical approximation.

2 Linear Algebra with incomplete information

The section formulates recovery from incomplete linear measurements as an information-game problem and derives Bayesian, variationally optimal approximations in the energy norm. With Q = A, the conditional estimator is the Galerkin best approximation, while measurement choice controls the worst-case error.

  • Recovery problem: Incomplete measurements cannot generally recover x exactly when the measurement rank m is smaller than n.The recovery problem assumes known measurements Φx and norm constraints on the unknown source or solution.
  • Bayesian recovery: The Bayesian approximation is the conditional expectation E[X|ΦX = y], obtained after replacing the deterministic source by a centered Gaussian vector.The conditional law has mean Ψy and covariance KΦ.
  • Variational recovery: For w ∈ R^m, Ψw uniquely minimizes the quadratic energy subject to the measurement constraints Φv = w.The estimator reproduces the measurements and its image is orthogonal to Ker(Φ) in the relevant scalar product.
  • Variational recovery: Theorem 2.3 makes ΨΦv the unique best approximation of v in the K^-1 norm from the measurement-induced approximation space.For the original solution, the same estimator minimizes error over all admissible measurement coefficients.
  • Energy-norm selection: With A symmetric positive definite and Q = A, Ψy is the Galerkin best approximation of x in the A-energy norm without prior knowledge of b.The approximation space is the image of A^-1Φ^T.
  • Measurement selection: Choosing Φ to measure the projection onto the first m eigenvectors of A gives the best measurement strategy in the stated eigenbasis analysis.The supplied passage introduces this comparison but its displayed error value is truncated.

3 Numerical homogenization and design of the interpolation operator in the continuous case

The continuous formulation uses a Gaussian-field information game to construct interpolation functions satisfying measurement constraints while minimizing the PDE energy. The resulting gamblets have optimal recovery properties, exponential decay, and localized approximation with preserved O(H) energy convergence.

  • Game formulation: The continuous interpolation operator is identified through a non-cooperative game in which source terms are opposed to solution recovery from finitely many measurements.The method extends the incomplete-information analysis to numerical homogenization.
  • Gaussian formulation: A centered Gaussian source is selected so that the resulting field has the PDE Green’s function as covariance and the differential operator as precision.This choice preserves the energy-based structure used to define the interpolation basis.
  • Gamblets: Gamblets are basis functions dual to the measurement functions, satisfying ∫Ω ψ_i φ_j = δ_i,j and representing bets on the PDE solution.They are obtained by conditioning the Gaussian solution on the measurements.
  • Variational characterization: Each gamblet is the unique minimizer of the PDE energy subject to its measurement constraints, and linear combinations inherit optimal variational recovery properties.The associated approximation u* is the best energy-norm approximation in the span of the gamblets or equivalent Green-operator images.
  • Optimal recovery: The gamblet approximation u* is the best energy-norm approximation of u in span{L^-1φ_i : i ∈ {1, . . . , m}}.The gamblet span equals the span of the corresponding inverse-operator images.
  • Localization: Gamblets decay exponentially away from their measurement supports, enabling localized computation on subdomains.Localization preserves the O(H) energy-norm convergence rate when subdomains have size O(H ln(1/H)).

4 Multiresolution operator decomposition

Nested games based on measurements at multiple resolutions yield an orthogonal multiresolution decomposition of the PDE and a near-linear complexity algorithm with a-priori error bounds.

  • Multiresolution operator decomposition: Hierarchical nested games use measurements at different resolutions to derive an orthogonal-across-subscales multiresolution decomposition.The same construction provides a multiresolution algorithm with a-priori error bounds.

4.1 Hierarchy of nested measurement functions

The hierarchy is built from nested measurement functions supported on nested partitions whose cells have controlled geometry and scale. Convexity is used only to sharpen or simplify constants, not as an essential requirement.

  • Hierarchy definition: An index tree organizes measurement functions across q resolution levels and records parent-child relationships between indices.The notation i(k) and I(k) identifies prefixes and level-wise index sets.
  • Nested measurements: The partitions are nested across levels, so every finer cell belongs to a parent cell at the preceding scale.This nesting induces a corresponding hierarchy of measurement functions.
  • Geometric assumptions: Convexity of the subdomains is not necessary for the stated results and primarily provides sharper or simpler constants.Without convexity, approximation bounds remain valid with a multiplicative factor involving π.
  • Normalization: The construction may assume level-wise constant measurement-function magnitudes without loss of generality by rescaling functions and transfer coefficients.The general case requires tracking dependence on the corresponding maximum magnitudes.

4.2 Hierarchy of nested gamblets and multiresolution approximations

Nested measurements generate nested approximation spaces whose gamblet bases provide best energy-norm approximations of the PDE solution. These spaces are characterized through L^-1 applied to measurement functions.

  • u(k) is the best energy-norm approximation of the PDE solution within V(k).
  • Nested measurement functions induce nested spaces V(k) ⊂ V(k+1).

4.3 Nested games and martingale/multiresolution decomposition

A hierarchy of nested measurement games produces multiresolution approximations that form a martingale and decompose the solution into orthogonal scale components. The resulting spaces and increments yield an orthogonal decomposition of H^1_0(Ω).

  • Nested measurements from coarse to fine produce approximations u(k) that form a martingale with independent increments.
  • Each u(k) is the finite element solution in V(k), and u(q+1) equals the full solution u.
  • The full solution is represented as the sum of its successive multiresolution increments.
  • The spaces decompose as V(1) ⊕a W(2) ⊕a ··· ⊕a W(q+1), with each increment u(k+1)−u(k) lying in W(k+1).

4.4 Interpolation and restriction matrices/operators

Nested approximation spaces induce restriction matrices and transpose interpolation operators between adjacent resolutions. Their entries represent bets on finer-scale solution values conditioned on coarser information.

  • Nested spaces define restriction matrices R(k,k+1) between adjacent measurement levels.
  • The transpose R(k+1,k) serves as the interpolation or prolongation matrix.
  • Matrix entries are interpreted as Player II’s best bets on finer-scale solution values.
  • The restriction identities follow from integrating the nested basis relations against finer-level measurement functions and conditioning on the coarser filtration.

4.5 Nested computation of the interpolation and stiffness matrices

The interpolation and stiffness matrices can be computed hierarchically across resolution levels using covariance, restriction, and prolongation relations. A variational characterization identifies prolongation coefficients as unique minimizers.

  • The covariance vector of measurements is Gaussian with a symmetric positive definite covariance matrix Θ(k).
  • Gamblets admit representation formulas, and Θ(k)^-1 equals the stiffness matrix A(k).
  • The hierarchy computes A(k) from A(k+1) through nested interpolation and restriction relations.
  • For any coarse vector b, prolongation produces the unique fine-scale coefficient vector minimizing the stated variational objective.
  • Restriction and prolongation satisfy inverse-type identities involving π(k,k+1), while Θ(k) is obtained by coarse-scale projection of Θ(k+1).

4.6 Multiresolution gamblets

The section constructs multiresolution gamblet bases by choosing local, interscale-detail matrices W^(k) whose rows span kernels of restriction maps and form bases of W^(k).

  • Multiresolution gamblets: Gamblets χ_i^(k) form a basis of the detail space W^(k), whose dimension is |I^(k)| − |I^(k−1)|.The basis vectors lie in the kernel of the restriction map and are linearly independent.
  • Multiresolution gamblets: The matrices W^(k) are sparse across parent blocks, enabling fast multiplication by W^(k) and its transpose.Nonzero entries are restricted to pairs with matching parent indices.
  • Multiresolution gamblets: Each W^(k) is built from block-local vectors that are linearly independent and orthogonal to the all-ones vector within each parent block.The construction separates detail degrees of freedom among descendants sharing the same parent index.
  • Multiresolution gamblets: The hierarchy admits a game-theoretic interpretation, with the gamblet construction linked to the nested partition shown in Figure 1.The passage identifies the interpretation but does not provide the figure’s full visual encoding.
  • Multiresolution gamblets: Construction 4.13 has complexity |I^(k−1)| × m_s^2 and satisfies W^(k)W^(k),T = J^(k), the identity on detail coordinates.Thus its rows are orthonormal in the Euclidean coefficient representation.

4.7 Multiresolution operator inversion

The operator inversion decomposes the finite-element solution into a coarse component and independent subband corrections obtained from recursively computed right-hand sides and decoupled systems.

  • Multiresolution operator inversion: For each k ≥ 2, the subband coefficient vector w^(k) solves a system involving the detail stiffness matrix B^(k).The section defines B^(k) as the stiffness matrix on the detail index set J^(k).
  • Multiresolution operator inversion: The right-hand-side vectors g^(k) are computed iteratively from the source term across the hierarchy.The construction therefore propagates forcing information through nested scales rather than assembling each scale independently.
  • Multiresolution operator inversion: The coarse coefficient vector U^(1) is obtained by solving the scale-one system A^(1)U^(1) = g^(1).This provides the coarse component used alongside the subband solves.
  • Multiresolution operator inversion: The solution at any scale is computed by solving the decoupled linear systems (4.26) and (4.27).The theorem identifies these systems as sufficient for computing the finite-element solution and its scale increments.
  • Multiresolution operator inversion: The scale increment satisfies u^(k) − u^(k−1) = P^(k)w^(k), yielding a hierarchical reconstruction of u^(k).The increment is represented in the detail basis and added to the preceding-scale solution.

4.8 Uniformly bounded condition numbers across subscales/subbands

The paper establishes an orthogonal subscale decomposition of −div(a∇) and bounds the conditioning of the resulting restricted operators using geometric properties of the hierarchy and coefficient ellipticity.

  • 4.8 Uniformly bounded condition numbers across subscales/subbands: The operator −div(a∇) decomposes into layered subbands, with uniformly bounded subband condition numbers when H^(k−1)/H^k is uniformly bounded.A geometric sequence of scales is given as a sufficient example.
  • 4.8 Uniformly bounded condition numbers across subscales/subbands: Theorem 4.16 bounds the energy-to-L2 behavior on V^(k) and on each detail space W^(k).The proof uses the decomposition V^(k) = V^(k−1) ⊕_a W^(k) together with approximation estimates.
  • 4.8 Uniformly bounded condition numbers across subscales/subbands: Theorem 4.17 derives condition-number bounds for the scale matrices and detail stiffness matrices from mesh scales, ellipticity, and partition parameters.The proof bounds the smallest and largest eigenvalues through the energy estimates for V^(k) and W^(k).
  • 4.8 Uniformly bounded condition numbers across subscales/subbands: Under Construction 4.13, Cond(W^(k)W^(k),T) = 1; under Construction 4.11, it is bounded by 2^d.The latter construction also receives a bound based on the local block sizes m_j.
  • 4.8 Uniformly bounded condition numbers across subscales/subbands: The matrices P^(k) are projections, and their action is analyzed through the decomposition of coefficient vectors into detail and coarse components.The analysis uses the orthogonality of Im(W^(k),T) and Im(π^(k,k−1)).

4.9 Well conditioned relaxation across subscales

The well-conditioned systems used for gamblets and subband solutions can be solved efficiently with iterative methods because their conditioning is independent of mesh resolution and coefficient regularity.

  • 4.9 Well conditioned relaxation across subscales: For geometric H^k, the systems for gamblets and subband solutions have uniformly bounded condition numbers.The bound is independent of mesh size or resolution and of the regularity of a(x).
  • 4.9 Well conditioned relaxation across subscales: These uniformly conditioned systems can be solved efficiently using iterative methods such as Conjugate Gradient.The section connects the conditioning result directly to iterative solution of the linear systems.
  • 4.9 Well conditioned relaxation across subscales: Conjugate Gradient reduces the error according to its standard condition-number-dependent iteration bound.The cited passage introduces the error criterion and the associated complexity dependence on Cond(A) and the number of nonzeros.

4.10 Hierarchical localization and error propagation across scales

The section localizes gamblet computations to recover near-linear complexity and controls how localization errors propagate across scales. With sufficiently large localization radii, the localized method preserves approximation accuracy and uniformly bounded conditioning.

  • Motivation: Dense multiresolution matrices must be truncated by localizing basis-function computations to achieve near-linear complexity.The localization relies on the exponential decay of gamblets and requires controlling error propagation across scales.
  • Hierarchical localization: Localization neighborhoods are defined through index sets of subdomains within radius ρ at each scale, and localized gamblets are constructed hierarchically.The localized spaces and basis functions are then used to define localized coarse-to-fine solution components.
  • Error propagation: ∥R∥2 ≤ CH^(d/2)e^(−ρk/C), showing that interscale truncation effects decay exponentially with localization radius.This bound is used to propagate localization errors from finer to coarser scales.
  • Localized systems: The localized matrices retain controlled conditioning, with Cond(A(1),loc) ≤ CH^−2 and Cond(B(k),loc) ≤ CH^−2−2d under the stated error conditions.Theorem 4.24 also bounds the resulting energy-norm solution errors in terms of the localization error E(k,χ).
  • Accuracy guarantee: If ρk satisfies the threshold in Theorem 4.27, then ∥u−u(k),loc∥a ≤ C(Hk+ϵ)∥g∥L2(Ω) and each localized interscale error is at most ϵ/2k^2∥g∥H−1(Ω).The same theorem guarantees the localized systems' condition-number bounds at all scales.

5 The algorithm, its implementation and complexity

Algorithm 2 constructs localized gamblets across nested resolutions and solves the resulting multiresolution systems with prescribed H1 accuracy. Its near-linear performance rests on nesting, bounded condition numbers, and exponential-decay-based truncation.

  • Initialization: The discretization uses a regular fine mesh with localized finite-element nodal basis elements and accommodates piecewise-constant rough coefficients.The method is also meshless/algebraic because it only requires specification of the basis elements.
  • Initialization: Resolutions form a geometric hierarchy, with nested index sets and aggregated measurement functions defining the multilevel construction.The example uses H = 1/2 and q = 6, with the finest level corresponding to the fine mesh.
  • Exact gamblet transform and multiresolution operator inversion: The exact algorithm decomposes the poorly conditioned fine-grid system into nested linear systems with uniformly bounded condition numbers and computes independent subband corrections.Subband solutions represent behavior over successive scale bands, and their summation recovers the multiresolution approximations.
  • Localized algorithm: Algorithm 2 localizes and truncates dense systems using exponential decay, preserving prescribed approximation accuracy while enabling near-linear computation.The three enabling properties are nesting, uniformly bounded condition numbers, and localization/truncation based on exponential decay.
  • Accuracy guarantees: The localized method provides H1 error bounds for the solution and subband corrections while retaining bounded condition-number estimates for its localized matrices.The global approximation satisfies ∥u − uloc∥H1_0(Ω) ≤ ϵ∥g∥H−1(Ω).
  • Complexity: The first solve and subsequent solves have near-linear complexity, while localized basis coefficients can be sparsified without loss in energy norm.Subsequent solves reuse gamblets and stiffness matrices rather than recomputing them.
Loading 1503.03467v5…