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

1.1 Scientific discovery as a decision theory problem

• Multigrid convergence can be severely affected by rough coefficients, and provably robust methods remain an open practical problem; identifying suitable interpolation operators is a central difficulty.Piecewise linear finite elements can perform arbitrarily badly, so interpolation bases must adapt to the coefficient microstructure.Source: §1.1 · p. 2, §1.2 · p. 3
• For epsilon ~ N^-1/d, the proposed method achieves epsilon-accuracy in H1 in O(N ln^(3d) N) operations; hierarchical matrices achieve epsilon-accuracy in L2 in O(N ln^(2d+8) N).For subsequent solves, the method achieves accuracy epsilon ≈ N^-1/d in H1 in O(N ln^(d+1) N) operations.Source: §1.1 · p. 3
• The method decomposes H1_0 into energy-orthogonal or near-orthogonal subspaces on which the PDE has uniformly bounded condition numbers.The introduction identifies this fast decomposition as the mechanism supporting the method’s complexity.Source: §1.1 · p. 3

1.2 Outline of the paper

• The interpolation problem is cast as a game over incomplete solution measurements; Gaussian-prior analysis yields deterministic gamblets with optimal recovery properties, exponential decay, and robustness to high contrast.The game has one player choose an admissible source term and another approximate the solution from measurements.Source: §1.2 · pp. 4–5
• Under a mixed strategy, nested measurements produce multiresolution approximations forming a martingale, while gamblets yield nested energy-orthogonal bases and sparse solution-space compression.The resulting decomposition consists of independent linear systems with uniformly bounded condition numbers.Source: §1.2 · pp. 5–6

1.3 On the correspondence between statistical inference and numerical approximation

• Bayesian numerical quadrature recovers the trapezoidal rule under a Brownian-motion prior and cubic spline interpolation under an integrated Brownian prior, illustrating the link between inference and approximation.Carefully chosen priors can also formulate classical min-max quadrature rules as Bayesian inference problems.Source: §1.3 · pp. 6–7

2.1 The recovery problem

• Identifying interpolation operators is equivalent to recovering a linear-system solution from m<n coarse measurements using known norm constraints on the source image.Because m<n, the measurements are not sufficient a priori to recover the exact solution.Source: §2.1 · pp. 7–8

2.2 Player I’s mixed strategy

• With a Gaussian prior for the unknown source, Bayesian recovery approximates x by E[X|ΦX=y]=Ψy, and the conditional solution is Gaussian with covariance KΦ.The conditional covariance KΦ has rank n−m.Source: §2.2 · p. 8

2.3 Variational/optimal recovery properties and approximation error

• Ψy uniquely minimizes ∥Av∥Q^-1 subject to Φv=y, and ΨΦv minimizes the K^-1-norm error over approximations Ψw.The conditional variance of solution coordinate i controls the accuracy of its deterministic approximation.Source: §2.3 · pp. 9–10

2.4 Energy norm estimates and selection of the prior

• When A is symmetric positive definite and Q=A, Ψy is the Galerkin approximation and best approximation in the energy norm within Im(Ψ)=Im(A^-1ΦT), without knowing b.The energy error equals the minimum over vectors in Im(A^-1ΦT).Source: §2.4 · p. 10

2.5 Impact and selection of the measurement matrix Φ

• For symmetric positive definite A with Q=A, Theorem 2.9 allows measurements identified for B with α1B≤A to be reused for A, with an increased error bound.The bound increases by a multiplicative factor α1^-1/2.Source: §2.5 · p. 11
• When λmin(a)=1, the finite-element error bound on ∥uh−uH∥a does not depend on the coefficient contrast λmax(a)/λmin(a).Source: §2.5 · p. 12

3.1 Information Game and Gamblets

• Conditioning a Gaussian solution field on measurements yields gamblets that form a basis for betting on the PDE solution from those measurements.The noise covariance is chosen as the differential operator L, which acts as the field’s precision operator.Source: §3.1 · pp. 13–14

3.2 Optimal recovery properties

• Each gamblet is the unique minimizer of energy subject to measurement constraints, and the resulting functions have optimal variational recovery properties.The optimization enforces the measurement constraints on the basis function.Source: §3.2 · p. 14

3.3 Optimal accuracy of the recovery

• The recovery u* is the best energy-norm approximation of the solution in the span of the measurement-source responses L^-1φ_i.The stated approximation space is also the span of the gamblets.Source: §3.3 · p. 15

3.4 Transfer property and selection of the measurement functions

• A transfer property supports accurate finite-dimensional measurement bases, including constructions using functions supported on individual subdomains.Convexity of the subdomains sharpens constants; without convexity, approximation bounds remain valid after multiplication by π.Source: §3.4 · pp. 15–17

3.5 Exponential decay of gamblets

• Gamblets decay exponentially with distance from their measurement supports, so their computation can be restricted to subdomains of diameter O(H ln 1/H) without changing the order of accuracy.This localization claim is given for the construction in which each measurement subdomain contains a ball of diameter δH.Source: §3.5 · p. 17

3.6 Localization of the basis elements

• Localized basis elements are defined by minimizing a quadratic problem on a patch and extending the resulting function by zero outside it.The patch is formed from subdomains intersecting a ball around a point in the measurement subdomain.Source: §3.6 · p. 19
• Localization to subdomains of size O(H ln(1/H)) preserves gamblets’ O(H) convergence rate in the energy norm.The result supports using these localized gamblets as basis functions in numerical homogenization.Source: §3.6 · p. 21
• The authors identify localized gamblets as basis functions for numerical homogenization [85].Source: §3.6 · p. 21
• The authors state that Section 4 will show how to compute gamblets hierarchically at near-linear complexity.Source: §3.6 · p. 21

4 Multiresolution operator decomposition

• Nested games gamble on approximating the PDE solution from measurements at different levels of resolution.Source: §4 · p. 21
• The hierarchical games derive a multiresolution decomposition of the PDE that is orthogonal across subscales.Source: §4 · p. 21
• The resulting multiresolution algorithm has near-linear complexity and a priori error bounds.Source: §4 · p. 21

4.1 Hierarchy of nested measurement functions

• An index tree of depth q is a finite set of q-tuples whose entries are chosen from positive-integer ranges that may depend on earlier entries.The notation i(k) denotes a tuple's first k entries, and I(k) collects those prefixes.Source: §4.1 · p. 22
• Construction 4.2 specifies scales 0 < Hq < ··· < H1 < 1 and nested partitions of Ω into Lipschitz convex subsets, each of diameter at most Hk and containing a ball of diameter δHk.Source: §4.1 · p. 22
• Nesting of the domain decomposition implies nesting of the measurement functions, with the stated interlevel entries zero for indices that are not children.The hierarchy associates each index at level k with its descendants at level k+1.Source: §4.1 · p. 22
• The measurement functions may be assumed to have constant norms across indices; otherwise, rescale them and the corresponding entries of π(k,k+1).The text notes that constants may then depend on the maximum domain measure.Source: §4.1 · p. 22

4.2 Hierarchy of nested gamblets and multiresolution approximations

• The mixed strategy for Player II replaces the source term g with a centered Gaussian field of covariance L = -div(a∇).Source: §4.2 · p. 22
• At each level k, u(k) is the best energy-norm approximation to the solution in V(k).Source: §4.2 · p. 23
• The spaces V(k) are nested, with V(k) = span{L^-1φ_i^(k) : i ∈ I(k)}.The theorem states V(k) ⊂ V(k+1) for k = 1, …, q.Source: §4.2 · p. 23

Nested games and martingales

• The coarse-to-fine sequence of game-based solution approximations forms a martingale with independent increments.The measurement σ-algebras form a filtration; the martingale property implies the Gaussian increments are independent.Source: §4.3 · pp. 23–24
• Martingale results can support fluctuation bounds and stopping strategies for numerical simulations.The paper identifies concentration inequalities and optimal stopping based on losses combining computation costs and imperfect decisions.Source: §4.3 · p. 24

Multiresolution decomposition

• The solution decomposes orthogonally into its coarsest approximation and successive increments, each increment solving the PDE in its corresponding detail space.At the finest level, u(q+1) = u.Source: §4.3 · p. 24

4.4 Interpolation and restriction matrices/operators

• Nested spaces V(k) imply the existence of a matrix R(k,k+1) connecting adjacent scales.The passage gives its dimensions as I(k) × I(k+1).Source: §4.4 · p. 25
• R(k,k+1) is the restriction matrix, and its transpose R(k+1,k) is the interpolation/prolongation matrix.Source: §4.4 · p. 25
• Theorem 4.7 characterizes R(k,k+1) as Player II’s best bet on the value of the relevant quantity.The supplied passage fragments do not include the full formula or specify that quantity.Source: §4.4 · p. 25

• Θ(k) is a symmetric positive-definite covariance matrix, and its inverse is the stiffness matrix A(k) of the elements ψ(k).Source: §4.5 · p. 25
• Theorem 4.8 enables hierarchical computation of A(k) from A(k+1) and characterizes R(k+1,k)b as a unique minimizer.Source: §4.5 · pp. 25–26
• The restriction and interpolation maps satisfy inverse-composition identities, while Θ(k) is obtained by projecting Θ(k+1) through π(k,k+1).Source: §4.5 · p. 26

4.6 Multiresolution gamblets

• The elements χ^(k)_i form a basis of W(k), whose members are characterized by having zero restriction to the preceding level.The basis claim follows from the kernel condition and the equality |J(k)| = |I(k)| − |I(k−1)|.Source: §4.6 · p. 26
• To enable fast multiplication by W(k) and its transpose, the construction sets entries to zero across different parent indices and uses m_s−1 independent vectors orthogonal to the all-ones vector within each block.The choice of W(k) is not unique.Source: §4.6 · p. 27
• Construction 4.13 has stated construction complexity |I(k−1)| × m_s^2 and gives W(k)W(k),T equal to the identity J(k).The identity is of size J(k) × J(k).Source: §4.6 · p. 28

4.7 Multiresolution operator inversion

• At each scale k, B(k) is the stiffness matrix, and the text states that it is positive, symmetric, and definite, with inverse B(k)^−1.The passage introduces B(k) as a J(k) × J(k) matrix.Source: §4.7 · p. 28
• The right-hand-side vector g(k) is defined from g on I(k) and can be computed iteratively.The supplied passage gives no further details of the iteration.Source: §4.7 · p. 29
• The solution of (1.1) at any scale can be computed by solving the decoupled linear systems (4.26) and (4.27).The passage identifies this as a direct consequence of Theorem 4.6.Source: §4.7 · p. 29
• For k = 2, ..., q, Theorem 4.15 relates consecutive finite element solutions through an identity for u(k) − u(k−1).Here u(k) is the finite element solution of (1.1) in V(k); the supplied passage truncates the identity after the equals sign.Source: §4.7 · p. 29

• The basis elements yield an orthogonal multiresolution decomposition, described as a generalization of wavelets for the PDE solution space H^1_0(Ω).Source: §4.8 · p. 29
• The decomposition of H^1_0(Ω) induces a subscale decomposition of −div(a∇) into layered subbands of increasing frequencies.Source: §4.8 · pp. 29–30
• If H_(k−1)/H_k is uniformly bounded, the condition number of −div(a∇) restricted to each W^(k) is uniformly bounded.A geometric sequence of scales is given as an example.Source: §4.8 · p. 30
• Under Construction 4.13, Cond(W^(k)W^(k),T) = 1; under Construction 4.11, it is at most max_j 3(m_j−1)m_j/2.Source: §4.8 · pp. 30, 33

4.9 Well conditioned relaxation across subscales

• For geometric scales or uniformly bounded Hk−1/Hk, the systems used to calculate gamblets χ(k) have uniformly bounded condition numbers.Source: §4.9 · p. 33
• The subband solutions u(1) and u(k+1) − u(k) also have condition numbers bounded independently of mesh size or resolution and coefficient regularity.Source: §4.9 · p. 33
• These systems can therefore be solved efficiently using iterative methods such as Conjugate Gradient (CG).Source: §4.9 · p. 33

4.10 Hierarchical localization and error propagation across scales

• Localization is needed to truncate the dense matrices B(k) and A(k) and achieve near-linear complexity; the error control relies on gamblet decay and uniform condition-number bounds.The subsection studies how localization and truncation errors propagate across scales to bound localization-domain sizes.Source: §4.10 · pp. 33–34
• Theorem 4.24 bounds solution error and localized-system conditioning in terms of localization errors E(k) and E(k, χ).It gives E(k, χ) ≤ C H^-d/2 E(k); under stated bounds on E(k, χ) and E(1), it controls Cond(B(k),loc), Cond(A(1),loc), and energy-norm solution errors.Source: §4.10 · p. 36
• Theorem 4.26 uses one-step propagation estimates and finest-scale decay to bound localization error E(k) across the hierarchy.Its proof obtains the bound by induction over scales.Source: §4.10 · p. 37
• With localization radii meeting Theorem 4.27’s lower-bound condition, localized solutions achieve ε-level energy accuracy while localized systems retain bounded condition numbers.The theorem also bounds the error in each scale increment by ε/(2k^2) times the source H^-1 norm, and bounds the difference between the full and localized level-k solutions by C(H^k + ε) times the source L2 norm.Source: §4.10 · p. 37

5.1 The initialisation of the algorithm

• Initialization builds a nested index tree with geometrically spaced resolutions to bound condition numbers, while the method itself requires only specified basis elements and is meshless.The finite-element presentation assumes piecewise-constant coefficients on the fine mesh and targets solution accuracy ϵ.Source: §5.1 · pp. 38–39

5.2 Exact gamblet transform and multiresolution operator inversion

• Algorithm 1 computes gamblets from fine to coarse scales and decomposes the solution into subbands; its nested systems have uniformly bounded condition numbers, unlike the original system.The approximations u(k) have error proportional to Hk for k<q, and zero error at the fine-grid level; the gamblet basis also permits energy-norm compression by zeroing small coefficients.Source: §5.2 · pp. 40, 42, 43, 44

5.3 Fast Gamblet transform/solve

• Algorithm 2 targets near-linear solution and transform costs by localizing or truncating Algorithm 1’s dense systems, with accuracy controlled by exponential decay.The truncation levels depend a priori on the desired accuracy ϵ.Source: §5.2 · p. 40, §5.3 · p. 45

5.4 Complexity vs accuracy of Algorithm 2 and choice of the localization radii ρk

• For direct-solver use, localization radii determine complexity and accuracy, and the algorithm is parallelizable across both space and subscales.The discrete analysis carries over the localization error bounds; it notes inverse decay for symmetric well-conditioned banded matrices in the first gamblet-construction step.Source: §5.4 · pp. 47–48
• With suitable localization radii, Algorithm 2 controls approximation error and keeps its localized system condition numbers uniformly bounded.Theorem 5.5 bounds the localized solution error and states that the algorithm achieves accuracy ϵ in a direct solve when the radii are chosen appropriately.Source: §5.4 · p. 48
• Theorem 5.6 permits inexact solves of Algorithm 2's localized systems, provided each solve meets its specified accuracy requirement.The theorem preserves Theorem 5.5's results under accuracy conditions for the localized systems and multiscale solves.Source: §5.4 · pp. 48–49
• Sparse, banded mass and stiffness matrices can be computed in O(N) time, while localization makes subsequent work depend on local neighborhoods.The passage gives O(N) for matrix computation and describes localized solve costs in terms of the radii and the number of coarse-scale indices.Source: §5.4 · pp. 49–50
• For repeated right-hand sides, gamblets and stiffness matrices are reused, so each subsequent solve costs the multiscale solve rather than recomputing them.The passage identifies this cost with Algorithm 2, Line 9.Source: §5.4 · p. 50
• For Lipschitz-continuous g, a coarse approximation needs O(H^-kd) evaluations and yields coefficient-computation cost independent of N.Representing the approximation on the full fine mesh still incurs a cost stated as O(N(ln 1/ϵ)^d), while restricting output to a subdomain lowers that cost.Source: §5.4 · p. 51
• For L2 approximations with Lipschitz data, coarse coefficients suffice; periodic or ergodic coefficients with mixing length Hp allow gamblets to be computed on representative cells.The cell-based reduction applies when the scale matches the coefficient period or mixing length, and the passage gives a first-solve cost for the resulting approximation.Source: §5.4 · pp. 51–52