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

arXiv:1503.03467v5 [math.NA] 10 Feb 2017

Houman Owhadi∗
February 13, 2017

Abstract
We introduce a near-linear complexity (geometric and meshless/algebraic) multigrid/multiresolution method for PDEs with rough (L∞ ) 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 H01 (Ω) (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 H01 (Ω) (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

1.1

Scientific discovery as a decision theory problem

The process of scientific discovery is oftentimes based on intuition, trial and error and
plain guesswork. This paper is motivated by the question of the existence of a rational
∗

California Institute of Technology, Computing & Mathematical Sciences , MC 9-94 Pasadena, CA
91125, owhadi@caltech.edu

1

decision framework that could be used to facilitate/guide this process, or turn it, to some
degree, into an algorithm. In exploring this question, we will consider the problem of
finding a method for solving (up to a pre-specified level of accuracy) PDEs with rough
(L∞ ) coefficients as fast as possible with the following prototypical PDE (and its possible
discretization over a fine mesh) as an example
(

− div a(x)∇u(x) = g(x) x ∈ Ω; g ∈ L2 (Ω), or g ∈ H −1 (Ω)
(1.1)
u = 0 on ∂Ω,
where Ω is a bounded subset of Rd (of arbitrary dimension d ∈ N∗ ) with piecewise
Lipschitz boundary, a is a symmetric, uniformly elliptic d × d matrix with entries in
L∞ (Ω) and such that for all x ∈ Ω and l ∈ Rd ,
λmin (a)|l|2 ≤ lT a(x)l ≤ λmax (a)|l|2 .

(1.2)

Although multigrid methods [41, 16, 50, 51, 101] are now well known as the fastest
for solving elliptic boundary-problems and have successfully been generalized to other
types of PDEs and computational problems [123], their convergence rate can be severely
affected by the lack of regularity of the coefficients [37, 114]. Furthermore, although significant progress has been achieved in the development of multigrid methods that are, to
some degree, robust with respect to meshsize and lack of smoothness (we refer in particular to algebraic multigrid [94], multilevel finite element splitting [124], hierarchical basis
multigrid [9, 24], multilevel preconditioning [106], stabilized hierarchical basis methods
[107, 109, 110], energy minimization [65, 114, 122, 121, 108] and homogenization based
methods [37, 34]), the design of multigrid methods that are provably robust with respect
to rough (L∞ ) coefficients has remained an open problem of practical importance [17].
Alternative hierarchical strategies for the resolution of (1.1) are (1) wavelet based
methods [18, 15, 3, 29, 38] (2) the Fast Multipole Method [49] and (3) Hierarchical matrices [52, 11]. Although methods based on (classical) wavelets achieve a multiresolution
compression of the solution space of (1.1) in L2 and although approximate wavelets
and approximate L2 projections can stabilize hierarchical basis methods [109, 110], their
applications to (1.1) are limited by the facts that (a) the underlying wavelets can perform arbitrarily badly [7] in their H01 (Ω) approximation of the solution space and (b)
the operator (1.1) does not preserve the orthogonality between subscales/subbands with
classical wavelets. The Fast Multipole Method and hierarchical matrices exploit the
property that sub-matrices of the inverse discrete operator are low rank away from the
diagonal. This low rank property can be rigorously proven for (1.1) (based on the approximation of its Green’s function by sums of products of harmonic functions [10]) and
leads to provable convergence (with rough coefficients), up to the pre-specified level of
accuracy  in L2 -norm, in O(N ln6 N ln2d+2 1 ) operations ([10] and [11, Thm. 2.33 and
Thm. 4.28]). Can the problem of finding a fast solver for (1.1) be, to some degree,
reformulated as an Uncertainty Quantification/Decision Theory problem that could, to
some degree, be solved as such in an automated fashion? Can discovery be computed?
Although these questions may seem unorthodox their answer appears to be positive: this
2

paper shows that this reformulation is possible and leads to a multigrid/multiresolution
method/algorithm solving (1.1), up to the pre-specified level of accuracy  in H 1 -norm
(i.e. finding uapp such that ku − uapp kH01 (Ω) ≤ kgkH −1 (Ω) for an arbitrary g decomposed

over N degrees of freedom), in O N ln3d max( 1 , N 1/d ) operations (for  ∼ N −1/d ,
the hierarchical matrix method achieves -accuracy in L2 norm in O(N ln2d+8 N ) operations and the proposed multiresolution method achieves -accuracy in H 1 norm in
O(N ln3d N ) operations). For subsequent solves (i.e. if (1.1) needs to be solved for
1
more than one g) then the proposed multiresolution method achieves accuracy  ≈ N − d
in H 1 -norm in O(N lnd+1 N ) operations (we refer to Subsection 5.4 and in particular
to Table 1 for a detailed complexity analysis of the proposed method, which can also
achieve sublinear complexity if one only requires L2 -approximations).
The core mechanism supporting the complexity of the method presented here is the
fast decomposition of H01 (Ω) into a direct sum of linear subspaces that are orthogonal
(or near-orthogonal) with respect to the energy scalar product and over which (1.1)
has uniformly bounded condition numbers. It is, to some degree, surprising that this
decomposition can be achieved in near linear complexity and not in the complexity of
an eigenspace decomposition. Naturally [86], this decomposition can be applied to the
fast simulation of the wave and parabolic equations associated to (1.1) or to its fast
diagonalization.
The essential step behind the automation of the discovery/design of scalable numerical solvers is the observation that fast computation requires repeated computation with
partial information (and limited resources) over hierarchies of levels of complexity and
the reformulation of this process as that of playing underlying hierarchies of adversarial
information games [111, 112].
Although the problem of finding a fast solver for (1.1) may appear disconnected
from that of finding statistical estimators or making decisions from data sampled from an
underlying unknown probability distribution, the proposed game theoretic reformulation
is, to some degree, analogous to the one developed in Wald’s Decision Theory [113],
evidently influenced by Von Neumann’s Game Theory [111, 112] (the generalization of
worst case Uncertainty Quantification analysis [83] to sample data/model uncertainty
requires an analogous game theoretic formulation [80], see also [79] for how the underlying
calculus could be used to guide the discovery of new Selberg identities). We also refer
to subsection 1.3 for a review of the correspondence between statistical inference and
numerical approximation.

1.2

Outline of the paper

The essential difficulty in generalizing the multigrid concept to PDEs with rough coefficients lies in the fact that the interpolation (downscaling) and restriction (upscaling) operators are, a priori, unknown. Indeed, in this situation, piecewise linear finite-elements
can perform arbitrarily badly [7] and the design of the interpolation operator requires
the identification of accurate basis elements adapted to the microstructure a(x).
This identification problem has also been the essential difficulty in numerical homoge-

3

nization [117, 6, 4, 18, 59, 33, 84, 17]. Although inspired by classical homogenization ideas
and concepts (such as oscillating test functions [68, 36, 35], cell problems/correctors and
effective coefficients [13, 90, 1, 72, 39, 46], harmonic coordinates [61, 6, 4, 76, 12, 2, 84],
compactness by compensation [100, 45, 67, 14]) an essential goal of numerical homogenization has been the numerical approximation of the solution space of (1.1) with arbitrary rough coefficients [84], i.e., in particular, without the assumptions found in classical
homogenization, such as scale separation, ergodicity at fine scales and -sequences of operators (otherwise the resulting method could lack robustness to rough coefficients, even
under the assumption that coefficients are stationary [8]). Furthermore, to envisage applications to multigrid methods, the computation of these basis functions must also be
provably localized [5, 85, 64, 48, 87, 58] and compatible with nesting strategies [87]. In
[77], it has been shown that this process of identification (of accurate basis elements for
numerical homogenization), could, in principle, be guided through its reformulation as a
Bayesian Inference problem in which the source term g in (1.1) is replaced by noise ξ and
one tries to estimate the value of the solution at a given point based on a finite number
of observations. In particular it was found that Rough Polyharmonic Splines [87] and
Polyharmonic Splines [54, 30, 31, 32] can be re-discovered as solutions of Gaussian filtering problems. This paper is inspired by the suggestion that this link between numerical
homogenization and Bayesian Inference (and the link between Numerical Quadrature
and Bayesian Inference [92, 27, 97, 74, 75]) are not coincidences but particular instances
of mixed strategies for underlying information games and that optimal or near optimal
methods could be obtained by identifying such games and their optimal strategies.
The process of identification of these games starts with the (Information Based Complexity [119]) notion that computation can only be done with partial information. For
instance, since the operator (1.1) is infinite dimensional, one cannot directly compute
with u ∈ H01 (Ω) but only with finite-dimensional features ofRu. An example
R of such finitedimensional features is the m-dimensional vector um := ( Ω uφ1 , . . . , Ω uφm ) obtained
by integrating the solution u of (1.1) against m test/measurement functions φi ∈ L2 (Ω).
However to achieve an accurate approximation of u through computation with um one
must fill the information gap between um and u (i.e. construct an interpolation operator
giving u as a function of um ). We will, therefore, reformulate the identification of this
interpolation operator as a non-cooperative (min max) game where Player I chooses the
source term g (1.1) in an admissible set/class (e.g. the unit ball of L2 (Ω)) and Player II
is shown um and must approximate u from these incomplete measurements. Using the
energy norm
Z
∇uT (x)a(x)∇u(x) dx,

kuk2a :=

(1.3)

Ω

to quantify the accuracy of the recovery and calling u∗ Player I’s bet (on the value of
u), the objective of Player I is to maximize the approximation error ku − u∗ ka , while
the objective of Player II is to minimize it. A remarkable result from Game Theory
(as developed by Von Neumann [111], Von Neumann and Morgenstern [112] and Nash
[70]) is that optimal strategies for deterministic zero sum finite games are mixed (i.e.
randomized) strategies. Although the information game described above is zero sum,
4

it is not finite. Nevertheless, as in Wald’s Decision Theory [113], under sufficient regularity conditions it can be made compact and therefore approximable by a finite game.
Therefore although the information game described above is purely deterministic (and
has no a priori connection to statistical estimation), under compactness (and continuity
of the loss function), the best strategy for Player I is to play at random by placing a
probability distribution πI on the set of candidates for g (and select g as a sample from
πI ) and the optimal strategy for Player II is to place a probability distribution πII on
the set of candidates for g and approximate the solution of (1.1) by the expectation
of
R
u (under πII used as a prior distribution) conditioned on the measurements Ω uφi .
Although the estimator employed by Player II may be called Bayesian, the game
described here is not (i.e. the choice of Player I might be distinct from that of Player II)
and Player II must solve a min max optimization problem over πI and πII to identify an
optimal prior distribution for the Bayesian estimator (a careful choice of the prior also
appears to be important due to the possible high sensitivity of posterior distributions
[81, 79, 82]). Although solving the min max problem over πI and πII may be one way of
determining the strategy of Player II, it will not be the method employed here. We will
instead analyze the error of Player II’s approximation as a function of Player II’s prior
and the source term g picked by Player I. Furthermore, to preserve the linearity of the
calculations we will restrict Player II’s decision space (the set of possible priors πII ) to
Gaussian priors on the source term g. Since the resulting analysis is independent of the
structure of (1.1) and solely depends on its linearity we will first perform this investigation, in Section 2, in the algebraic framework of linear systems of equations, identify
Player II’s optimal mixed strategy and show that it is characterized by deterministic
optimal recovery and accuracy properties. The mixed strategy identified in Section 2
will then be applied in 3 to the numerical homogenization of (1.1) and the discovery of
interpolation interpolators. In particular, it will be shown that the resulting elementary
gambles form a set of deterministic basis functions (gamblets) characterized by (1) optimal recovery and accuracy properties (2) exponential decay (enabling their localized
computation) (3) robustness to high contrast.
To compute fast, the game presented above must not be limited to filling the information gap between um ∈ Rm and u ∈ H01 (Ω). This game must be played (and repeated)
over hierarchies of levels of complexity (e.g. one must fill information gaps between R4
and R16 , then R16 and R64 , etc...). We will therefore, in Section 4, consider the (hierarchical) game where Player I chooses the r.h.s of (1.1) and Player II must (iteratively)
gamble on the value of its solution based on a hierarchy of nested measurements of u
(from coarse to fine measurements). Under Player II’s mixed strategy (identified in Section 2 and used in Section 3), the resulting sequence of multi-resolution approximations
forms a martingale. Conditioning and the independence of martingale increments lead
to the hierarchy of nested interpolation operators and to the multiresolution orthogonal
decomposition of (1.1) into independent linear systems of uniformly bounded condition
numbers. The resulting elementary gambles (gamblets) (1) form a hierarchy of nested
basis functions leading to the orthogonal decomposition (in the scalar product of the
energy norm) of H01 (Ω) (2) enable the sparse compression of the solution space of (1.1)
5

(3) can be computed and stored in near-linear complexity by solving a nesting of linear
systems with uniformly bounded condition numbers (4) enable the computation of the
solution of (1.1) (or its hyperbolic or parabolic analogues) in near-linear complexity.
The implementation and complexity of the algorithm are discussed in Section 5 with
numerical illustrations.

1.3

On the correspondence between statistical inference and numerical
approximation

As exposed by Diaconis [27], the investigation of the correspondence between statistical inference and numerical approximation can be traced back to Poincaré’s course in
Probability Theory [92]. It is useful to recall Diaconis’ compelling example [27] as an
illustration of this conection. Let f : [0, 1] → R beRa given function and assume that we
1
are interested in the numerical approximation of 0 f (t) dt. The Bayesian approach to
this quadrature problem is to (1) Put a prior (probability distribution) on continuous
functions C[0, 1] (2) Calculate f at x1 , Rx2 , . . . , xn (to obtain the data (f (x1 ), . . . , f (xn )))
1
(3) Compute a posterior (4) Estimate 0 f (t) dt by the Bayes rule. If the prior on C[0, 1]
is that of a Brownian
Motion (i.e. f (t) = Bt where Bt is a Brownian motion and B0 is

normal), then E f (x) f (x1 ), . . . , f (xn ) is the piecewise linear interpolation of f between
the points x1 , . . . , xn and one re-discovers the trapezoidal quadrature rule.
R t If the prior
on C[0, 1] is that
of
the
first
integral
of
a
Brownian
Motion
(i.e.
f
(t)
∼
0 Bs ds) then


the posterior E f (x) f (x1 ), . . . , f (xn ) is the cubic spline interpolant and integrating k
times yields splines of order 2k + 1.
Subsequent to Poincaré’s early discovery [92], Sul’din [102] and (in particular) Larkin
[62] initiated the systematic investigation of the correspondence between conditioning
Gaussian measures/processes and numerical approximation. As noted by Larkin [62], despite Sard’s introduction of probabilistic concepts in the theory of linear approximation
[95], and Kimeldorf and Wahba’s exposition [60] of the correspondence between Bayesian
estimation and spline smoothing/interpolation, the application of probabilistic concepts
and techniques to numerical integration/approximation “attracted little attention among
numerical analysts” (perhaps due to the counterintuitive nature of the process of randomizing a known function). However, a natural framework for understanding this process
of randomization can be found in the pioneering works of Woźniakowski [118], Packel
[88], and Traub, Wasilkowski, and Woźniakowski [104] on Information Based Complexity [71, 119], the branch of computational complexity that studies the complexity of
approximating continuous mathematical operations with discrete and finite ones up a to
specified level of accuracy. Indeed the concept that numerical implementation requires
computation with partial information and limited resources emerges naturally from Information Based Complexity, where it is also augmented by concepts of contaminated
and priced information associated with, for example, truncation errors and the cost of
numerical operations. In this framework, the performance of an algorithm operating on
incomplete information can be analysed in the usual worst case setting or the average
case (randomized) setting [93, 73] with respect to the missing information. Although
the measure of probability (on the solution space) employed in the average case setting
6

may be arbitrary, as observed by Packel [88], if that measure is chosen carefully (as the
solution of a game theoretic problem) then the average case setting can be interpreted
as lifting a (worst case) min max problem (where saddle points of pure strategies do
not, in general, exist) to a min max problem over mixed (randomized) strategies (where
saddle points do exist [111, 112]). As exposed by Diaconis [27] (see also Shaw [97]) the
randomized setting also establishes a correspondence between Numerical Analysis and
Bayesian Inference providing a natural framework for the statistical description of numerical errors (in which confidence intervals can be derived from posterior distributions).
Furthermore [89, 27], classical min max numerical quadrature rules can be formulated
as solutions of Bayesian inference problems with carefully chosen priors [27] and, as
shown by O’Hagan [74, 75], this correspondence can be exploited to discover new and
useful numerical quadratures rules. As envisioned by Skilling [99], by placing a (carefully chosen) probability distribution on the solution space of an ODE and conditioning
on quadrature points, one obtains a posterior distribution on the solution whose mean
may coincide with classical numerical integrators such as Runge-Kutta methods [96]. As
shown in [23] the statistical approach is particularly well suited for chaotic dynamical
systems for which deterministic worst case error bounds may provide little information.
While in [99, 96, 23] the probability distribution is directly placed on the solutions space,
for PDEs [77] argues that the prior distribution must be placed on source terms (or on
the image space of an integro-differential operator) and propagated/filtered through the
inverse operator to reflect the structure of the solution space. In particular [77] shows
that this process of filtering noise with the inverse operator, when combined with conditioning, produces accurate finite-element basis functions for the solution space whose
deterministic worst case errors can be bounded by standard deviation errors using the
reproducing kernel structure of the covariance function of the filtered Gaussian field. As
already witnessed in [23, 96, 77, 56, 55, 19, 25], it is natural to expect that the possibilities offered by combining numerical uncertainties/errors with model uncertainties/errors
in a unified framework will stimulate a resurgence of the statistical inference approach
to numerical analysis.

2
2.1

Linear Algebra with incomplete information
The recovery problem

The problem of identifying interpolation operators for (1.1) is equivalent (after discretization or in the algebraic setting) to that of recovering or approximating the solution of
a linear system of equations from an incomplete set of measurements (coarse variables)
given known norm constraints on the image of the solution.
Let n ≥ 2 and A be a known real invertible n × n matrix. Let b be an unknown
element of Rn . Our purpose is to approximate the solution x of
Ax = b

7

(2.1)

based on the information that (1) x solves
Φx = y,

(2.2)

where Φ (the measurement matrix) is a known, rank m, m × n real matrix such that
m < n and y (the measurement vector) is a known vector of Rm , and (2) the norm
bT T −1 b of b is known or bounded by a known constant (e.g., bT T −1 b ≤ 1), where T −1
is a known positive definite n × n matrix (with T −1 being the identity matrix as a
prototypical example). Observe that since m < n, the measurements (2.2) are, a priori,
not sufficient to recover the exact value x.
As described in Section 1, by formulating this recovery problem as a (non-cooperative)
information game (where Player I chooses b and Player II chooses an approximation x∗
of x based on the observation Φx), one (Player II) is naturally lead to search for mixed
strategy in the Bayesian class by placing a prior distribution on b. The purpose of this
section is to analyze the resulting approximation error and select the prior distribution
accordingly. To preserve the linearity (i.e. simplicity and computational efficiency) of
calculations we will restrict Player II’s decision space to Gaussian priors.

2.2

Player I’s mixed strategy

We will therefore, in the first step of the analysis, replace b in (2.1) by ξ, a centered
Gaussian vector of Rn with covariance matrix Q (which may be distinct from T ) and
consider the following stochastic linear system
AX = ξ .

(2.3)

The Bayesian answer (a mixed strategy for Player II) to the recovery problem of Section
2 is to approximate x by the conditional expectation E[X|ΦX = y].
Theorem 2.1. The solution X of (2.3) is a centered Gaussian vector of Rn with covariance matrix
K = A−1 Q(A−1 )T .
(2.4)
Furthermore, X conditioned on the value ΦX = y is a Gaussian vector of Rn with mean
E[X|ΦX = y] = Ψy, and of covariance matrix K Φ , where Ψ is the n × m matrix
Ψ := KΦT (ΦKΦT )−1 ,

(2.5)

and K Φ is the rank n − m positive n × n symmetric matrix defined by K Φ := K − ΨΦK.
Proof. (2.4) simply follows from X = A−1 ξ. Since X is a Gaussian vector,
 E[X|ΦX 2=
y] = Ψy where Ψ is a n × m matrix minimizing
the
mean
squared
error
E
|X − M ΦX|


2
over all n×m matrices M . We have E |X −M ΦX| = Trace[K]+Trace[M ΦKΦT M T ]−
2 Trace[ΦKM ] whose minimum is achieved for M = Ψ as defined by (2.5). The covariance matrix
obtained by observing that for v ∈ Rn ,
 T of X Tgiven ΦX
 = Ty is then
T
Φ
2
T
v K v = E |v X − v ΨΦX| = v Kv − v ΨΦKv.
8

2.3

Variational/optimal recovery properties and approximation error

For a n × n symmetric positive definite matrix M let ·, ·
Rn

Rn ,

M

be the (scalar) product

uT M v

1

2
the
on
defined by: for u, v ∈
u, v M :=
and write kvkM := v, v M
corresponding norm. When M is the identity matrix then we write u, v and kvk the
corresponding scalar product and norm. For a linear subspace V of Rn we write PV,M
for the orthogonal projections onto V with respect to the scalar product ·, · M . For a
(possibly rectangular) matrix B we write Im(B) the image (range) of B and Ker(B) the
null space of B. For an integer n, let In be the n × n identity matrix.

Theorem 2.2. For w ∈ Rm , Ψw is the unique minimizer of the following quadratic
problem
(
Minimize
v, v K −1
(2.6)
Subject to Φv = w and v ∈ Rn .
In particular, v = Ψy, the Bayesian approximation of the solution of (2.1), is the unique
minimizer of kAvkQ−1 under the measurement constraints Φv = y. Furthermore, it also
holds true that (1) ΦΨ = Im (2) Im(Ψ) is the orthogonal complement of Ker(Φ) with
respect to the product ·, · K −1 and (3) ΨΦ = PIm(KΦT ),K −1 and In − ΨΦ = PKer(Φ),K −1 .
Proof. First observe that (2.5) implies that ΦΨ = Im where Im is the identity m×m matrix. Therefore Φ(Ψw) = w. Note that (2.5) implies that for all z ∈ Rm , Ψz, v K −1 =
−1
z T ΦKΦT
Φv. Therefore if v ∈ Ker(Φ) then Ψz, v K −1 = 0 for all z ∈ Rm . Conversely if Ψz, v K −1 = 0 for all z ∈ Rm then v must belong to Ker(Φ). Since the
dimension of Im(Ψ) is m and that of Ker(Φ) is n − m we conclude that Im(Ψ) is the
orthogonal complement Ker(Φ) with respect to the product ·, · K −1 and in particular,
Ψw, v

K −1

= 0,

∀w ∈ Rm and ∀v ∈ Rn such that Φv = 0 .

(2.7)

Let w ∈ Rm and v ∈ Rn such that Φv = w. Since Ψw − v ∈ Ker(Φ), it follows from (2.7)
that v, v K −1 = Ψw, Ψw K −1 + v − Ψw, v − Ψw K −1 . Therefore Ψw is the unique
minimizer of v, v K −1 over all v ∈ Rn such that Φv = w. Now consider f ∈ Rn , since
Im(Ψ) = Im(KΦT ) and Im(Ψ) is the orthogonal complement of Ker(Φ) with respect to
the product ·, · K −1 , there exists a unique z ∈ Rm and a unique g ∈ Ker(Φ) such that
f = KΦT z + g. Since ΨΦ = KΦT (ΦKΦT )−1 Φ, it follows that ΨΦf = KΦT z and (In −
ΨΦ)f = g. We conclude by observing that g = PKer(Φ),K −1 f .
Theorem 2.3. For v ∈ Rn , w∗ = Φv is the unique minimizer of kv − ΨwkK −1 over
all w ∈ Rm . In particular, kv − ΨΦvkK −1 = minz∈Rm kv − KΦT zkK −1 and if x is the
solution of the original equation (2.1), then kx − ΨykK −1 = minw∈Rm kx − ΨwkK −1 =
minz∈Rm kx − KΦT zkK −1 .
Proof. The proof follows by observing that v −ΨΦv belongs to the null space of Φ which,
from Theorem 2.2, is the orthogonal complement of the image of Ψ with respect to the
scalar product defining the norm k · kK −1 . Observe also that the image of Ψ is equal to
that of KΦT .
9

Remark 2.4. Observe that, from Theorem 2.2, v − ΨΦv spans the null space of Φ, and
2
2
kvk2K −1 = v − ΨΦv K −1 + ΨΦv K −1 . Therefore if D is a symmetric positive definite
n × n matrix then supv∈Rn v − ΨΦv D /kvkK −1 = supv∈Rn , Φv=0 kvkD /kvkK −1 . In particular, if x is the solution of (2.1) and y the vector in (2.2), then x − Ψy D /kbkQ−1 ≤
supv∈Rn , Φv=0 kvkD /kvkK −1 and the right hand side is the smallest constant for which
the inequality holds (for all b).
Remark 2.5. A simple calculation (based on the reproducing Kernel property v, K·,i K −1 =
vi ) shows that if x is the solution of (2.1) and y the vector in (2.2), then
1
Φ 2 kbk
(x − Ψy)i ≤ Ki,i
Q−1 , i.e. the variance of the ith entry of the solution of the
stochastic system (2.3) conditioned on ΦX = y, controls the accuracy of the approximation of the ith entry of the solution of the deterministic system (2.1). In that sense, the
role of K Φ is analogous to that of the power function in radial basis function interpolation
[116, 40] and that of the Kriging function [120] in geostatistics [69].

2.4

Energy norm estimates and selection of the prior

We will from now on assume that A is symmetric positive definite. Observe that in this
situation the energy norm k · kA is of practical significance for quantifying the approxi1
mation error and Theorem 2.3 leads to the estimate kx − ΨykK −1 = minz∈Rm kQ− 2 b −
1
1
1
Q− 2 A 2 K 2 ΦT zk which simplifies to the energy norm estimate expressed by Corollary
(2.6) under the choice Q = A (note that K −1 = A under that choice).
Corollary 2.6. If A is symmetric positive definite and Q = A, then for v ∈ Rn ,
kv − ΨΦvkA = minz∈Rm kv − A−1 ΦT zkA . Therefore, if x is the solution of (2.1) and y
the vector in (2.2), then kx − ΨykA = minw∈Rm kx − ΨwkA = minz∈Rm kx − A−1 ΦT zkA .
In particular
1
1
kx − ΨykA = minm kA− 2 b − A− 2 ΦT zk .
(2.8)
z∈R

Remark 2.7. Therefore, according to Corollary 2.6, if Q = A, then Ψy is the Galerkin
approximation of x, i.e. the best approximation of x in k · kA -norm in the image of
Ψ (which is equal to the image of A−1 ΦT ). This is interesting because Ψy is obtained
without the prior knowledge of b.
Corollary 2.6 and Remark 2.7 motivate us to select Q = A as the covariance matrix
of the Gaussian prior distribution (mixed strategy of Player II).

2.5

Impact and selection of the measurement matrix Φ

It is natural to wonder how good this recovery strategy is (under the choice Q = A)
compared to the best possible function of y and how the approximation error is impacted
by the measurement matrix Φ. If the energy norm is used to quantify accuracy, then
the recovery problem can be expressed as finding the function θ of the measurements
y minimizing the (worst case) approximation error inf θ supkbk≤1 kx − θ(y)kA /kbk with
10

x = A−1 b and y = ΦA−1 b. Writing 0 < λ1 (A) ≤ · · · ≤ λn (A), the eigenvalues of A
in increasing order, and a1 , . . . , an , the corresponding eigenvectors, it is easy to obtain
that (1) the best choice for Φ would correspond to measuring the p
projection of x on
span{a1 , . . . , am } and would lead to the worst approximation error 1/ λm+1 and (2) the
worst choice would correspond to measuring the projection
√ of x on a subspace orthogonal
to a1 and would lead to the worst approximation
error
1/
λ1 . Under the decision Q = A
p
the minimal value of (2.8) is also 1/ √λm+1 and achieved for Im(ΦT ) = span{a1 , . . . , am }
and the maximal value of (2.8) is 1/ λ1 and achieved when Im(ΦT ) is orthogonal to a1 .
The following theorem, which is a direct application of (2.8) and the estimate derived in
[53, p. 10] (see also [66]), shows that, the subset of measurement matrices that are not
nearly optimal is of small measure if the rows of ΦT are sampled independently on the
unit sphere of Rn .
Theorem 2.8. If Φ is a n × m matrix with i.i.d. N (0, 1) (Gaussian) entries, Q = A,
x is the solution of the original equation (2.1), p
and 2 ≤ p then with probability at least
√
√
1 − 3p−p , kx − ΨykA /kbk ≤ (1 + 9 m + p n)/ λm+1 .
Although the randomization of the measurement matrix [42, 63, 43, 78] can be an
efficient strategy in compressed sensing [105, 21, 20, 28, 44, 22] and in Singular Value
Decomposition/Low Rank approximation [53], we will not use this strategy here because
the design of the interpolation operator presents the (added) difficulty of approximating
the eigenvectors associated with the smallest eigenvalues of A rather than those associated with the largest ones. Furthermore, Ψ has to be computed efficiently and the
dependence of the approximation constant in Theorem 2.8 on n and m can be problematic if sharp convergence estimates are to be obtained. We will instead select the
measurement matrix based on the transfer property introduced in [14] and given in a
discrete context in the following theorem.
Theorem 2.9. If A is symmetric positive definite, Q = A and x is the solution of the
original equation (2.1), then for any symmetric positive definite matrix B, we have
r
r
v T Bv
v T Bv
T
infn
min
min kb−ΦT zkB −1 (2.9)
kb−Φ
zk
−1 ≤ kx−ΨykA ≤ sup
B
v∈R
v T Av z∈Rm
v T Av z∈Rm
v∈Rn
Proof. Corollary 2.6 implies that if x is the solution of the original equation (2.1), then
kx − ΨykA = minz∈Rm kb − ΦT zkA−1 . We finish the proof by observing that if A and B
are symmetric positive definite matrices such that α1 B ≤ A ≤ α2 B for some constants
α1 , α2 > 0 then α2−1 B −1 ≤ A−1 ≤ α1−1 B −1 .
Therefore according to Theorem 2.9, once a good measurement matrix Φ has been
identified for a symmetric positive definite matrix B such that α1 B ≤ A, the same
measurement matrix can be used for A at the cost of an increase of the bound on the
−1/2
error by the multiplicative factor α1 . As a prototypical example, one may consider a
(stiffness) matrix A obtained from a finite element discretization of the PDE (1.1) and

11

B may be the stiffness matrix of the finite element discretization of the Laplace Dirichlet
PDE
− ∆u0 (x) = g(x) on Ω with u0 = 0 on ∂Ω,
(2.10)
obtained from the same finite-elements (e.g. piecewise-linear nodal basis functions over
the same fine mesh Th ). Using the energy norm (1.3), Theorem 2.9 and Remark 2.7
imply the following proposition
Proposition 2.10. Let uh (resp. u0h ) be the finite element approximation of the solution
u of (1.1) (resp. the solution u0 of (2.10)) over the finite nodal elements of Th . Let uH
(resp. u0H ) be the finite element approximation of the solution u of (1.1) (resp. the
solution u0 of (2.10)) over linear space spanned by the rows of A−1 ΦT (resp. over the
linear space spanned by the rows of B −1 ΦT ). It holds true that
1
1
p
ku0h − u0H kH01 (Ω) ≤ kuh − uH ka ≤ p
ku0h − u0H kH01 (Ω)
λmax (a)
λmin (a)

(2.11)

Observe that the right hand side of (2.11) does not depend on λmax (a), therefore if
λmin (a) = 1, then the error bound on kuh − uH ka does not depend on the contrast of a
(i.e. λmax (a)/λmin (a)).

3

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

We will now generalize the results and continue the analysis of Section 2 in the continuous case and design the interpolation operator for (1.1) in the context of numerical
homogenization.

3.1

Information Game and Gamblets

As in Section 2 we will identify the interpolation operator (that will be used for the
multigrid algorithm) through a non cooperative game formulation where Player I chooses
the source term g (1.1) and Player II tries
to approximate the solution u of (1.1) based
R
on a finite number of measurements ( Ω uφi )1≤i≤m obtained from linearly independent
test functions φi ∈ L2 (Ω). As in Section 2, this game formulation, motivates the search
for a mixed strategy for Player II that can be expressed by replacing the source term g
with noise ξ. We will therefore consider the following SPDE


(
− div a(x)∇v(x) = ξ(x) x ∈ Ω;
(3.1)
v = 0 on ∂Ω,
where Ω and a are the domain and conductivity of (1.1). As in Section 2, to preserve the
computational efficiency of the interpolation operator we will assume that ξ is a centered
Gaussian field on Ω. The decision space of Player II is therefore the covariance function
of ξ. Write L the differential operator − div(a∇) with zero Dirichlet boundary condition
12

mapping H01 (Ω) onto H −1 (Ω). Motivated by the analysis (Remark 2.7) of Subsection 2.4
(which can be reproduced in the continuous case) we will select the covariance function
1
of
R ξ (Player II’s decision) to be L. Therefore, under that choice, for all
R f ∈ H0 (Ω),2
Ω f (x)ξ(x) dx is a Gaussian random variable with mean 0 and variance Ω f Lf = kf ka
where kf ka is the energy norm of f defined in (1.3). Introducing the scalar product on
H01 (Ω) defined by
Z
v, w

a

(∇v)T a∇w ,

:=

(3.2)

Ω

1
recall that if (e1 , e2 , . . .) is an orthonormal
P∞ basis of (H0 (Ω), k · ka ) diagonalizing L, then ξ
can formally be represented as ξ = i=1 (Lei )Xi (where the Xi are i.i.d. N (0, 1) random
variables) and, therefore, ξ can also
H01 (Ω)
R linear isometryPmapping
P∞ be identified as the
∞
onto a Gaussian space and f = i=1 f, ei a ei onto Ω f (x)ξ(x) dx = i=1 f, ei a Xi .
Observe also that [77], if ξ 0 is White Noise on Ω (i.e. a Gaussian field with covariance
1
function δ(x − y)) then ξ can be represented as ξ = L− 2 ξ 0 . Furthermore [77, Prop. 3.1]
the solution of (3.1) is Gaussian field with covariance function G(x, y) (where G is the
Green’s function of the PDE (1.1), i.e. LG(x, y) = δ(x−y) with
R G(x, y) = 0 for y ∈ ∂Ω).
Let F be the σ-algebra generated by the random variables Ω v(x)φi for i ∈ {1, . . . , m}
(with v solution of (3.1)). We will identify the interpolation basis elements by conditioning the Rsolution of (3.1) on F. Observe that the covariance matrix of the measurement
vector ( Ω v(x)φi )1≤i≤m is the m × m symmetric matrix Θ defined by
Z
Θi,j :=
φi (x)G(x, y)φj (y) dx dy
(3.3)
Ω2

Note thatPfor l ∈ Rm , lT Θl = kwk2a where w is the solution of (1.1) with right hand
side g = m
i=1 li φi . Therefore (since the test functions φi are linearly independent) Θ
is positive definite and we will write Θ−1 its inverse. Write δi,j the Kronecker’s delta
(δi,i = 1 and δi,j = 0 for i 6= j).
Theorem 3.1. Let v be the solution of (3.1). It holds true that
Z
m

 X
E v(x) F =
ψi (x) v(y)φi (y) dy

(3.4)

Ω

i=1

where the functions ψi ∈ H01 (Ω) are defined by
Z
h
i
ψi (x) := E v(x)
v(y)φj (y) dy = δi,j , j ∈ {1, . . . , m}

(3.5)

Ω

and admit the following representation formula
Z
m
X
−1
ψi (x) =
Θi,j
G(x, y)φj (y) dy .
j=1

(3.6)

Ω

Furthermore, the distribution of v conditioned on FPis that of a Gaussian field with mean
(3.4) and covariance function Γ(x, y) = G(x, y) + m
i,j=1 ψi (x)ψj (y)Θi,j
R
R
Pm
Pm
− i=1 ψi (x) Ω G(y, z)φi (z) dz − i=1 ψi (y) Ω G(x, z)φi (z) dz .
13

Proof. The proof is similar to that of [77, Thm. 3.5]. Theh identification of the covari

ance function follows from the expansion of Γ(x, y) = E v(x) − E v(x) F v(y) −

i
E v(y) F . Note that (3.6) proves that ψi ∈ H01 (Ω).
Since, according to (3.5) and the discussion preceding (3.1), each ψi Ris an elementary
gamble (bet) on value of the solution of (1.1) given the information Ω φj u = δi,j for
j = 1, . . . , m we will refer to the basis functions (ψi )1≤i≤m , as gamblets. According to
(3.4), once gamblets have been identified, theyR form a basis for betting on the value of
the solution of (1.1) given the measurements ( Ω φj u)1≤i≤m .

3.2

Optimal recovery properties

Although gamblets admit the representation formula (3.6), we will not use it for their
practical (numerical) computation. Instead we will work with variational properties
inherited from the conditioning of the Gaussian field v. To guide our intuition, note that
since L is the precision function (inverse of the covariance
function) of v, the conditional
R
expectation of v can be identified by minimizing Ω ψLψ given measurements constraints.
This observation motivates us to consider, for i ∈ {1, . . . , m}, the following quadratic
optimization problem
(
Minimize
kψka
R
(3.7)
Subject to ψ ∈ H01 (Ω) and Ω φj ψ = δi,j for j = 1, . . . , m
where kψka is the energy norm of ψ defined in (1.3).
The following theorem shows that (3.7) can be used to identify ψi and that gamblets
are characterized by optimal (variational) recovery properties.
Theorem 3.2. It holds true that (1) The optimization problem
Pm(3.7) admits a unique
m
minimizer ψi defined by (3.5)R and (3.6) (2) For w ∈ R ,
i=1 wi ψi is the unique
minimizer of kψka subject to Ω ψ(x)φj (x) = wj for j ∈ {1, . . . , m} and (3) (using the
scalar product defined in (3.2)) ψi , ψj a = Θ−1
i,j .
P
m
Proof. Let w ∈ Rm
R and ψw = i=1 wi ψi with ψi defined as in (3.6). The definition
of Θ implies that Ω ψw (x)φj (x) = wj for j ∈ {1, . . . , m}. Furthermore we obtain by
R
P
−1
integration by parts that for all ϕ ∈ H01 (Ω), ψw , ϕ a = m
i,j=1 wi Θi,j Ω φj ϕ. Therefore,
R
if ψ ∈ H01 (Ω) is such that Ω ψ(x)φj (x) = wj for j ∈ {1, . . . , m} then ψw , ψ − ψw a = 0
and
kψk2a = kψw k2a + kψ − ψw k2a
(3.8)
which finishes the proof of optimality of ψi and ψw .

3.3

Optimal accuracy of the recovery

Define
∗

u (x) :=

m
X

Z
ψi (x)

u(y)φi (y) dy
Ω

i=1

14

(3.9)

where u is the solution of (1.1) and ψi are the gamblets defined by (3.5) and (3.6).
Note
u∗ corresponds to Player II’s bet on the value of u given the measurements
R
( Ω u(y)φi (y) dy)1≤i≤m . In particular, if v is the solution of (3.1) then
Z
Z


∗
u(y)φi (y) dy
(3.10)
v(y)φi (y) dy =
u (x) = E v(x)
Ω

Ω

For φ ∈ H −1 (Ω) write L−1 φ the solution of (1.1) with g = φ. The following Theorem
shows that u∗ is the best approximation (in energy norm) of u in span{L−1 φi : i ∈
{1, . . . , m}}.
Theorem 3.3. Let u be the solution of (1.1), u∗ defined in (3.9) and (3.10). It holds
true that
ku − u∗ ka =
inf
ku − ψka
(3.11)
ψ∈span{L−1 φi :i∈{1,...,m}}

−1 φ : i ∈ {1, . . . , m}} = span{ψ , . . . , ψ } and (3.11)
Proof. By Theorem 3.1 span{L
i
1
m
R
follows from the fact that Ω (u − u∗ )φj = 0 for all j implies that u − u∗ is orthogonal to
span{ψ1 , . . . , ψm } with respect to the scalar product ·, · a .

3.4

Transfer property and selection of the measurement functions

We will now select the measurement (test) functions φi by extending the result of Proposition 2.10 to the continuous case. For V , a finite dimensional linear subspace of H −1 (Ω),
define
(div a∇)−1 V := span{(div a∇)−1 φ : φ ∈ V }.
(3.12)
where (div a∇)−1 φ is the solution of (1.1) with g = −φ. Similarly define ∆−1 V :=
span{∆−1 φ : φ ∈ V } where ∆−1 φ is the solution of (2.10) with g = −φ.
Proposition 3.4. If u and u0 are the solutions of (1.1) and (2.10) (with the same r.h.s.
g) and V is a finite dimensional linear subspace of H −1 (Ω), then
1
1
p
inf ku0 −vkH01 (Ω) ≤
inf
ku−vka ≤ p
inf ku0 −vkH01 (Ω)
−1 V
−1 V
−1 V
v∈∆
v∈(div
a∇)
v∈∆
λmax (a)
λmin (a)
(3.13)
Proof. Write G the Green’s function of (1.1) and G∗ the
R Green’s function of (2.10).
−1
2
Observe that for f ∈ V and v = (div a∇) f , ku − vka = Ω2 (g(x) − f (x))G(x, y)(g(y) −
f (y)) dx dy. The monotonicity
of Green’s function as a quadratic form (see for
R
R instance
1
[12, Lemma 4.13]), implies Ω2 (g(x)−f (x))G(x, y)(g(y)−f (y)) dx dy ≤ λmin (a) Ω2 (g(x)−
f (x))G∗ (x, y)(g(y) − f (y)) dx dy (with a similar inequality on the l.h.s.) which concludes
the proof.
This extension, which is also directly related to the transfer property of the flux-norm
(introduced in [14] and generalized in [103], see also [115]), allows us to select accurate
finite dimensional bases for the approximation of the solution space of (1.1).
15

Construction 3.5. Let (τi )1≤i≤m be a partition of Ω such that each τi is Lipschitz,
2
convex and of diameter at most H. Let (φi )1≤i≤m be elements
R of L (Ω) such that for
each i, the support of φi is contained in the closure of τi and τi φi 6= 0.
Proposition 3.6. Let (φi )1≤i≤m be the elements of Construction 3.5 and let u be the
solution of (1.1). If V = span{φi : 1 ≤ i ≤ m} then
inf

v∈(div a∇)−1 V

ku − vka ≤ CHkgkL2 (Ω)

p
−1 
with C = π λmin (a)
1 + max1≤i≤m

1
φ2
|τi | τi i
R
1
( |τ | τ φi )2
i
i

R

1

)2



(3.14)

(writing |τi | the volume of τi ).

Proof. Using Proposition 3.4 it is sufficient to complete proof when a is the constant
Pm iden−1 V . Note that ∆v =
tity matrix. Let u0 be the solution
of
(2.10)
and
v
∈
∆
i=1 ci φi ,
R
R
R
Pm
0
2
0
therefore ku − vkH 1 (Ω) = − Ω (u − v)(g − i=1 ci φi ). Taking ci = τi g/ τi φi we obR
P0
0
2
tain that τi (g − m
j=1 cj φj ) = 0 and, writing |τi | the volume of τi , ku − vkH01 (Ω) =
R
Pm R
P
− i=1 τi (u0 − v − |τ1i | τi (u0 − v))(g − m
j=1 cj φj ) which by Poincaré’s inequality (see
Pm R
0
[91] for the optimal constant 1/π used here) lead to ku0 − vk2H 1 (Ω) ≤ H
i=1
π
τi |∇(u −
0
1
1 R
P
2 2 . Therefore, by using Cauchy-Schwartz inequality and simc
φ
)
v)|2 2 τi (g − m
j
j
j=1
Pm
plifying, ku0 − vkH01 (Ω) ≤ H
in τi
i=1 ci φi kL2 (Ω) . Now, since each φi has support
π kg −
R
R
1
2
2
R
φ
φ
P
Pm
|τi | τi i
2 R τi i ≤ kgk2
R
we have k i=1 ci φi k2L2 (Ω) = m
, which
i=1 ( τi g) (
L2 (Ω) max1≤i≤m ( 1
φ )2
φ )2
τi

i

concludes the proof.

τi

τi

i

The value of the constant C in Proposition 3.6 motivates us to modify Construction
3.5 as follows.
Construction 3.7. Let (φi )1≤i≤m be the elements constructed in 3.5 under the additional assumptions that (a) each φi is equal to one on τi and zero elsewhere (b) there
exists δ ∈ (0, 1) such that for each i ∈ {1, . . . , m}, τi contains a ball of diameter δH.
Let (φi )1≤i≤m be as in Construction 3.7. Note that the additional
assumption (a)
p
implies that the constant C in Proposition 3.6 is equal to 2/(π λmin (a)). Assumption
(b) will be used for localization purposes in subsections 3.5 and 3.6 (and is not required
for Theorem 3.8). The following theorem is a direct consequence of Proposition 3.6 and
Theorem 3.3.
Theorem 3.8. If u is the solution of (1.1) and (ψi )m
i=1 are the gamblets identified in
(3.5), (3.7) and (3.6) then
2
ku − vka ≤ p
HkgkL2 (Ω)
v∈span{ψ1 ,...,ψm }
π λmin (a)
inf

(3.15)

and the minimum in the l.h.s of (3.15) is achieved for v = u∗ defined in (3.9) and (3.10).

16

Remark 3.9. The assumption of convexity of the subdomains τi is only used to derive
sharper constants via Poincaré’s inequality for convex domains (without it, approximation error bounds remain valid after multiplication by π). Similarly, the transfer property
can be used to derive constructions that are distinct from 3.5 and 3.7.
Remark 3.10. Gamblets defined via the constrained energy minimization problems (3.7)
are analogous to the energy minimizing bases of [65, 114, 122, 121] and in particular
[108]. However they form a different set of basis functions when global constraints are
added: the (total) energy minimizing
bases of [65, 114, 122, 121,
P
P 108] are defined by
2
minimizing the total energy i kψi ka subject to the constraint i ψi (x) = 1 related to
the local preservation of constants. Numerical experiments
[122] suggest that total en√
ergy minimizing basis functions could lead to a O( H) convergence rate (with rough
coefficients). Note that (3.7) is also analogous to the constrained minimization problems
associated with Polyharmonic Splines [54, 30, 31, 32, 87], which can be recovered with a
Gaussian prior (on ξ) with covariance function δ(x − y) (corresponding to exciting (3.1)
with white noise). We suspect that the basis functions obtained in the orthogonal decomposition of [64] can also be recovered via the variational formulation (3.7) by identifying
the null space of the Clement quasi-interpolation operator with that of appropriately chosen measurement functions φi .

3.5

Exponential decay of gamblets

Theorems 3.2 and 3.3 show that the gamblets ψi have optimal recovery properties analogous to the discrete case of Theorem 2.2 and Corollary 2.6. However one may wonder
why one should compute these gamblets rather than the elements (div a∇)−1 φi since
they span the same linear space (by the representation formula (3.6)). The answer lies
in the fact that each gamblet ψi decays exponentially as a function of the distance from
the support of φi and its computation can therefore be localized to a subdomain of
diameter O(H ln H1 ) without impacting the order of accuracy (3.15). Consider the construction 3.7. Let ψi be defined as in Theorem 3.2 and let xi be an element of τi . Write
B(x, r) the ball of center x and radius r.
Theorem 3.11. Exponential decay of the basis elements ψi . It holds true that
Z
Z
r
T
1− lH
(∇ψi ) a∇ψi ≤ e
(∇ψi )T a∇ψi
(3.16)
Ω∩(B(xi ,r))c

with l = 1 + (e/π)

Ω

p
3
λmax (a)/λmin (a)(1 + 2 2 (2/δ)1+d/2 ) (where e is Euler’s number).

Proof. Let k, l ∈ N∗ and i ∈ {1, . . . , m}. Let S0 be the union of all the domains τj that
are contained in the closure of B(xi , klH) ∩ Ω, let S1 be the union of all the domains τj
that are contained in the closure of (B(xi , (k + 1)lH))c ∩ Ω and let S ∗ = S0c ∩ S1c ∩ Ω (be
the union of the remaining elements τj not contained in S0 or S1 ). Let η be the function
on Ω defined by η(x) = dist(x, S0 )/(dist(x, S0 )+dist(x, S1 )). Observe that (1) 0 ≤ η ≤ 1
1
(2) η is equal to zero on S0 (3) η is equal to one on S1 (4) k∇ηkL∞ (Ω) ≤ lH
. Observe that
17

R−

R

R

Ω ηψi div(a∇ψi ) = Ω ∇(ηψi )aψi
T
S1 (∇ψi ) a∇ψi ≤ I1 + I2 with

I1 = k∇ηkL∞ (Ω)

=

X Z
τj ⊂S ∗

R

T
Ω η(∇ψi ) a∇ψi

ψi2

Z

1

2

S∗

τj

+

R

T
Ω ψi (∇η) a∇ψi .

(∇ψi )T a∇ψi

1 p
2
λmax (a)

Therefore

(3.17)

R
and I2 = − Ω ηψi div(a∇ψi ). By (3.6), − div(a∇ψi ) is piecewise constant and equal to
R
Θ−1
i,j on τj . By the constraints of (3.7) τj ψi = 0 for i 6= j. Therefore (writing ηj the
volume average of η over τj ) we have
X Z
1 X
kψi kL2 (τj ) k div(a∇ψi )kL2 (τj ) . (3.18)
I2 ≤ −
(η − ηj )ψi div(a∇ψi ) ≤
l
∗
∗ τj
τj ⊂S

τj ⊂S1 ∪S

We will now need the following lemma
Lemma 3.12. If v ∈ span{ψ1 , . . . , ψm } then
1

k div(a∇v)kL2 (Ω) ≤ H −1 kvka (λmax (a)25+d /δ 2+d ) 2
(3.19)
P
P
−1
that − div(a∇v) = m
Proof. Let c ∈ Rm and v = m
i=1 ci Θi,j in τj
i=1 ci ψi . Observing
P
m
and using the decomposition k div(a∇v)k2L2 (Ω) = i=1 k div(a∇v)k2L2 (τj ) , we obtain that
k div(a∇v)k2L2 (Ω) =

m X
m
X
2
(
ci Θ−1
i,j ) |τj |

(3.20)

j=1 i=1

Furthermore,
v can be decomposed over τj as v = v1 +v2 , where v1 solves − div(a∇v1 ) =
Pm
−1
2 ) = 0 in τj with
i=1 ci Θi,j in τj with v1 = 0 on ∂τj , and v2 solves − div(a∇v
R
R
2
T
v2 = v on ∂τj . Using the notation |ξ|a = ξ aξ, observe that τj |∇v|2a = τj |∇v1 |2a +
R
R
R
Pm
−1
2
2
i=1 ci Θi,j τj v1 . Writing Gj the Green’s funcτj |∇v2 |a . Furthermore, τj |∇v1 |a =
tion
of the
operator −
R div(a∇·) with Dirichlet boundary condition on ∂τj , note that
R
Pm
−1
c
Θ
)
v
=
(
i
1
i=1
i,j
τj2 Gj (x, y) dx dy. Using the monotonicity of the Green’s function
τj
R
as a quadratic form (as in the proof of Proposition 3.4), we have τ 2 Gj (x, y) dx dy ≥
j
R
1
∗ (x, y) dx dy where G∗ is the Green’s function of the operator −∆ with
G
2
j
j
λmax (a) τj
R
Dirichlet boundary condition on ∂τj . Recall that 2 τj G∗j (x, y) dy is the mean exit time
(from τj ) of a Brownian motion started from x and the mean exit time of a Brownian
motion started from x to exit a ball of center x and radius r is Rr2 (see for instance
[12]). Since τj contains a ball of diameter δH, it follows that 2 τ 2 G∗j (x, y) dx dy ≥
j

(δH/4)2+d Vd (where Vd is the volume of the d-dimensional unit ball). Therefore (after
using |τj | ≤ Vd (H/2)d and simplification),
Z
m
X
2
2 2+d
|∇v1 |2a ≥ (
ci Θ−1
/(25+d λmax (a)),
(3.21)
i,j ) |τj |H δ
τj

i=1

which finishes the proof after taking the sum over j.
18

R
Now observe that since τj ψi = 0 for i 6= j, we obtain, using Poincaré’s inequality (with the optimal constant of [91]), that kψi kL2 (τj ) ≤ k∇ψi kL2 (τj ) H/π. Therefore,
combining (3.17), (3.18) and the result of Lemma 3.12, we obtain after simplification
Z
Z
3
1p
T
1+d/2
(∇ψi ) a∇ψi ≤
(∇ψi )T a∇ψi (3.22)
λmax (a)/λmin (a)(1 + 2 2 (2/δ)
)
πl
S1
S∗
p
3
1+d/2 ) and enlarging the integration domain
2
Taking l ≥ πe λmax (a)/λmin (a)(1
R + 2 (2/δ)
R
on the right hand side we obtain S1 (∇ψi )T a∇ψi ≤ e−1 S ∗ ∪S1 (∇ψi )T a∇ψi . We conclude
the proof via straightforward iteration on k.

3.6

Localization of the basis elements

Theorem 3.11 allows us to localize the construction of basis elements ψi as follows. For
r > 0 let Sr be the union of the subdomains τj intersecting B(xi , r) (recall that xi is an
element of τi ) and let ψiloc,r be the minimizer of the following quadratic problem
(
R
T
Minimize
Sr (∇ψ) a∇ψ
R
(3.23)
Subject to ψ ∈ H01 (Sr ) and Sr φj ψ = δi,j for j such that τj ⊂ Sr .
We will naturally identify ψiloc,r with its extension to H01 (Ω) by setting ψiloc,r = 0 outside
of Sr . From now on, to simplify the expression of constants, we will assume without loss
of generality that the domain is rescaled so that diam(Ω) ≤ 1.
Theorem 3.13. It holds true that
r

kψi − ψiloc,r ka ≤ Ce− 2lH ,
where l is defined in Theorem 3.11, C = (λmax (a)/
and Vd is the volume of the d-dimensional unit ball.

(3.24)

p
√
d
λmin (a))H − 2 −2 22d+9 /( Vd δ d+2 )

Proof. We will need the following lemma.
Lemma 3.14. It holds true that
d

kψi ka ≤ (Hδ)− 2 −1

p
3
1
λmax (a)2 2 d+2 (Vd )− 2

(3.25)

where Vd is the volume of the d-dimensional unit ball, and,
| ψi , ψj

ri,j

| ≤ e− 2lH H −2−d λmax (a)2
a

5d+11
2

/(Vd δ d+2 )

(3.26)

where l is the constant of Theorem 3.11 and ri,j is the distance between τi and τj .
Proof. Since τi contains a ball B(xi , δH/2) of center xi ∈ τi and diameter δH/2, there
exists a piece-wise differentiable function η, equal to 1 on B(xi , δH/4), equal to R0 on
4
(B(xi , δH/2))c and such that 0 ≤ η ≤ 1 with k∇ηkL∞ (Ω) ≤ Hδ
. Since ψ = η/( τi η)
satisfies the constrains of the minimization problem (3.7) we have kψi ka ≤ kψka , which
19

proves (3.25). Theorem 3.2 implies that ψi , ψj a = Θ−1
i,j . Observing that − div(a∇ψi )
−1
is piecewise constant and equal to Θi,j on τj and applying (3.21) (with v = ψi and using
R
R
2
2
τj |∇v1 |a ≤ τj |∇v|a ), we obtain that
Z
1
 1 −1
−1
5+d
2+d
2
(3.27)
(∇ψi )T a∇ψi 2 .
|Θi,j | ≤ λmax (a)2 /(δ |τj |) H
τj

which leads to (3.26) by the exponential decay obtained in Theorem 3.11 and (3.25).
Let us now prove Theorem 3.13. Let S0 be the union of the subdomains τj not
contained in Sr and let S1 be the union of the subdomains τj that are at distance at
least H from S0 (for S0 = ∅ the proof is trivial, so we may assume that S0 6= ∅, similarly
it is no restriction to assume that S1 6= ∅). Let η be the function on Ω defined by
η(x) = dist(x, S0 )/(dist(x, S0 ) + dist(x, S1 )). Observe that is a piecewise differentiable
function on Ω such that (1) η is equal to one on S1 and zero on S0 (2) k∇ηkL∞ (Ω) ≤ H1
and (3) 0 ≤ η ≤ 1. Since ψiloc,r satisfies the constraints of (3.7), we have from (3.8),
(3.28)
kψi − ψiloc,r k2a = kψiloc,r k2a − kψi k2a .
R
R
Let ψki,r be the minimizer of Sr (∇ψ)T a∇ψ subject to ψ ∈ H01 (Sr ) and Sr φj ψ = δk,j
R
Pm
i,r
i,r
i,r
:=
for τj ⊂ Sr . Write wj = Ω ηψi φj . Let ψw
j=1 wj ψj . Noting that ψw =
P
ψiloc,r + τj ⊂S ∗ wj ψji,r , where S ∗ is the union of τj ⊂ Sr not contained in S1 , and using
R
i,r T
i,r
property (3) of Theorem 3.2 (with Θi,−1
k,k0 = Sr (∇ψk ) a∇ψk0 ) it follows that
X
X i,−1
i,r 2
kψw
ka = kψiloc,r k2a + k
wj ψji,r k2a + 2
Θi,j wj .
(3.29)
τj ⊂S ∗

τj ⊂S ∗

i,r
Noting that ηψi ∈ H01 (Sr ), Theorem 3.2 implies that kψw
ka ≤ kηψi ka , which, combined
P
with (3.29) and (3.28) leads to kψi − ψiloc,r k2a ≤ kηψi k2a − kψi k2a − 2 τj ⊂S ∗ Θi,−1
i,j wj and
R
2
2
T
(using kηψi ka − kψi ka ≤ S∗ ∇(ηψi ) a∇(ηψi ))
Z
X i,−1
kψi − ψiloc,r k2a ≤
∇(ηψi )T a∇(ηψi ) + 2|
Θi,j wj | .
(3.30)
S∗

τj ⊂S ∗

R
R
(a) R
2
Now observe that 12 S∗ ∇(ηψi )T a∇(ηψi ) ≤ Ω∩(B(xi ,r−2H))c (∇ψi )T a∇ψi + λmax
2
S ∗ |ψi | .
H
R
R
P
Applying Poincaré’s inequality we obtain S ∗ |ψi |2 ≤ π12 H 2 τj ⊂S ∗ τj |∇ψi |2 (since
R
R
R
H2
∗ ), and
2 ≤
|ψ
|
(∇ψi )T a∇ψi . Comψ
=
0
for
τ
⊂
S
∗
i
j
i
2
τj
S
π λmin (a) Ω∩(B(xi ,r−2H))c
bining these equations with the exponential decay of Theorem 3.11 we deduce
Z


r−2H
∇(ηψi )T a∇(ηψi ) ≤ 2 1 + λmax (a)/ π 2 λmin (a) e1− lH kψi k2a .
(3.31)
S∗

Similarly, using Cauchy-Schwartz and Poincaré inequalities we have for τj ⊂ S ∗ ,
1
1 R
1 p
|wj | ≤ |τj | 2 kψi kL2 (τj ) ≤ |τj | 2 ( τj (∇ψi )T a(∇ψi )) 2 / λmin (a) and
20

1
P
1 R
i,−1 2
T
2
2
Θi,−1
τj ⊂S ∗ (Θi,j ) |τj ||
i,j wj | ≤ |
S ∗ (∇ψi ) a(∇ψi )/λmin (a) . Using (3.27)
1
1
R
P
1
i,r T
i,r  2
5+d /δ 2+d 2 H −1
2
2
,
we obtain that | τj ⊂S ∗ (Θi,−1
i,j ) |τj || ≤ λmax (a)2
S ∗ (∇ψi ) a∇ψi

|

P

τj ⊂S ∗

which by the exponential decay of Theorem (3.11) (and kψi ka ≤ kψii,r ka ) leads to
|

X

Θi,−1
i,j wj | ≤

τj ⊂S ∗

λmax (a)25+d  21 −1 i,r 2 1− r−2H
lH
H kψi ka e
.
λmin (a)δ 2+d

(3.32)

Using (3.25) to bound kψii,r ka and combining (3.32) with (3.31) and (3.30) concludes
the proof.
The following theorem shows that gamblets preserve the O(H) rate of convergence (in
energy norm) after localization to sub-domains of size O(H ln(1/H)). They can therefore
be used as localized basis functions in numerical homogenization [5, 85, 64, 87]. Section
4 will show that they can also be computed hierarchically at near linear complexity.
Theorem 3.15. Let u be the solution of (1.1) and (ψ1loc,r )1≤i≤m the localized gamblets
identified in (3.23), then for r ≥ H(C1 ln H1 + C2 ) we have
inf

loc,r
v∈span{ψ1loc,r ,...,ψm
}

1
HkgkL2 (Ω) .
ku − vka ≤ p
λmin (a)

The constants are C1 = (d + 4)l and C2 = 2l ln



3

λmax (a) 2 2 d+11
λmin (a) δ d+2

(3.33)



where l is the constant of
P
loc,r R
Theorem 3.13. Furthermore, the inequality (3.33) is achieved for v = m
i=1 ψi
Ω uφi .
R
Pm
P
loc,r
with ci = Ω uφi . Theorem 3.8 implies
Proof. Let v1 := m
i ψi and v2 =
i=1 ci ψi
i=1 cp
that ku − v1 ka ≤ 2/(π λmin (a))HkgkL2 (Ω) . Observe that ku − v2 ka ≤ ku − v1 ka +
P
kv1 − v2 ka and kv1 − v2 ka ≤ maxi kψi − ψiloc,r ka m
Poincaré’s inequality
i=1 |ci |. Using P
R
m
kukL2 (Ω) ≤ diam(Ω)k∇ukL2 (Ω) (with diam(Ω) ≤ 1) we obtain
i=1 |ci | ≤ Ω |u| ≤
1

kgkL2 (Ω) 2−d/2 Vd2 /λmin (a). We conclude using Theorem 3.13 to bound maxi kψi −ψiloc,r ka .

4

Multiresolution operator decomposition

Building on the analysis of Section 3, we will now gamble on the approximation of the
solution of (1.1) based on measurements performed at different levels of resolution. The
resulting hierarchical (and nested) games will then be used to derive a multiresolution
decomposition of (1.1) (orthogonal across subscales) and a near-linear complexity multiresolution algorithm with a priori error bounds.

4.1

Hierarchy of nested measurement functions

In order to define the hierarchy of games we will first define a hierarchy of nested measurement functions.
21

Definition 4.1. We say that I is an index tree of depth q if it is a finite set of qtuples of the form i = (i1 , . . . , iq ) with 1 ≤ i1 ≤ m0 and 1 ≤ ij ≤ m(i1 ,...,ij−1 ) for
j ≥ 2, where m0 and m(i1 ,...,ij−1 ) are strictly positive integers. For 1 ≤ k ≤ q and
i = (i1 , . . . , iq ) ∈ I, we write i(k) := (i1 , . . . , ik ) and I (k) := {i(k) : i ∈ I}. For
0
k ≤ k 0 ≤ q and j = (j1 , . . . , jk0 ) ∈ I (k ) we write j (k) := (j1 , . . . , jk ). For i ∈ I (k) and
0
0
k ≤ k 0 ≤ q we write i(k,k ) the set of elements j ∈ I (k ) such that j (k) = i.
Construction 4.2. Let I be an index tree of depth q. Let δ ∈ (0, 1) and 0 < Hq < · · · <
(k)
H1 < 1. Let (τi , k ∈ {1, . . . , q}, i ∈ I (k) ) be a collection of subsets of Ω such that (1)
(k)
(k)
for 1 ≤ k ≤ q, (τi , i ∈ I (k) ) is a partition of Ω such that each τi is a Lipschitz, convex
subset of Ω of diameter at most Hk and contains a ball of diameter δHk (2) the sequence
(k)
(k+1)
of partitions is nested, i.e. for k ∈ {1, . . . , q − 1} and i ∈ I (k) , τi := ∪j∈i(k,k+1) τj
.
(k)

As in Remark 3.9, the assumption of convexity of the subdomains τi is not necessary
to the results presented here and is only used to derive sharper/simpler constants. Let
(k)
(k)
(k)
(k)
(k)
φi be the indicator function of the set τi (i.e. φi = 1 if x ∈ τi and φi = 0
(k)
if x 6∈ τi ). Note that the nesting of the domain decomposition implies that of the
measurement functions, i.e. for k ∈ {1, . . . , q − 1} and i ∈ I (k) ,
X
(k)
(k,k+1) (k+1)
φi =
πi,j
φj
(4.1)
j∈I (k+1)
(k,k+1)

where π (k,k+1) is the I (k) × I (k+1) matrix defined by πi,j
(k,k+1)

πi,j

= 1 if j ∈ i(k,k+1) and
(k)

= 0 if j 6∈ i(k,k+1) . We will assume without loss of generality that kφi k2L2 (Ω) =

(k)

(k)

|τi | is constant in i (for the general case, rescale/renormalize each φi and the entries of
π (k,k+1) by the corresponding multiplicative factors, we will keep track of the dependence
(k)
(k)
of some of the constants on maxi,j |τi |/|τj |).

4.2

Hierarchy of nested gamblets and multiresolution approximations

Let us now consider the problem of recovering the solution of (1.1) based on the nested
R
(k)
measurements ( Ω uφi )i∈I (k) for k ∈ {1, . . . , q}. As in Section 3 we are lead to investigate the mixed strategy (for Player II) expressed by replacing the source term g with
a centered Gaussian field with covariance function L = − div(a∇). Under that mixed
strategy, Player II’s bet on the value of the solution of (1.1), given the measurements
R
(k)
( Ω u(y)φi (y) dy)i∈I (k) , is (see Subsection 3.3)
Z
X (k)
(k)
u(k) (x) :=
ψi (x) u(y)φi (y) dy,
(4.2)
Ω

i∈I (k)

(k)

where (see Theorem 3.2), for k ∈ {1, . . . , q} and i ∈ I (k) , ψi is the minimizer of
(
Minimize
kψka
(4.3)
R (k)
Subject to ψ ∈ H01 (Ω) and Ω φj ψ = δi,j for j ∈ I (k) .
22

Define V(q+1) := H01 (Ω) and, for k ∈ {1, . . . , q},
(k)

V(k) := span{ψi

| i ∈ I (k) }.

(k)

(4.4)

(k)

By Theorem 3.1 span{ψi | i ∈ I (k) } = span{L−1 φi | i ∈ I (k) }, and the nesting
(4.1) of the measurement functions implies the nesting of the spaces V(k) . The following
theorem is (which is a direct application of theorems 3.3 and 3.8) shows that u(k) is the
best (energy norm) approximation of the solution of (1.1) in V(k) .
Theorem 4.3. It holds true that (1) for k ∈ {1, . . . , q}, V(k) ⊂ V(k+1) and V(k) =
(k)
span{L−1 φi | i ∈ I (k) } and (2) If u is the solution of (1.1) and u(k) defined in (4.2)
then
2
ku − u(k) ka = inf ku − vka ≤ p
Hk kgkL2 (Ω)
(4.5)
v∈V(k)
π λmin (a)

4.3

Nested games and martingale/multiresolution decomposition

As in Section 3 we consider the mixed strategy (for Player II) expressed by replacing
the source term g with a centered Gaussian field with covariance function L. Under
this mixed strategy, Player II’s bet (4.2) on the value of the solution of (1.1), given
R
(k)
the measurements ( Ω u(y)φi (y) dy)i∈I (k) , can also be obtained by conditioning the
solution v of the SPDE (3.1) (see (3.10)), i.e.
Z
Z
i
h
(k)
(k)
(k)
(4.6)
u (x) = E v(x)
v(y)φi (y) dy =
u(y)φi (y) dy, i ∈ I (k)
Ω

Ω

(k)
ψi
R

Furthermore, each gamblet
(1.1) given the measurements
(k)

ψi

represents Player II’s bet on the value of the solution of
(k)
Ω u(y)φj (y) dy = δi,j , i.e.
i
h Z
(k)
(4.7)
=E v
v(y)φj (y) dy = δi,j , j ∈ I (k)
Ω

Now consider the nesting of non-cooperative games where Player I chooses g in (1.1)
R
(k)
and Player II is shown the measurements ( Ω uφi )i∈I (k) , step by step, in a hierarchical
manner, from coarse (k = 1) to fine (k = q) and must, at each step k of the game, gamble
on the value of solution u. The following theorem and (4.6) show that the resulting
sequence of approximations u(k) form the realization of a martingale with independent
increments.
R
(k)
Theorem 4.4. Let Fk be the σ-algebra generated by the random variables ( Ω v(x)φi )i∈I (k)
and
Z
X (k)


(k)
(k)
v (x) := E v(x) Fk =
ψi (x) v(y)φi (y) dy
(4.8)
i∈I (k)

Ω

It holds true that (1) F1 , . . . , Fq forms a filtration, i.e. Fk ⊂ Fk+1 (2) For x ∈ Ω, v (k)
 (x)
is a martingale with respect to the filtration (Fk )k≥1 , i.e. v (k) (x) = E v (k+1) (x) Fk (3)
v (1) and the increments (v (k+1) − v (k) )k≥1 are independent Gaussian fields.
23

Proof. The nesting (4.1) of the measurement functions implies Fk ⊂ Fk+1 and (F
 k )k≥1

(k)
(k)
is therefore filtration. The fact that v
is a martingale follows from v = E v Fk .
Since v (1) and the increments (v (k+1) − v (k) )k≥1 are Gaussian fields belonging to the
same Gaussian space their independence is equivalent to zero covariance,
h which follows
 (1) (k+1)

(k)
from the martingale property, i.e. for k ≥ 1 E v (v
− v ) = E E v (1) (v (k+1) −
i
h
i




v (k) ) Fk = E v (1) E (v (k+1) −v (k) ) Fk = 0 and for k > j ≥ 1, E (v (j+1) −v (j) )(v (k+1) −
h


i
v (k) ) = E (v (j+1) − v (j) )E (v (k+1) − v (k) ) Fk = 0.
Remark 4.5. Theorem 4.4 enables the application of classical results concerning martingales to the numerical analysis of v (k) (and u(k) ). In particular (1) Martingale (concentration) inequalities can be used to control the fluctuations of v (k) (2) Optimal stopping
times can be used to derive optimal strategies for stopping numerical simulations based
on loss functions mixing computation costs with the cost of imperfect decisions (3) Tak(k)
ing q = ∞ in the construction of the basis elements ψi (with a sequence Hk decreasing
0) and
using the martingale convergence theorem imply that, for all ϕ ∈ C0∞ (Ω),
Rtowards
R
(k) ϕ →
1
Ωv
Ω vϕ as k → ∞ (a.s. and in L ).
The independence of the increments v (k+1) −v (k) is related to the following orthogonal
multiresolution decomposition of the operator (1.1). For V(k) defined as in (4.4) and for
k ∈ {2, . . . , q + 1} let W(k) be the orthogonal complement of V(k−1) within V(k) with
respect to the scalar product ·, · a . Write ⊕a the orthogonal direct sum with respect to
the scalar product ·, · a . Note that by Theorem 4.3, u(k) defined by (4.2) is the finite
element solution of (1.1) in V(k) (in particular we will write u(q+1) = u).
Theorem 4.6. It holds true that (1) For k ∈ {2, . . . , q + 1},
V(k) = V(1) ⊕a W(2) ⊕a · · · ⊕a W(k) ,

(4.9)

(2) for k ∈ {1, . . . , q}, u(k+1) − u(k) belongs to W(k+1) and
u = u(1) + (u(2) − u(1) ) + · · · + (u(q) − u(q−1) ) + (u − u(q) )

(4.10)

is the orthogonal decomposition of u in H01 (Ω) = V(1) ⊕a W(2) ⊕a · · · ⊕a W(q) ⊕a W(q+1) ,
and (3) u(k+1) − u(k) is the finite element solution of (1.1) in W(k+1) .
Proof. Observe that since the V(k) are nested (Theorem 4.3) u(k+1) − u(k) belongs to
V(k+1) . Furthermore (by Property (1) of Theorem 4.3 and integration by parts), for
R
(k)
(k)
i ∈ I (k) , u(k+1) − u(k) , ψi a belongs to span{ Ω (u(k+1) − u(k) )φi | i ∈ I (k) }. Finally,
(4.2), the constraints of (4.3) and the nesting property (4.1) imply that for i ∈ I (k) ,
R
R
P
(k,k+1) R
(k+1)
(k)
(k+1) − u(k) )φ(k) =
− Ω uφi = 0 which implies that
j∈I (k+1) πi,j
i
Ω (u
Ω uφj
u(k+1) − u(k) belongs to W(k+1) .

24

4.4

Interpolation and restriction matrices/operators

Since the spaces V(k) are nested there exists a I (k) × I (k+1) matrix R(k,k+1) such that
for 1 ≤ k ≤ q − 1 and i ∈ I (k)
X
(k)
(k,k+1) (k+1)
ψi =
Ri,j
ψj
(4.11)
j∈I (k+1)

We will refer to R(k,k+1) as the restriction matrix and to its transpose R(k+1,k) :=
(R(k,k+1) )T as the interpolation/prolongation matrix. The following theorem shows that
R
(k,k+1)
(k+1)
(see Figure 1) Ri,j
is Player II’s best bet on the value of Ω uφj
given the inforR
(k)
(k)
mation that Ω uφs = δi,s , s ∈ I ).
Theorem 4.7. It holds true that for i ∈ I (k) and j ∈ I (k+1) ,
R

R (k) (k+1)
R
(k,k+1)
(k+1)
(k)
= E Ω v(y)φj
(y) dy Ω v(y)φl (y) dy = δi,l , l ∈ I (k) .
Ri,j
= Ω ψi φj
(k+1)

Proof. The first equality is obtained by integrating (4.11) against φj

and using the

(k+1)
ψj

constraints satisfied by
in (4.3). For the second equality, observe that since Fk
is a filtration we can replace v in the representation formula (4.7) by v (k) (as defined by
the r.h.s. of (4.8)) and obtain
R

R
P
(k)
(k+1)
(k+1)
(k)
(x)E Ω v(y)φj
(y) dy Ω v(y)φl (y) dy = δi,l , l ∈ I (k)
ψi (x) =
j∈I (k+1) ψj
which corresponds to (4.11).

4.5

Nested computation of the interpolation and stiffness matrices

R
(k)
Let v be the solution of (3.1). Observe that ( Ω v(x)φi )i∈I (k) is a Gaussian vector with
(symmetric, positive definite) covariance matrix Θ(k) defined by for i, j ∈ I (k) ,
Z
(k)
(k)
(k)
Θi,j :=
φi (x)G(x, y)φj (y) dx dy .
(4.12)
Ω2

As in (3.3), Θ(k) is invertible and we write Θ(k),−1 its inverse. Observe that, as in
(k)
Theorem 3.2, ψi admits the following representation formula
X (k),−1 Z
(k)
(k)
ψi (x) =
Θi,j
G(x, y)φj (y) dy
(4.13)
Ω

j∈I (k)

Observe that, as in Theorem 3.2, Θ(k),−1 = A(k) where A(k) is the (symmetric,
(k)
positive definite) stiffness matrix of the elements ψi , i.e., for i, j ∈ I (k) ,
(k)

(k)

(k)

Ai,j := ψi , ψj

a

(4.14)

Write π (k+1,k) the transpose of the matrix π (k,k+1) (defined below (4.1)) and I (k)
the I (k) × I (k) identity matrix. The following theorem enables the hierarchical/nested
computation of A(k) from A(k+1) .
25

Theorem 4.8. For b ∈ RI

(k)

, R(k+1,k) b is the (unique) minimizer c ∈ RI
(
Minimize
cT A(k+1) c
Subject to π (k,k+1) c = b

(k+1)

of
(4.15)

Furthermore R(k,k+1) π (k+1,k) = π (k,k+1) R(k+1,k) = I (k) , R(k,k+1) = A(k) π (k,k+1) Θ(k+1) ,
Θ(k) = π (k,k+1) Θ(k+1) π (k+1,k) and
A(k) = R(k,k+1) A(k+1) R(k+1,k) .
R (k) (k)
Proof. Using the decompositions (4.11) and (4.1) in Ω φj ψi = δi,j leads to

(4.16)

(k)

R(k,k+1) π (k+1,k) = I (k) . Using (4.13) and (4.1) to expand ψi Theorem 4.7 leads to
(k)
(k)
R(k,k+1) = A(k) π (k,k+1) Θ(k+1) . Using (4.1) to expand φi and φj in (4.12) leads to
(k)

(k)

Θ(k) = π (k,k+1) Θ(k+1) π (k+1,k) . Using (4.11) to expand ψi and ψj in (4.14) leads to
P
(k)
(k)
(4.16). Let b ∈ RI . Theorem 3.2 implies that i∈I (k) bi ψi is the unique minimizer
R
(k)
of kvk2a subject to v ∈ H01 (Ω) and Ω φj v = bj for j ∈ I (k) . Since V(k) ⊂ V(k+1)
and since the minimizer is in V(k) , the minimization over v ∈ H01 (Ω) can be reduced to
P
(k+1)
v ∈ V(k+1) of the form v = i∈I (k+1) ci ψi
, which after using (4.1) to expand the
R (k)
constraint Ω φj v = bj , corresponds to (4.15).

4.6

Multiresolution gamblets

The interpolation and restriction operators are sufficient to derive a multigrid method
for solving (1.1). To design a multiresolution algorithm we need to continue the analysis
and identify basis functions for the subspaces W(k) . For k = 2, . . . , q let J (k) be the
finite set of k-tuples of the form i = (i1 , . . . , ik ) with 1 ≤ i1 ≤ m0 , 1 ≤ ij ≤ m(i1 ,...,ij−1 )
for 2 ≤ j ≤ k − 1 and 1 ≤ ik ≤ m(i1 ,...,ik−1 ) − 1. Where the integers m· are the same
as those defining the index tree I. For a matrix M write Im(M ) and Ker(M ) its image
and kernel.
Lemma 4.9. For k = 2, . . . , q let W (k) be a J (k) × I (k) matrix such that Im(W (k),T ) =
(k)
Ker(π (k−1,k) ). It holds true that the elements (χi )i∈J (k) ∈ V(k) defined as
X
(k)
(k) (k)
χi :=
Wi,j ψj
(4.17)
j∈I (k)

form a basis of W(k) .
(k−1)

Proof. Since LV(k−1) = span{φi
| i ∈ I (k−1) }, w ∈ V(k) belongs to W(k) if and only if
R (k−1)
(k)
w = 0 for all j ∈ I (k−1) , which, taking w = χi and using (4.1), translates into
Ω φj
(π (k−1,k) W (k),T )j,i = 0. Writing |J (k) | the number of elements of J (k) (which is equal to
the dimension of W(k) ), observe that |J (k) | = |I (k) | − |I (k−1) |. Therefore Im(W (k),T ) =
(k)
Ker(π (k−1,k) ) also implies that the |J (k) | elements χi are linearly independent and,
therefore, form a basis of W(k) .
26

(k−1)

(k)

Remark 4.10. Observe that since 0 = ψi
, χj a = (R(k−1,k) A(k) W (k),T )i,j , it also
holds true that Im(W (k),T ) = Ker(R(k−1,k) A(k) ) and Im(A(k) W (k),T ) = Ker(R(k−1,k) ) .
From now on we choose, for each k ∈ {2, . . . , q}, a J (k) × I (k) matrix W (k) as in
Lemma 4.9. This choice is not unique and to enable fast multiplication by W (k) (or its
(k)
transpose) we require that for (j, i) ∈ J (k) × I (k) , Wj,i = 0 if j (k−1) 6= i(k−1) . Therefore,
the construction of W (k) requires, for each s ∈ I (k−1) , to specify a number ms − 1 of
(k)
(k)
ms -dimensional vectors W(s,1),(s,·) , . . . , W(s,ms −1),(s,·) that are linearly independent and
orthogonal to the ms -dimensional vector (1, 1, . . . , 1, 1). We propose two simple constructions.
(k)

Construction 4.11. For k ∈ {2, . . . , q}, choose W (k) such (1) Wj,i = 0 for (j, i) ∈
J (k) × I (k) with j (k−1) 6= i(k−1) and (2) for s ∈ I (k−1) , t ∈ {1, . . . , ms − 1} and t0 ∈
(k)
{1, . . . , ms }, W(s,t),(s,t0 ) = δt,t0 − δt+1,t0 .
For k ∈ {2, . . . , q} and i = (i1 , . . . , ik−1 , ik ) ∈ J (k) define i+ := (i1 , . . . , ik−1 , ik + 1)
and observe that under construction 4.11,
(k)

χi

(k)

= ψi

(k)

− ψi+

(4.18)

whose game-theoretic interpretation is provided in Figure 1.

(k)

(k)

Figure 1: If (τs , s ∈ I (k) ) is a nested rectangular partitionR of Ω then (a) ψi is Player
II’s best bet on the value of the solution u of (1.1) given τ (k) u = δi,j for j ∈ I (k) (b)
j
R
(k)
(k,k+1)
χi is Player II’s best bet on u given τ (k) u = δi,j − δi+ ,j for j ∈ I (k) (c) Ri,j
is
j
R
R
Player II’s best bet on τ (k+1) u given τ (k) u = δi,j for j ∈ I (k) .
j

j

For the second construction we need the following lemma whose proof is trivial.
Lemma 4.12. Let U (n) be the sequence of n × n matrices defined (1) for n = 2 by
(2)
(2)
(n+1)
(n)
U1,· = (1, −1) and U2,· = (1, 1) and (2) iteratively for n ≥ 2 by Ui,j
= Ui,j for
27

(n+1)

(n+1)

1 ≤ i, j ≤ n, Un+1,j = 1 for 1 ≤ j ≤ n + 1, Ui,n+1 = 0 for 1 ≤ i ≤ n − 1 and
(n+1)

(n)

Un,n+1 = −n. Then for n ≥ 2, the rows of U (n) are orthogonal, Un,j = 1 for 1 ≤ j ≤ n
and we write Ū (n) the corresponding orthonormal matrix obtained by renormalizing the
rows of U (n) .
(k)

Construction 4.13. For k ∈ {2, . . . , q}, choose W (k) such (1) Wj,i = 0 for (j, i) ∈
J (k) × I (k) with j (k−1) 6= i(k−1) and (2) for s ∈ I (k−1) and t ∈ {1, . . . , ms − 1} and
(k)
(m )
t0 ∈ {1, . . . , ms } W(s,t),(s,t0 ) = Ūt,t0 s (where Ū (ms ) is defined in Lemma 4.12).
Observe that under Construction 4.13 (1) the complexity of constructing W (k) is
× m2s and (2) W (k) W (k),T = J (k) where J (k) is the J (k) × J (k) identity matrix.

|I (k−1) |

4.7

Multiresolution operator inversion
(1)

(k)

We will now use the basis functions ψi and χi to perform the multiresolution inversion
(k)
(k)
(k)
of (1.1). Let B (k) be the J (k) × J (k) (stiffness)matrix Bi,j = χi , χj a and observe
that
B (k) = W (k) A(k) W (k),T
(4.19)
Observe that B (k) is positive, symmetric, definite and write B (k),−1 its inverse. Let
π̄ (k,k+1) be the I (k) × I (k+1) matrix defined by
(k,k+1)

π̄i,j

(k,k+1)

= πi,j

/(π (k,k+1) π (k+1,k) )i,i

(4.20)

Using the notations of Definition 4.1 note that (π (k,k+1) π (k+1,k) )i,i = mi . Let D(k,k−1)
be the J (k) × I (k−1) matrix defined as
D(k,k−1) := −B (k),−1 W (k) A(k) π̄ (k,k−1)

(4.21)

and write D(k−1,k) := D(k,k−1),T its transpose.
Theorem 4.14. It holds true that for k ∈ {1, . . . , q − 1} and i ∈ I (k) ,
X
X
(k)
(k,k+1) (k+1)
(k,k+1) (k+1)
ψi =
π̄i,l
ψl
+
Di,j
χj
.
l∈I (k+1)

(4.22)

j∈J (k+1)

In particular,
Proof. For s ∈ I (k)
I (k) }. Let x ∈ RI

(k)

R(k,k+1) = π̄ (k,k+1) + D(k,k+1) W (k+1)
(4.23)
P
(k)
(k,k+1) (k+1)
(k)
write ψ̄s := l∈I (k+1) π̄s,l
ψl
and V̄(k) := span{ψ̄s | s ∈
, y ∈ RJ

(k+1)

and
X
ψ=
xs ψ̄s(k) +
s∈I (k)

X
j∈J (k+1)

28

(k+1)

y j χj

.

(4.24)

(k)

If ψ = 0 then integrating ψ against φi

for i ∈ I (k) (and observing that

R

(k) (k)
Ω φi ψ̄s

=

(k)
(k+1)
δi,s ) implies x = 0 and y = 0. Therefore the elements ψ̄s , χj
form a basis for
(k)
(k+1)
(k+1)
(k)
(k+1)
V̄ + W
. Observing that dim(V
) = dim(V̄ ) + dim(W
) we deduce that
(k)

V(k+1) = V̄(k) + W(k+1) . Therefore, since V(k) ⊂ V(k+1) , ψi can be decomposed as in
R (k) (k)
(4.24). The constraints Ω φs ψi = δi,s lead to xs = δi,s . The orthogonality between
(k+1)

ψ and W(k+1) leads to the equations ψ, χj
= 0 for j ∈ J (k+1) , i.e.
a
P
P
(k,k+1)
(k+1)
(k+1)
(k+1)
(k+1)
ψl
, χj
+ j 0 ∈J (k+1) yj 0 χj 0
, χj
= 0, which transl∈I (k+1) π̄i,l
a
a
(k+1,k)

lates into W (k+1) A(k+1) π̄·,i
+ B (k+1) y, that is (4.22). Plugging (4.17) in (4.22) and
comparing with (4.11) leads to (4.23).
Let g be the r.h.s of (1.1). For k ∈ {1, . . . , q} let g (k) be the |I (k) |-dimensional vector
R (k)
(k)
defined by gi = Ω ψi g for i ∈ I (k) . Observe that g (k) can be computed iteratively
using
g (k) = R(k,k+1) g (k+1) .
(4.25)
For k ∈ {2, . . . , q}, let w(k) be |J (k) |-dimensional vector defined as the solution of
B (k) w(k) = W (k) g (k)

(4.26)

Furthermore let U (1) be the |I (1) |-dimensional vector defined as the solution of
A(1) U (1) = g (1)

(4.27)

According to following theorem, which is a direct consequence of Theorem 4.6, the
solution of (1.1) can be computed at any scale by solving the decoupled linear systems
(4.26) and (4.27).
Theorem 4.15. For k ∈ {2, . . . , q}, let u(k) be the finite element solution of (1.1) in
P
(k) (k)
V(k) . It holds true that u(k) − u(k−1) = i∈J (k) wi χi and, in particular,
(k)

u

=

X

(1) (1)
Ui ψi

+

(k0 ) (k0 )
χi

wi

(4.28)

k0 =2 i∈J (k0 )

i∈I (1)

4.8

k
X
X

Uniformly bounded condition numbers across subscales/subbands
(k)

Taking q = ∞ in Theorem 4.6, the construction of the basis elements ψi
multiresolution orthogonal decomposition,
∞

H01 (Ω) = V(1) ⊕a W(i) .

leads to the

(4.29)

i=2

(k)

(k)

In that sense the basis elements ψi and χi could be seen as a generalization of wavelets
to the orthogonal decomposition of H01 (Ω) (rather than L2 (Ω)) adapted to the solution
space of the PDE (1.1). We will now show that this orthogonal decomposition induces
29

a subscale decomposition of the operator − div(a∇) into layered subbands of increasing
frequencies. Moreover the condition number of the operator − div(a∇) restricted to
each subspace W(k) will be shown to be uniformly bounded if Hk−1 /Hk is uniformly
bounded (e.g. if Hk is a geometric sequence). Write H0 := 1 and let δ be defined as in
Construction 4.2.
Theorem 4.16. If k ∈ {1, . . . , q} and v ∈ V(k) then
δ 1+d/2
kvka
p
Hk ≤
.
5/2+d/2
k div(a∇v)kL2 (Ω)
λmax (a)2

(4.30)

If (k = 1 and v ∈ V(1) ) or (k ∈ {2, . . . , q} and v ∈ W(k) ) then
kvka
1
Hk−1 .
≤p
k div(a∇v)kL2 (Ω)
λmin (a)

(4.31)

Proof. (4.30) is a direct consequence of Lemma 3.12. For k = 1 (4.31) is a simple
consequence of Poincaré’s inequality. Let k ∈ {2, . . . , q}. V(k) = V(k−1) ⊕a W(k) and
Theorem 4.3 imply
kvka
kv−v 0 ka
supv∈W(k) k div(a∇v)k
≤
sup
≤ √ 2
Hk−1 .
(k) inf v 0 ∈V(k−1)
v∈V
k
div(a∇v)k
2
2
L (Ω)

L (Ω)

π

λmin (a)

Write |c| the Euclidean norm of c and for k ∈ {1, . . . , q} let
k
γk := inf
(k)
¯
c∈RI

P

(k) 2
i∈I (k) ci φi kL2 (Ω)
|c|2

k
and γ̄k := sup
c∈RI

(k) 2
i∈I (k) ci φi kL2 (Ω)
|c|2

P

(k)

(k)

(4.32)
(k)

Write |τ | the volume of a set τ and note that γ̄k ≤ maxi∈I (k) |τi | and γk ≥ mini∈I (k) |τi |,
¯
therefore γ̄k /γk ≤ δ −d .
p
p
¯
For a given matrix M , write Cond(M ) := λmax (M T M )/ λmin (M T M ) its condition number.
Theorem 4.17. It holds true that
1 λmax (a)25+d
,
H12 λmin (a)δ 2+2d

(4.33)

Hk−1 2 λmax (a) 2 211+2d
Cond(W (k) W (k),T ) .
Hk
λmin (a) δ 4+7d π 2

(4.34)

Cond(A(1) ) ≤
and for k ∈ {2, . . . , q},
Cond(B (k) ) ≤

Furthermore, Cond(W (k) W (k),T ) = 1 under Construction 4.13 and Cond(W (k) W (k),T ) ≤
2d
2 Hk−1 /(δHk )
under Construction 4.11.

30

P
(k)
. Write v =
i∈I (k) ci ψi . Observing that
P
(k)
kvk2a = cT A(k) c and k div(a∇v)k2L2 (Ω) = k i∈I (k) (A(k) c)i φi k2L2 (Ω) ≥ γk |A(k) c|2 , (4.30)
¯
implies that γk Hk2 δ 2+d /(λmax (a)25+d ) ≤ cT A(k) c/|A(k) c|2 , which, after taking the mini¯ to (for k ≥ 1)
mum in c leads

Proof. Let k ∈ {1, . . . , q} and c ∈ RI

(k)

λmax (A(k) ) ≤ λmax (a)25+d /(Hk2 δ 2+d γk ),
¯
and for k ≥ 2 (using (4.19))

(4.35)

λmax (B (k) ) ≤ λmax (W (k) W (k),T )λmax (a)25+d /(Hk2 δ 2+d γk ) .
¯

(4.36)

Similarly for k = 1 (4.31) leads to λmin (A(1) ) ≥ λmin (a)/γ̄1 . Now let us consider
P
(k)
(k) (k)
k ∈ {2, . . . , q} and c ∈ RJ . Write w = i∈J (k) , j∈I (k) ci Wi,j ψj . (4.17) and (4.19)
imply that kwk2a = cT B (k) c and (using (4.32))
P
(k)
k div(a∇w)kL2 (Ω) = k i∈J (k) , j∈I (k) (A(k) W (k),T c)j φj kL2 (Ω) ≤ γ̄k |A(k) W (k),T c|2 . Observing that w ∈ W(k) , (4.31) implies that
c = B (k),−1 y for y ∈ RJ

(k)

we deduce that

cT B (k) c
2 . Taking
≤ γ̄k λmin1 (a) Hk−1
cT W (k) (A(k) )2 W (k),T c
y T B (k),−1 y
2 . Writing
≤ γ̄k λmin1 (a) Hk−1
(k)
|A W (k),T B (k),−1 y|2

N (k) = −A(k) W (k),T B (k),−1 , we have obtained that

2
γ̄k λmax (N (k),T N (k) ) ≤ λmin (B (k) ) .
λmin (a)/ Hk−1

(4.37)

For k ∈ {2, . . . , q} let P (k) := π (k,k−1) R(k−1,k) . Using R(k−1,k) = A(k−1) π (k−1,k) Θ(k) and
π (k−1,k) Θ(k) π (k,k−1) = Θ(k−1) (Theorem 4.8) we obtain that (P (k) )2 = P (k) , i.e. P (k) is
a projection. Write kP (k) kKer(π(k−1,k) ) := supx∈Ker(π(k−1,k) ) |P (k) x|/|x|.
Lemma 4.18. It holds true that for k ∈ {2, . . . , q},
λmax (N

(k),T

N

(k)

)≤

1 + kP (k) k2Ker(π(k−1,k) )
λmin (W (k) W (k),T )

.

(4.38)


(k)
Proof. Since Im(W (k),T ) and Im(π (k,k−1) ) are orthogonal and dim(RI ) = dim Im(W (k),T ) +

(k)
(k)
(k−1)
dim Im(π (k,k−1) ) , for x ∈ RI there exists a unique y ∈ RJ
and z ∈ RI
such
(k),T
(k,k−1)
2
(k),T
2
(k,k−1)
2
(k)
that x = W
y+π
z and |x| = |W
y| + |π
z| . Observe that W x =
W (k) W (k),T y (since W (k) π (k,k−1) = 0) and R(k−1,k) x = R(k−1,k) W (k),T y + z (since
R(k−1,k) π (k,k−1) = I (k−1) from Theorem 4.8). Therefore, |x|2 = |W (k),T y|2 + |P (k) (x −
(k)
W (k),T y)|2 with y = (W (k) W (k),T )−1 W (k) x. Let v ∈ RJ . Taking x = A(k) W (k),T v and
observing that P (k) x = 0 (since R(k−1,k) A(k) W (k),T = 0 from the ·, · a -orthogonality
between V(k−1) and W(k) ) leads to |A(k) W (k),T v|2 = |W (k),T y|2 +|P (k) W (k),T y|2 with y =
(W (k) W (k),T )−1 B (k) v. Therefore |A(k) W (k),T v|2 ≤ (1+kP (k) k2Ker(π(k−1,k) ) ) λ

|B (k) v|2
(k) W (k),T ) ,

min (W

which concludes the proof after taking v = B (k),−1 v 0 and maximizing the l.h.s. over
|v 0 | = 1.

31

Lemma 4.19. Writing kM k2 := supx |M x|/x the spectral norm, we have
kP (k) k2Ker(π(k−1,k) ) ≤ kπ (k,k−1) A(k−1) π (k−1,k) k2

xT Θ(k) x
xT x
x∈Ker(π (k−1,k) )
sup

(4.39)

Proof. Let x ∈ Ker(π (k−1,k) ). Using P (k) = π (k,k−1) A(k−1) π (k−1,k) Θ(k) we obtain that
1
1
|P (k) x| = kπ (k,k−1) A(k−1) π (k−1,k) (Θ(k) ) 2 k2 |(Θ(k) ) 2 x|. Observing that for
1
1
M = π (k−1,k) (Θ(k) ) 2 we have M M T = Θ(k−1) and for N = π (k,k−1) A(k−1) π (k−1,k) (Θ(k) ) 2
we have λmax (N T N ) = λmax (N N T ) we deduce
1
kπ (k,k−1) A(k−1) π (k−1,k) (Θ(k) ) 2 k22 = kπ (k,k−1) A(k−1) π (k−1,k) k2 and conclude by taking the
supremum over x ∈ Ker(π (k−1,k) ).
Lemma 4.20. It holds true that
γ̄k2
xT Θ(k) x
2
≤
H
k−1
xT x
γk π 2 λmin (a)
x∈Ker(π (k−1,k) )
¯
sup

(4.40)

P
(k)
(k)
Proof. Let y ∈ RJ
and α ∈ R. Let x = αW (k),T y. Write φ = i∈I (k) xi φi and
ψ = (− div(a∇·))−1 φ. Observe that kψk2a = xT Θ(k) x ≥ αy T W (k) Θ(k) W (k),T y. UsR (k) (k)
(k)
ing Ω φi φl = 0 for i 6= l and selecting α = kφi k−2
(assuming, without loss
L2 (Ω)
(k)

(k)

of generality, that kφi k2L2 (Ω) = |τi | is constant in i, for the general case, rescale
R
(k)
(k−1)
=
each φi by a multiplicative constant) we obtain that for j ∈ I (k−1) , Ω φφj
R
P
(k) 2
(k−1,k)
= (π (k−1,k) W (k),T y)j = 0. Therefore, since kψk2a = Ω φψ,
i∈I (k) xi kφi kL2 (Ω) πj,i
R
(k−1)
we have for ψ 0 ∈ span{φi
| i ∈ I (k−1) } kψk2a = Ω φ(ψ − ψ 0 ) ≤ kφkL2 (Ω) kψ − ψ 0 kL2 (Ω) .
P
(k−1)
(k−1) 2
(k−1) R
/kφi
kL2 (Ω) we obtain (via Poincaré and
Choosing ψ 0 = i∈I (k−1) φi
Ω ψφi
Cauchy-Schwartz inequalities as in the proof of Proposition 3.6) that kψ − ψ 0 kL2 (Ω) ≤
p
p
Hk−1 kψka /(π λmin (a)) and deduce kψka ≤ Hk−1 kφkL2 (Ω) /(π λmin (a)). Observing
that kφk2L2 (Ω) ≤ |x|2 γ̄k and γk ≤ α−1 ≤ γ̄k we summarize and obtain that
2 ¯|x|2 γ̄ 2 /(π 2 λ
2
(k),T y|2 γ̄ 2 /(γ π 2 λ
y T W (k) Θ(k) W (k),T y ≤ Hk−1
min (a)), which
min (a)) ≤ Hk−1 |W
k
k
k
¯
(k−1,k)
(k),T
concludes the proof of the lemma (since Ker(π
) = Im(W
)).
Observing that kπ (k,k−1) A(k−1) π (k−1,k) k2 ≤ λmax (π (k,k−1) π (k−1,k) )λmax (A(k−1) ) and
using (4.35), we derive from lemmas 4.19 and 4.20 that
γ̄k2 25+d λmax (a)
.
(4.41)
γk γk−1 δ 2+d π 2 λmin (a)
¯ ¯
(k,k−1)
(k−1,k)
Observing that π
π
is block-diagonal and using
the notations of Definition 4.1
Pm
j
we have λmax (π (k,k−1) π (k−1,k) ) = maxj∈I (k−1) supx∈Rmj | i=1
xi |2 /|x|2 = maxj∈I (k−1) mj .
kP (k) k2Ker(π(k−1,k) ) ≤ λmax (π (k,k−1) π (k−1,k) )

(k−1)

Noting that a set τj
(k)

sets τi

(k−1)

can contain at most (maxj∈I (k−1) |τj

(k)

|)/(mini∈I (k) |τi |) sub-

we have
d
max mj ≤ Hk−1 /(δHk )

j∈I (k−1)

32

(4.42)

d
and conclude that λmax (π (k,k−1) π (k−1,k) ) ≤ Hk−1 /(δHk ) . Therefore (4.37) and Lemma
4.18 imply, after simplification, that
λmin (B (k) ) ≥

d
2+2d π 2 λ
λmin (a)
min (a)
(k)
(k),T Hk γk−1 γk δ
.
λ
(W
W
)
¯
¯
min
2
2
d
6+d
Hk−1 γ̄k
Hk−1 γ̄k 2 λmax (a)

(4.43)

Recalling that γ̄k /γk ≤ δ −d , using γ̄k /γk−1 ≤ Hkd /(Hk−1 δ)d , and summarizing we conclude the proof of¯ (4.33) and (4.34).¯ Recall that under construction 4.13 we have
W (k) W (k),T = J (k) which implies Cond(W (k) W (k),T ) = 1. Under construction 4.11,
W (k) W (k),T is block diagonal with for j ∈ I (k−1) , diagonal blocks corresponding to
(mj − 1) × (mj − 1) matrices M (mj −1) such that (1) for n = 1 and x ∈ R, xT M (1) x = 2x2
(2) x = x2 + (x − x )2 + x2 and (3) for for n ≥ 3,
(2) for n = 2 and x ∈ R2 , xT M P
2
1
1
2
2 + x2 . Note that for all n ≥ 1,
and x ∈ Rn , xT M (n) x = x21 + n−2
(x
−
x
)
i
i+1
n
i=1
λmax (M (n) ) ≤ 3. Furthermore, for n ≥ 3 (n ≤ 2 is trivial), introducing the variables
y2 = x2 − x1 , . . . , yn = xn − xn−1 we obtain that xT M (n) x = x21 + y22 + · · · + yn2 + x2n
and |x|2 = x21 + (x1 + y2 )2 + · · · + (x1 + y2 + · · · + yn )2 ≤ (xT M (n) x)n(n + 1)/2.
Therefore, λmin (M (n) ) ≥ 2/(n(n + 1)). We conclude that under construction 4.11
Cond(W (k) W (k),T ) ≤ maxj∈I (k−1) 3(mj − 1)mj /2 and bound mj as in (4.42).

4.9

Well conditioned relaxation across subscales

If Hk is a geometric sequence or if Hk−1 /Hk is uniformly bounded, then, by Theorem
(4.17), the linear systems ((4.26) and (4.27)) entering in the calculation of the gamblets
(k)
(k)
χi (and therefore ψi ) and the subband/subscale solutions u(1) and u(k+1) − u(k)
have uniformly bounded condition numbers (in particular, these condition numbers are
bounded independently from mesh size/resolution and the regularity of a(x)). Therefore
these systems can be solved efficiently using iterative methods. One such methods is the
Conjugate Gradient (CG) method [57]. Recall [98] that the application of the CG method
to a linear system Ax = b (where A is a n × n symmetric positive definite matrix) with
initial guess x(0) , yields√a sequence of approximations x(l) satisfying (writing |e|2A :=


Cond(A)−1 l
eT Ae) |x − x(l) |A ≤ 2 √
|x − x(0) |A where Cond(A) := λmax (A)/λmin (A).
Cond(A)+1

Recall [98] also that the maximum number of iterations
p required to2 reduce the error by
1
(l)
(0)
a factor  (|x − x |A ≤ |x − x |A ) is bounded by 2 Cond(A) ln  and has complexity
p
(number of required arithmetic operations) O( Cond(A)NA ) (writing NA the number
of non-zero entries of A).

4.10

Hierarchical localization and error propagation across scales

Although the multi-resolution decomposition presented in this section leads to well conditioned linear systems, the resulting matrices B (k) and A(k) are dense and to achieve
near-linear complexity in the resolution of (1.1) these matrices must be truncated by
(k)
(k)
localizing the computation of the basis functions ψi (and therefore χi ). The approximation error induced by these localization/truncation steps is controlled by the
33

exponential decay of gamblets and the uniform bound on the condition numbers of the
matrices B (k) . To make this control explicit and derive a bound the size of the localization sub-domains we need to quantify the propagation of truncation/localization errors
across scales and this is the purpose of this subsection.
For k ∈ {1, . . . , q}, ρ ≥ 1 and i ∈ I (k) we define (1) iρ as the subset of indices
(k)
(k)
j ∈ I (k) whose corresponding subdomains τj are at distance at most Hk ρ from τi
(k)

and (2) Sρi := ∪j∈iρ τj . Let ρ1 , . . . , ρq ≥ 1. For k ∈ {1, . . . , q − 1} and i ∈ I (k) ,
write iρ,k+1 as the subset of indices j ∈ I (k+1) such that j (k) ∈ iρ . For i ∈ I (q) let
(q+1),loc
(k),loc
Vi
:= H01 (Sρi q ). For k ∈ {1, . . . , q} and i ∈ I (k) , let ψi
be the minimizer of
Minimize kψka subject to ψ ∈

(k+1),loc
Vi

Z
and
Ω

(k)

ψφj

= δi,j for j ∈ iρk

(4.44)

(k+1),loc

where for k ∈ {1, . . . , q − 1} and i ∈ I (k) , Vi
is defined (via induction) by
(k+1),loc
(k+1),loc
ρ
,k+1
Vi
:= span{ψj
|j∈i k
}.
From now on we will assume that Hk = H k for some H ∈ (0, 1) (or simply that Hk
is uniformly bounded from below and above by H k ). To simplify the presentation, we
will also, from now on, write C any constant that depends only d, Ω, λmin (a), λmax (a), δ
(e.g., 2Cλmax (a) will still be written C). The following theorem allows us to control the
localization error propagation across scales. For k ∈ {1, . . . , q}, let A(k),loc be the I (k) ×
(k),loc
(k),loc
(k),loc
I (k) matrix defined by Ai,j
:= ψi
, ψj
and let E(k) be the (localization)
a
1

P
(k)
(k),loc 2 2
error E(k) :=
ka .
i∈I (k) kψi − ψi
Theorem 4.21. For k ∈ {1, . . . , q − 1}, we have
d

d

E(k) ≤ CH − 2 E(k + 1) + Ce−ρk /C H 2 −(k+1)(d+1) .

(4.45)

Proof. We will need the following lemma summarizing and simplifying some results
obtained in Theorem 4.16 when Hk = H k .
Lemma 4.22. Let Hk = H k and W (k) be as in Construction
4.11 or Construction 4.13.
√
It holds true that for k ∈ {q, . . . , 2} (1) kW (k) k2 ≤ 3 (2) 1/λmin (W (k) W (k),T ) ≤ CH −2d
d
(3) kπ̄ (k−1,k) k2 ≤ CH 2 (4) kπ (k−1,k) k2 ≤ CH −d/2 (5) kR(k−1,k) k2 ≤ CH −d/2 (6)
Cond(B (k) ) ≤ CH −2−2d (7) λmax (B (k) ) ≤ CH −k(2+d)
(8) 1/λmin (B (k) ) ≤ CH k(2+d)−2−2d . Furthermore, (9) Cond(A(1) ) ≤ CH −2 (10) 1/λmin (A(1) ) ≤
CH d and for k ∈ {1, . . . , q} (11) λmax (A(k) ) ≤ CH −k(2+d) .
Proof. From the proof of Theorem 4.17 we have (1) and 1/λmin (W (k) W (k),T ) ≤
(k−1,k)
maxi∈I (k−1) (mi −1)mi /2, which implies (2). For (3), noting that π̄i,j
= 0 if j (k−1) 6= i
d
√
(k−1,k)
and π̄i,j
= 1/mi otherwise, we have kπ̄ (k−1,k) k2 = maxi∈I (k−1) 1/ mi ≤ CH 2 . (4)
follows from λmax (π (k,k−1) π (k−1,k) ) = maxi∈I (k−1) mi ≤ CH −d . Let us now prove (5).
Using (4.23), (4.21) and defining N (k) = −A(k) W (k),T B (k),−1 as in Lemma 4.18, we have
R(k−1,k) = π̄ (k−1,k) + π̄ (k−1,k) N (k) W (k) , which leads to kR(k−1,k) k2 ≤ kπ̄ (k−1,k) k2 (1 +
34

kN (k) k2 kW (k) k2 ). Using Lemma
we obtain that λmax (N (k),T N (k) ) ≤
 4.18 and(k)(4.41)
(k,k−1)
(k−1,k)
d
(k),T
1 + Cλmax (π
π
)H /λmin (W W
) and therefore kN (k) k2 ≤ CH −d .
Summarizing we have obtained (5). (6), (7), (8) and (11) follow from Theorem 4.16
and in particular (4.36), (4.43) and (4.35). See (4.33) and the proof of Theorem 4.17 for
(9) and (10).
We will also need the following lemma.
Lemma 4.23. Let k ∈ {1, . . . , q − 1} and let R be the I (k) × I (k+1) matrix defined by
(k,k+1)
and Ri,j = 0 for j ∈ iρk ,k+1 and Ri,j = Ri,j
for j ∈ I (k+1) /iρk ,k+1 . It holds true that
kRk2 ≤ CH d/2 e−ρk /C .
P
(k,k+1) 2
| with |I (k) | ≤
Proof. Observe that kRk22 ≤ |I (k) | maxi∈I (k) j∈I (k+1) /iρk ,k+1 |Ri,j
CHk−d . Let i ∈ I (k) . Using Theorem 4.7 and Cauchy-Schwartz inequality we have
P
(k,k+1)
(k)
(k+1)
(k,k+1) 2
| ≤
|Ri,j
| ≤ kψi kL2 (τ (k+1) ) kφj
kL2 (τ (k+1) ) . Therefore j∈I (k+1) /iρk ,k+1 |Ri,j
j
j
P
P
(k)
(k)
2
2
d
CHk+1
j∈I (k+1) /iρk ,k+1 kψi kL2 (τ (k) ) =
j∈I (k+1) /iρk ,k+1 kψi kL2 (τ (k+1) ) . Observe that
j
j
R
P
(k) 2
(k)
= 0 for s 6= i we obtain from Poincaré’s ins∈I (k) /iρk kψi kL2 (τ (k) ) . Since τs(k) ψi
s
P
(k)
(k)
(k,k+1) 2
equality that kψi kL2 (τs ) ≤ Ck∇ψi kL2 (τ (k) ) Hk . Therefore j∈I (k+1) /iρk ,k+1 |Ri,j
| ≤
s
P
(k)
2
d H2
CHk+1
s∈I (k) /iρk k∇ψi kL2 (τ (k) ) . Using Theorem 3.11 we obtain that
k
s
P
(k) 2
−C −1 ρk kψ (k) k2 . Using (3.25) we have kψ (k) k2 ≤ CH −d−2 ,
k∇ψ
k
≤
Ce
ρ
(k)
(k)
a
a
s∈I /i k
i
i
i
k
L2 (τs )
P
−1 ρ
(k,k+1) 2
d
−C
k.
therefore j∈I (k+1) /iρk ,k+1 |Ri,j
| ≤ CH e
Let us now prove Theorem 4.21. We obtain by induction (using the constraints
(k),loc
in (4.44)) that for k ∈ {1, . . . , q} and i ∈ I (k) , ψi
satisfies the constraints of
(4.3). Moreover (3.8) implies that if ψ satisfies the constraints of (4.3) then kψk2a =
(k)
(k)
(k),loc
kψi k2a + kψ − ψi k2a . Therefore, for k ∈ {2, . . . , q − 1}, ψi
is also the miniP
(k)
(k+1),loc
mizer of kψ − ψi ka over functions ψ of the form ψ = j∈iρk ,k+1 cj ψj
satisfying
P
(k,k+1)
(k+1),loc
the constraints of (4.44). Thus, writing ψ ∗ :=
ψj
, we have
j∈iρk ,k+1 Ri,j
(k),loc

(k)

(k)

(since ψ ∗ satisfies the constraints of (4.44)) kψi
− ψi ka ≤ kψ ∗ − ψi ka . Write
P
P
(k,k+1) (k+1),loc
(k,k+1) (k+1),loc
ψ1 := j∈I (k+1) Ri,j
ψj
and ψ2 := j∈I (k+1) /iρk ,k+1 Ri,j
ψj
. Observ(k),loc

(k)

(k)

ing that ψ ∗ = ψ1 −ψ2 we deduce that kψi
−ψi k2a ≤ 2kψ1 −ψi k2a +2kψ2 k2a and after
2
P
P
(k,k+1)
(k+1)
summing over i, E(k) ≤ 2(I1 + I2 ) with I1 = i∈I (k) k j∈I (k+1) Ri,j
(ψj
−
P
P
(k+1),loc 2
(k,k+1) (k+1),loc 2
ψj
)ka and I2 =
ψj
ka . Writing S the
i∈I (k) k
j∈I (k+1) /iρk ,k+1 Ri,j
(k+1)

I (k+1) ×I (k+1) symmetric positive matrix with entries Si,j = ψi

1
2

(k+1),loc
ψj
,
a

(k+1),loc

−ψi

(k+1)

, ψj

−

note that I1 = Trace[R(k,k+1) SR(k+1,k) ]. Writing S the matrix square root
of S, observe that for a matrix U , using the cyclic property of the trace, Trace[U SU T ] =
1
1
Trace[S 2 U T U S 2 ] ≤ λmax (U T U ) Trace[S], which (observing that Trace[S] = (E(k + 1))2
2
and λmax (U T U ) = kU k22 ) implies I1 ≤ kR(k,k+1) k22 E(k + 1) . Therefore (using Lemma
35

√
d
4.22) we have I1 ≤ CH 2 E(k + 1). Let us now bound I2 . Let R be defined as
(k+1),loc
(k+1),loc
(k+1),loc
in Lemma 4.23. Noting that ψi
, ψj
= Ai,j
we have (as above)
a
I2 = Trace[RA(k+1),loc RT ] ≤ λmax (RT R) Trace[A(k+1),loc ]. Summarizing and using
p
d
d
+ 1) + CH 2 e−ρk /C Trace[A(k+1),loc ]. ObLemma 4.23 p
we deduce that E(k) ≤ CH 2 E(k p
serving that Trace[A(k+1),loc ] ≤ E(k + 1) + Trace[A(k+1) ] and using Trace[A(k+1) ] ≤
(k+1) 2
ka

−d
CHk+1
maxi∈I (k+1) kψi
proof of the theorem.

(k+1)

and (Lemma 3.14) kψi

− d −1

2
, we conclude the
ka ≤ CHk+1

(k),loc

Let u(1),loc be the finite element solution of (1.1) in V(1),loc := span{ψj
| j ∈
(1)
(k)
I }. For k ∈ {2, . . . , q}, let W
be defined as in Construction 4.11 or Construction
P
(k),loc
(k) (k),loc
4.13. For i ∈ J (k) , let χi
:= j∈I (k) Wi,j ψj
. For k ∈ {2, . . . , q} let u(k),loc −
(k),loc

u(k−1),loc be the finite element solution of (1.1) in W(k),loc := span{χj
| j ∈ J (k) }.
P
For k ∈ {2, . . . , q}, write u(k),loc := u(1),loc + kj=2 (u(j),loc − u(j−1),loc ). Let B (k),loc be
(k),loc

(k),loc

(k),loc

. Observe that B (k),loc =
P
(k)
(k),loc 2  21
ka .
W (k) A(k),loc W (k),T . Write for k ∈ {2, . . . , q}, E(k, χ) :=
j∈J (k) kχj − χj
The following theorem allows us to control the effect of the localization error on the
approximation of the solution of (1.1).
the J (k) × J (k) matrix defined by Bi,j

:= χi

, χj

a

Theorem 4.24. It holds true that for k ∈ {2, . . . , q} (1) E(k, χ) ≤ CH −d/2 E(k). Furthermore for k ∈ {2, . . . , q} and E(k, χ) ≤ C −1 H −k(1+d/2)+1+d we have
(2) Cond(B (k),loc ) ≤ CH −2−2d , and (3) ku(k) − u(k−1) − (u(k),loc − u(k−1),loc )ka ≤
CE(k, χ)kgkH −1 (Ω) H k(1+d/2)−3d−3 . Similarly for E(1) ≤ C −1 H −d/2 , we have
(4) Cond(A(1),loc ) ≤ CH −2 , and (5) ku(1) − u(1),loc ka ≤ CE(1)kgkH −1 (Ω) H −2+d/2 .
Proof. We will need the following lemma.
Lemma 4.25. Let χ1 , . . . , χm be linearly independent elements of H01 (Ω). LetP
χ01 , . . . , χ0m
m
1
be another set of linearly independent elements of H0 (Ω). Write E :=
i=1 kχi −
1

χ0i k2a 2 . Let B (resp. B 0 ) be the m × m matrix defined by Bi,j = χi , χj a (resp.
0
Bi,j
= χ0i , χ0j a ). Let um (resp. u0m ) be the solution of (1.1) in span{χi | i =
p
1, . . . , m} (resp. span{χ0i | i = 1, . . . , m}). It holds
true that for E ≤ λmin (B)/2
p
(1) Cond(B 0 ) ≤ 8 Cond(B) (2) kB − B 0 k2 ≤ 3 λmax (B)E (3) kB −1 − (B 0 )−1 k2 ≤
p
−2
Cond(B)
12 λmax (B) λmin (B) E and (4) kum − u0m ka ≤ CEkgkH −1 (Ω) √
.
λmin (B)

p
p
P
Proof. For (1) observe that
λmax (B 0 ) = sup|x|=1 k m
xi χ0i ka ≤
λmax (B) + E
i=1
p
p
Pm
0
0
λmin (B) − E. For (2) observe that for
and λmin (B ) = inf |x|=1 k i=1 xi χi ka ≥
Pm
Pm
0
m
T
x, y ∈ R with |x| = |y| = 1 we have y (B − B 0 )x =
i=1 yi (χi − χi ),
i=1 xi χi a −
p
p
Pm
Pm
0
0
0
(2) and
i=1 yi χi ,
i=1 xi (χi − χi ) a ≤ ( λmax (B ) + λmax (B))E. (3) follows from
P
kB −1 − (B 0 )−1P
k2 ≤ kB − B 0 k2 / λmin (B)λmin (B 0 ) . ForR(4) observe that um = m
i=1 wi χi
m
0
0
0
−1
0
0
−1
(resp. um = i=1 wi χi ) where w = B b with bi = Ω gχi (resp. w = (B ) b0 with
p
R
b0i = Ω gχ0i ). Therefore kum − u0m ka ≤ |w|E + |w − w0 | λmax (B). w − w0 = B −1 (b − b0 ) −
36

0
0
0
B −1 (B − B 0 )w0 leads
H −1 (Ω) E + kB − B k2 |w |)/λmin (B). Using (2),
Pm to |w −2w | ≤ C(kgk
2
2
2
λmin (B)|w| ≤ k i=1 wi χi ka ≤ kuka ≤ CkgkH −1 (Ω) , and λmin (B 0 )|w0 |2 ≤ Ckgk2H −1 (Ω)
we conclude the proof of (4) after simplification.
P
(k)
(k),loc
(k)
(k)
(k),loc
Let us now prove Theorem 4.24. Using χj − χj
= i∈I (k) Wj,i (ψi − ψi
)
2 P
(k)
(k−1)
(k−1)
and noting that Wj,i = 0 for i
6= j
we have E(k, χ) ≤ j∈J (k)
2
P
(k) 2 P
(k)
(k),loc 2 
ka . Therefore, E(k, χ) ≤
i∈I (k) (Wj,i )
i∈I (k) ,i(k−1) =j (k−1) kψi − ψi
2
P
(k)
E(k) maxi∈I (k−1) mi maxj∈J (k) i∈I (k) (Wj,i )2 . Observing that (see (4.42))
P
(k)
maxi∈I (k−1) mi ≤ 1/(Hδ)d and i∈I (k) (Wj,i )2 ≤ λmax (W (k) W (k),T ) ≤ 3 (see Lemma
1

4.22) we conclude that (1) holds true with C = (3/δ d ) 2 . (2) and (3) are a direct
application of lemmas 4.25 and 4.22. For (3), observe that u(k) − u(k−1) (resp. u(k),loc −
(k),loc
u(k−1),loc ) is the finite element solution of (1.1) in W(k) (resp. W(k),loc := span{χj
|
(k)
j ∈ J }). The proof of (4) and (5) is similar to that of (2) and (3).
P
d
d
d
Theorem 4.26. Let k ∈ {1, . . . , q}. We have, E(k) ≤ C qj=k e−ρj /C C j−k H − 2 +k 2 −j3 2 .
Proof. By Theorem 4.21, for k ∈ {1, . . . , q −1}, we have E(k) ≤ ak +bk E(k +1) with ak =
d
d
Ce−ρk /C H 2 −(k+1)(d+1) and bk = CH − 2 . Therefore we obtain by induction that E(k) ≤
ak + bk ak+1 + bk bk+1 ak+2 + · · · + bk · · · bq−2 aq−1 + bk · · · bq−1 E(q). Using Theorem 3.13
we have E(q) ≤ CH −d/2−q(2+d/2) e−ρq /C and obtain the result after simplification.

1
Theorem 4.27. Let  ∈ (0, 1). It holds true that if ρk ≥ C (1 + ln(1/H)
) ln H1k + ln 1
for k ∈ {1, . . . , q} then (1) for k ∈ {1, . . . , q} we have ku(k) − u(k),loc ka ≤ kgkH −1 (Ω)
and ku − u(k),loc ka ≤ C(H k + )kgkL2 (Ω) (2) Cond(A(1),loc ) ≤ CH −2 , and for k ∈
{2, . . . , q} we have (3) Cond(B (k),loc ) ≤ CH −2−2d and (4) ku(k) − u(k−1) − (u(k),loc −
u(k−1),loc )ka ≤ 2k 2 kgkH −1 (Ω) .
Proof. Theorems 4.3 and 4.24 imply that the results of Theorem 4.27 hold true if for
k ∈ {1, . . . , q} E(k) ≤ C −1 H −k(1+d/2)+7d/2+3 /k 2 . Using Theorem 4.26 we deduce that
the results of Theorem 4.27 hold true if for k ∈ {1, . . . , q} and k ≤ j ≤ q we have
d
d
d
Ce−ρj /C C j−k H − 2 +k 2 −j3 2 ≤ H −k(1+d/2)+7d/2+3 /(k 2 j 2 ). We conclude after simplification.

5
5.1

The algorithm, its implementation and complexity
The initialisation of the algorithm

To describe the practical implementation of the algorithm we consider the (finite-element)
discretized version of (1.1). Let Th be a regular fine mesh discretization of Ω of resolution
h with 0 < h  1. Let N be the set of interior nodes zi and N = |N | be the number of
interior nodes (N = O(h−d )) of Th . Write (ϕi )i∈N a set of regular nodal basis elements

37

(of H01 (Ω)) constructed from Th such that for each i ∈ N , support(ϕi ) ⊂ B(zi , C0 h) and
for y ∈ RN ,
X
(5.1)
γhd |y|2 ≤ k
yi ϕi k2L2 (Ω) ≤ γ̄hd |y|2
¯
i∈N
for some constants γ, γ̄, C0 ≈ O(1). In addition to (5.1) the regularity of the finite
¯
elements is used to ensure
the availability of the inverse Poincaré inequality
k∇vkL2 (Ω) ≤ C1 h−1 kvkL2 (Ω)

(5.2)

for v ∈ span{ϕi | i ∈ N } and some constant C1 ≈ O(1), used to generalize the proof of
Theorem 4.16 P
to the discrete case.
Given g = i∈N gi ϕi we want to find u ∈ span{ϕi | i ∈ N } such that for all j ∈ N ,
Z
ϕj g for all j ∈ N
(5.3)
ϕj , u a =
Ω

In practical applications a is naturally assumed to be piecewise constant over the fine
mesh (e.g. of constant value in each triangle or square of Th ) and one purpose of the
algorithm is the fast resolution of the linear system (5.3) up to accuracy  ∈ (0, 1).

Figure 2: The (fine) mesh Th , a (in log10 scale) and u.

Example 5.1. We will illustrate the presentation of the algorithm with a numerical
example in which Th is a square grid of mesh size h = (1 + 2q )−1 with q = 6 and
64 × 64 interior nodes (Figure 2). a is piecewise constant
on each square of Th and


Q6
j
i
k
given by a(x) = k=1 1 + 0.5 cos 2 π( 2q +1 + 2q +1 )
1 + 0.5 sin 2k π( 2qj+1 − 3 2q i+1 )
j
j+1
for x ∈ [ 2q i+1 , 2i+1
q +1 ) × [ 2q +1 , 2q +1 ). The contrast of a (i.e., when a is scalar, the ratio
between its maximum and minimum value) is 1866. The finite-element discretization
(5.3) is obtained using continuous nodal basis elements ϕi spanned by {1, x1 , x2 , x1 x2 } in
each square of Th . Writing zP
for our
i the positions of the interior nodes of Th , we choose, 
numerical example, g(x) = i∈N cos(3zi,1 + zi,2 ) + sin(3zi,2 ) + sin(7zi,1 − 5zi,2 ) ϕi (x).

38

Figure 3: I (1) , I (2) and I (3) .

The first step of the proposed algorithm is the construction of the index tree I
of Definition 4.2 describing the domain decomposition of Definition 4.2. To ensure
a uniform bound on the condition numbers of the stiffness matrices (4.19) one must
select the resolutions Hk to form a geometric sequence (or simply such that Hk−1 /Hk
is uniformly bounded), i.e. Hk = H k for some H ∈ (0, 1) (for our numerical example
H = 1/2, q = 6 and we identify I (k) as the indices of the interior nodes of a square
grid of resolution (1 + 2k )−1 as illustrated in Figure 3). In this construction H q = h
(q)
corresponds to the resolution of the fine mesh and each subset τi (i ∈ I (q) ) contains
one and only one element of N (interior node of the fine mesh). Using this one to one
correspondence we use the elements of I = I q to (re)label the nodal elements (ϕi )i∈N
(k)
as (ϕi )i∈I . The measurement functions (φi )i∈I (k) are then identified (1) by selecting
(q)
φi = ϕi for i ∈ I (q) and (2) via the nested aggregation (4.1) of the nodal elements (as
P
P
(k)
(k,k+1) (k+1)
(k+1)
commonly done in AMG), i.e. φi = j∈I (k+1) πi,j
φj
= j∈i(k,k+1) φj
for
(k)
k ∈ {1, . . . , q − 1} and i ∈ I .
Remark 5.1. We refer to Figure 4 for an illustration of these measurement functions
(k)
for our numerical example. Note that the support of each φi is only approximatively
(k)
(k)
(and not exactly) τi and that the φi are only approximate set functions (and not
exact ones). This does not affect the design, accuracy and localization of the algorithm presented here because the frame inequalities (4.32), and the Poincaré inequalities
P
(k)
(k)
k i∈I (k) xi φi kH −1 (Ω) ≤ C H k−1 kφi kL2 (Ω) for x ∈ Ker(π (k−1,k) ), hold true. Indeed,
(5.1) and Construction 4.2 imply that the frame inequalities (4.32) with γ̄k ≤ γ̄δ −d and
γk ≥ γδ d , and the Poincaré inequalities are regularity/homogeneity conditions on the
¯mesh ¯and the aggregated elements. Although a fine mesh has been used to facilitate
the presentation of the algorithm, the proposed method is meshless (it only requires the
specification of the basis elements (ϕi )i∈I ).

39

Figure 4: The functions φki with k ∈ {1, . . . , q} and q = 6.

5.2

Exact gamblet transform and multiresolution operator inversion

The near-linear complexity of the proposed multi-resolution algorithm (Algorithm 2) is
based on three properties (i) nesting (ii) uniformly bounded condition numbers (iii) localization/truncation based on exponential decay. Truncation/localization levels/subsets
are, a priori, functions of the desired level of accuracy  ∈ (0, 1) in approximating the
solution of (5.3) and to distinguish the implementation of localization/truncation (and
its consequences) we will first describe this algorithm in its zero approximation error
version (i.e.  = 0 and without using localization/truncation, Algorithm 1). Although
this error-free version (Algorithm 1) performs the decomposition of the resolution of
the linear system (5.3) (whose condition number is of the order of h−d−2  1) into the
resolutions of a nesting of linear systems with uniformly bounded condition numbers,
it is not of near linear complexity due to the presence of dense matrices. Algorithm 2
achieves near-linear complexity by truncating/localizing the dense matrices appearing
in Algorithm 1 (-accuracy is ensured using the off-diagonal exponential decay of these
dense matrices). Let us now describe Algorithm 1 in detail. Lines 1 and 2 correspond to
the computation of the (sparse) mass and stiffness matrices of (5.3). Line 4 corresponds
(q)
to the calculation of level q gamblets ψi defined as the minimizer of kψka subject to
R
(q)
= δi,j and ψ ∈ span{ϕl | l ∈ I}, note that since the number of constraints
Ω ψφj
R
(q)
is equal to the number of degrees of freedom of ψ, and since Ω ϕl φj = Ml,j , level q
gamblets do not depend on a and are obtained by inverting the mass matrix in Line 3
40

Algorithm 1 Exact Gamblet transform/solve.
R
1: For i, j ∈ I (q) , Mi,j = Ω ϕi ϕj
R
2: For i, j ∈ I (q) , Ai,j = Ω (∇ϕi )T a∇ϕj
3: Compute M −1
P
(q)
−1
4: For i ∈ I (q) , ψi =
j∈I (q) Mi,j ϕj
(q)

(q)

5:

For i ∈ I (q) , gi

6:

For i, j ∈ I (q) , Ai,j = Ω (∇ψi )T a∇ψj
for k = q to 2 do
B (k) = W (k) A(k) W (k),T
w(k) = B (k),−1 W (k) g (k)
P
(k)
(k) (k)
For i ∈ J (k) , χi = j∈I (k) Wi,j ψj
P
(k) (k)
u(k) − u(k−1) = i∈J (k) wi χi
D(k,k−1) = −B (k),−1 W (k) A(k) π̄ (k,k−1)
R(k−1,k) = π̄ (k−1,k) + D(k−1,k) W (k)
A(k−1) = R(k−1,k) A(k) R(k,k−1)
P
(k−1)
(k−1,k) (k)
For i ∈ I (k−1) , ψi
= j∈I (k) Ri,j
ψj
(k−1)
(k−1,k)
(k)
g
=R
g
end for
U (1) = A(1),−1 g (1)
P
(1) (1)
u(1) = i∈I (1) Ui ψi
u = u(1) + (u(2) − u(1) ) + · · · + (u(q) − u(q−1) )

7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:

= gi
(q)

// gi
R

(q)

(q)

// Mass matrix
// Stiffness matrix
// Mass matrix inversion
// Level q gamblets
R (q)
P
= Ω ψi g with g = i∈I (q) gi ϕi
// A(q) = M −1 AM −1,T
// Eq. (4.19)
// Eq. (4.26)
// Eq. (4.17)
// Thm. 4.15
// Eq. (4.21)
// Eq. (4.23)
// Eq. (4.16)
// Eq. (4.11)
// Eq. (4.25)
// Eq. (4.27)
// Thm. 4.15
// Thm. 4.6 with u = u(q)

(note that by (5.1), the mass matrix is of O(1) condition number). Although not done
(q)
here, one can also initialize the algorithm (and its fast version) with ψi = ϕi (which
P
−1 (q)
is equivalent to using j∈I (q) Mi,j
ϕj as level q measurement functions). Line 5 corresponds to initialization of the vector g (q) introduced above (4.25). Line 6 corresponds
to the initialization of the stiffness matrix A(q) introduced in (4.14). The core of the
algorithm is the nested computation performed (iteratively from k = q down to k = 2)
(k)
in lines 8 to 16. Note that this nested computation takes A(k) , g (k) and (ψi )i∈I (k) as
(k−1)
inputs and produces (1) A(k−1) , g (k−1) and (ψi
)i∈I (k) as outputs for the next itera(k)
(k)
(k−1)
tion and (2) the subband u − u
of the solution and subband gamblets (χi )i∈J (k)
(which, do not need to be explicitly computed/stored since Line 11 is equivalent to
P
(k)
(k)
u(k) − u(k−1) = i∈I (k) (W (k),T w(k) )i ψi ). Note also that the gamblets (ψi )i∈I (k) and
(k)

(χi )i∈J (k) can be stored and displayed using the hierarchical structure (4.11). Through
this section and the remaining part of the paper we assume that the matrices W (k) are
obtained as in Construction 4.11 or 4.13. Note that the number of non-zero entries of
π (k−1,k) and W (k) is O(|I (k) |) (proportional to H −k in our numerical example). Lines 9
corresponds to solving the well conditioned linear system B (k) w(k) = W (k) g (k) and the
|I (k−1) | well conditioned linear systems B (k) D(k,k−1) = −W (k) A(k) π (k,k−1) . Note that by
41

Theorem 4.17 the matrices B (k) have uniformly bounded condition numbers and these
linear systems can be solved efficiently using iterative methods (such as the Conjugate
Gradient method recalled in Subsection (4.9)). u(1) is computed in lines 19 and 20 (recall that A(1) is also of uniformly bounded condition number) and the last step of the
algorithm, is to obtain u via simple addition of the subband/subscale solution u(1) and
(u(k) − u(k−1) )2≤k≤q . Observe that the operating diagram of Algorithm 1 is not a V or
(k)
W but an inverted pyramid (or a comb). More precisely, the basis functions ψi are
computed hierarchically from fine to coarse scales. Furthermore as soon as the elements
(k)
ψi have been computed, they can be applied (independently from the other scales) to
the computation of u(k) − u(k−1) (the projection of u onto W(k) corresponding to the
bandwidth [H k , H k−1 ]).

Figure 5: The basis elements ψik with k ∈ {1, . . . , 6}.

(k)

Example 5.2. We refer to figures 5 and 7 for an illustration of the gamblets ψi and
(k)
χj corresponding to Example 5.1 with W (k) defined by Construction 4.11. We refer
(k)

to Figure 6 for an illustration of the exponential decay of the gamblets ψi . We refer
to Figure 8 for an illustration of the condition numbers of A(k) and B (k) (with W (k)
still defined by Construction 4.11). Observe that the bound on the condition numbers
42

Figure 6: Exponential Decay.

Figure 7: The basis elements ψi1 and χki with k ∈ {2, . . . , 6}.

of B (k) depends on the contrast and the saturation of that bound occurs for smaller
values of k under low contrast. We refer to Figure 9 for an illustration of the subband
solutions u(1) , u(2) − u(1) , . . . , u(q) − u(q−1) corresponding to Example 5.1. Observe that
these (subband) solutions form a multiresolution decomposition of u as a sum of functions
characterising the behavior of u at subscales [H, 1], [H 2 , H],. . . ,[H q , H q−1 ]. Once the
components u(1) , u(2) − u(1) ,. . . , and u(q) − u(q−1) have been computed one obtains, via
simple summation, u(1) , . . . , u(q) , the finite-element approximation of u at resolutions
H, H 2 , . . . , H q illustrated in Figure 11. As described in Theorem 4.3 the error of the
approximation of u by u(k) is proportional to H k for k ∈ {1, . . . , q − 1}. For k = q, as

43

Figure 8: Condition numbers of A(k) and B (k) .

Figure 9: u(1) , u(2) − u(1) ,. . . , and u(q) − u(q−1) .

illustrated in Figure 11, this approximation error drops down to zero because there is no
(q)
gap between H q and the fine mesh (i.e., ψi and ϕi span the same linear space in the
discrete case). Moreover, as illustrated in Figure 10, the representation of u in the basis
(k)

(1)

formed by the functions

ψi

(1)
kψi ka

and

χj

(k)
kχj ka

is sparse. Therefore, as illustrated in Figure

11 one can compress u, in this basis, by setting the smallest coefficients to zero without
loss in energy norm.

44

(1)

Figure 10: The coefficients of u in the expansion u =

(1) ψi
i ci kψ (1) k

P

i

5.3

a

+

Pq

k=2

(k)

(k) χj
j cj kχ(k) k

P

j

.

a

Fast Gamblet transform/solve

Algorithm 2 achieves near linear complexity (1) in approximating the solution of (5.3)
to a given level of accuracy  and (2) in performing an approximate Gamblet transform
(sufficient to achieve that level of accuracy). This fast algorithm is obtained by localizing/truncating the linear systems corresponding to lines 3 and 12 in Algorithm 1. We
define these localization/truncation steps as follows. For k ∈ {1, . . . , q} and i ∈ I (k)
define iρ as in Subsection 4.10 (i.e. as the subset of indices j ∈ I (k) whose corresponding
(k)
(k)
subdomains τj are at distance at most Hk ρ from τi ).
(i,ρq )

Definition 5.2. For i ∈ I (q) , let M (i,ρq ) be the iρq × iρq matrix defined by Ml,j

= Ml,j
(i,ρ
)
for l, j ∈ iρq . Let e(i,ρq ) be the |iρq |-dimensional vector defined by ej q = δj,i for j ∈ iρq .
Let y (i,ρq ) be the |iρq |-dimensional vector solution of M (i,ρq ) y (i,ρq ) = e(i,ρq ) . We define the
−1,ρq

solution M·,i
2) as the

−1,ρq

of the localized linear system M i,ρq M·,i

iρq -vector

given by

−1,loc
Mj,i

=

(i,ρ )
yj q

for j ∈

= δ·,i (Line 3 of Algorithm

iρq .

(q),loc

Note that the associated gamblet ψi
(Line 4 of Algorithm
2) is also the solution
R
of the problem of finding ψ ∈ span{ϕj | j ∈ iρq } such that Ω ψϕj = δi,j for j ∈ iρq (i.e.
(q)

localizing the computation of the gamblet ψi to a subdomain of size Hq ρq ). Line 5
R (q),loc
(q),loc
(q),loc
can be replaced by gi
= gi without loss of accuracy (gi
= Ω ψi
g simplifies
45

Figure 11: u(1) , . . . , u(q) . Relative approximation error in energy norm in log10 scale.
(1)
(2)
(q)
Compression of u over the basis functions ψi , χi , . . . , χi by setting 99% of the smallest coefficients to zero in the decomposition of Figure 10.

the presentation of the analysis). Line 12 of Algorithm 2 is defined in a similar way as
follows.
Definition 5.3. Let k ∈ {2, . . . , q} and B be the positive definite J (k) × J (k) matrix
B (k),loc computed in Line 8 of Algorithm 2. For i ∈ I (k−1) , let ρ = ρk−1 and let iχ
be the subset of indices j ∈ J (k) such that j (k−1) ∈ iρ (recall that if j = (i1 , . . . , ik ) ∈
J (k) then j (k−1) := (j1 , . . . , jk−1 ) ∈ I (k−1) ). B (i,ρ) be the iχ × iχ matrix defined by
(i,ρ)
(i,ρ)
Bl,j = Bl,j for l, j ∈ iχ . Let b(i,ρ) be the |iχ |-dimensional vector defined by bj
=
(k)
(k),loc
(k,k−1)
χ
(i,ρ)
χ
−(W A
π̄
)j,i for j ∈ i . Let y
be the |i |-dimensional vector solution
of B (i,ρ) y (i,ρ) = b(i,ρ) . We define the solution D(k,k−1),loc of the localized linear system
Inv(B (k),loc D(k,k−1),loc = −W (k) A(k),loc π̄ (k,k−1) , ρk−1 ) as the J (k) × I (k−1) sparse matrix
(k,k−1),loc
(k,k−1),loc
(i,ρ)
given by Dj,i
= 0 for j 6∈ iχ and Dj,i
= yj
for j ∈ iχ . D(k−1,k),loc (Line
13 of Algorithm 2) is then defined as the transpose of D(k,k−1),loc .
46

Algorithm 2 Fast Gamblet transform/solve.
R
1: For i, j ∈ I (q) , Mi,j = Ω ϕi ϕj
// O(N )
R
2: For i, j ∈ I (q) , Ai,j = Ω (∇ϕi )T a∇ϕj
// O(N )

−1,ρ
q
i,ρ
d
3: M q M·,i
= δ·,i
// Def. 5.2, Thm. 5.6, O N ρq ln max( 1 , q)
P
−1,ρ
(q),loc
4: For i ∈ I (q) , ψi
= j∈iρq Mj,i q ϕj
// O(N ρdq )
R
(q),loc
(q),loc
5: For i ∈ I (q) , gi
= Ω ψi
g
// O(N ρdq )
R
(q),loc
(q),loc T
(q),loc
6: For i, j ∈ I (q) , Ai,j
= Ω (∇ψi
) a∇ψj
// O(N ρ2d
q )
7: for k = q to 2 do
8:
B (k),loc = W (k) A(k),loc W (k),T
// O(|I (k) |ρdk )
9:
w(k),loc = (B (k),loc )−1 W (k) g (k),loc
// Thm. 5.6, O(|I (k) |ρdk ln 1 )
P
(k),loc
(k)
(k),loc
10:
For i ∈ J (k) , χi
= j∈I (k) Wi,j ψj
// O(|I (k) |ρdk )
P
(k),loc (k),loc
11:
u(k),loc − u(k−1),loc = i∈J (k) wi
χi
// O(N ρdk )
12:
Inv(B (k),loc D(k,k−1),loc = −W (k) A(k),loc π̄ (k,k−1) , ρk−1 )
// Def. 5.3, Thm. 5.6,
O(|I (k) |ρdk−1 ρdk ln 1 )
13:
R(k−1,k),loc = π̄ (k−1,k) + D(k−1,k),loc W (k)
// Def. 5.3, O(|I (k−1) |ρdk−1 )
d
14:
A(k−1),loc = R(k−1,k),loc A(k),loc R(k,k−1),loc
// O(|I (k−1) |ρ2d
k−1 ρk )
P
(k−1),loc
(k−1,k),loc (k),loc
15:
For i ∈ I (k−1) , ψi
= j∈I (k) Ri,j
ψj
// O(|I (k−1) |ρdk−1 )
16:
g (k−1),loc = R(k−1,k),loc g (k),loc
// O(|I (k−1) |ρdk−1 )
17: end for
18: U (1),loc = A(1),loc,−1 g (1),loc
// O(|I (1) |ρd1 ln 1 )
P
(1)
(1)
19: u(1),loc =
// O(N ρd1 )
i∈I (1) Ui ψi
loc
(1),loc
(2),loc
(1),loc
(q),loc
(q−1),loc
20: u
=u
+ (u
−u
) + · · · + (u
−u
)
// O(N q)

Remark 5.4. Definition 5.3 (Line 8 of Algorithm 2) is equivalent to localizing the
(k−1)
computation of each gamblet ψi
to a subdomain of size Hk−1 ρk−1 , i.e., the gamblet
(k−1),loc
ψi
computed in Line 15 of Algorithm 2 is the solution of (1) the problem of
P
(k−1,k) (k),loc
(k),loc
| j (k−1) ∈ iρk−1 }
finding ψ in the affine space
ψj
+ span{χj
j∈I (k) π̄i,j
such that ψ is ·, ·

(k),loc

a

(k),loc

minimizing kψka in span{ψl
for j ∈ iρk−1 .

5.4

| j (k−1) ∈ iρk−1 }, and (2) the problem of
R (k−1)
∈ iρk−1 } subject to constraints Ω φj
ψ = δi,j

orthogonal to span{χj
| l(k−1)

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

The sizes of the localization radii ρk (and therefore the complexity of Algorithm 2)
depend on whether Algorithm 2 is used as a pre-conditioner (as done with AMG) or
as a direct solver. Although it is natural to expect the complexity of Algorithm 2
to be significantly smaller if used as pre-conditioner (since pre-conditioning requires
lower accuracy and therefore smaller localization radii) we will restrict our analysis
47

and presentation to using Algorithm 2 as a direct solver. Note that, when used as
a direct solver, Algorithm 2 is parallel both in space (via localization) and in band(k−1),loc
with/subscale (subscales can be computed independently from each other and ψi
and u(k),loc − u(k−1),loc can be resolved in parallel). We will base our analysis on the
results of Subsection 4.10 and in particular Theorem 4.27. Although obtained in a
continuous setting, these results can be generalized to the discrete setting without difficulty. Two small differences are worth mentioning. (1) In this discrete setting, an
alternative approach for obtaining localization error bounds in the first step of the al(q),loc
gorithm (the computation of the localized gamblets ψi
) is to use the exponential
decay property of the inverse of symmetric well-conditioned banded matrices [26]: since
M is banded and of uniformly bounded condition number [26] (see also [11, Thm 4.10])

(q) (q)
−1
implies that Mi,j
decays like exp − dist(τi , τj )/C which guarantees that the bound
E(q) ≤ CH −d/2−q(2+d/2) e−ρq /C (used in Theorem 4.26) remains valid in the discrete
setting. (2) Since the basis functions ϕi are not exact set functions, neither are the
R (k),loc (k)
(k)
resulting aggregates φi . This implies that, in the discrete setting, Ω ψi
φj is not
(k)

necessarily equal to zero if τj is adjacent to Sρi k (with j 6∈ iρk , using the notation of
Subsection 4.10). This, however does not prevent the generalization of the proof because
R (k),loc (k)
(k)
the value of Ω ψi
φj (when τj is adjacent to Sρi k ) can be controlled via the exponential decay of the basis functions (e.g. as done in the proof of Theorem 3.13). We will
summarize this generalization in the following theorem (where the constant C depends
on the constants C1 , C0 , γ̄ and γ associated with the finite elements (ϕi ) in (5.1), in
¯ δ).
addition to d, Ω, λmin (a), λmax (a),
Theorem 5.5. Let u be the solution of the discrete system (5.3). Let u(1),loc , u(k),loc −
u(k−1),loc , uloc , A(k),loc and B (k),loc be the outputs of Algorithm 2. Let u(1) and u(k) −
(k−1) be the outputs of Algorithm 1. For k ∈ {2, . . . , q}, write u(k),loc := u(1),loc +
u
Pk
1
1
(j),loc −u(j−1),loc ). Let  ∈ (0, 1). It holds true that if ρ ≥ C (1+
k
j=2 (u
ln(1/H) ) ln H k +

ln 1 for k ∈ {1, . . . , q} then (1) for k ∈ {1, . . . , q − 1} we have ku(k) − u(k),loc kH01 (Ω) ≤
kgkH −1 (Ω) and ku(k) − u(k),loc kH01 (Ω) ≤ C(H k + )kgkL2 (Ω) (2) Cond(A(1),loc ) ≤ CH −2 ,
and for k ∈ {2, . . . , q} we have (3) Cond(B (k),loc ) ≤ CH −2−2d and (4) ku(k) − u(k−1) −
(u(k),loc − u(k−1),loc )kH01 (Ω) ≤ 2k 2 kgkH −1 (Ω) . Finally, (6) ku − uloc kH01 (Ω) ≤ kgkH −1 (Ω) .
Therefore, according to Theorem 5.5 if the localization radii ρk are chosen so that
ρk = O ln max(1/, 1/Hk ) for k ∈ {1, . . . , q} then the condition numbers of the matrices
B (k),loc and A(1),loc remain uniformly bounded and the algorithm achieves accuracy  in
a direct solve. The following theorem shows that the linear systems appearing in lines
3, 9 and 12 of Algorithm 2 do not need to be solved exactly and provide bounds on the
accuracy requirements (to simplify notations, we will from now on drop the superscripts
of the vectors y and b appearing in definitions 5.2 and 5.3).
1
Theorem 5.6. The results of Theorem 5.5 remain true if (1) ρk ≥ C (1+ ln(1/H)
) ln H1k +

ln 1 for k ∈ {1, . . . , q} (2) For each i ∈ I (q) the localized linear system M i,ρq y = δ·,i
of Definition 5.2 and Line 3 of Algorithm 2 is solved up to accuracy |y − y app |M i,ρq ≤

48

C −1 H 7d/2+3 /q 2 (using the notations of Subsection 4.9, i.e. |e|2A := eT Ae, and writing
y app the approximation of y) (3) For k ∈ {2, . . . , q} and each i ∈ I (k−1) , the localized
linear system B (i,ρ) y = b of Definition 5.3 and Line 12 of Algorithm 2 is solved up
to accuracy |y − y app |B (i,ρ) ≤ C −1 H −k+7d/2+4 /(k − 1)2 . (4) For k ∈ {2, . . . , q} the
linear system B (k),loc y = W (k) g (k),loc of Line 9 of Algorithm 2 is solved up to accuracy
|y − y app |B (k),loc ≤ kgkH −1 (Ω) /(2q).
Proof. From the proof of Theorem (4.27) we need E(k) ≤ C −1 H −k(1+d/2)+7d/2+3 /k 2
for k ∈ {1, . . . , q}. By the inverse Poincaré inequality (5.2) this inequality is satisfied
(q)
(q),loc
for k = q for kψi − ψi
kL2 (Ω) ≤ C −1 H 7d/2+3 /q 2 for each i ∈ I (q) , which by
i,ρ
the definition of M q and Line 4 of Algorithm 2 leads to (2). For k ∈ {2, . . . , q}
the inequality E(k − 1) ≤ C −1 H −(k−1)(1+d/2)+7d/2+3 /(k − 1)2 is satisfied if for i ∈
(k−1)
(k−1),loc
I (k−1) , kψi
− ψi
ka ≤ C −1 H −(k−1)+7d/2+3 /(k − 1)2 . Using the notations of
P
P
(k−1),loc
(k−1,k) (k),loc
(k−1,k),loc (k),loc
=
ψj
+ j∈iχ Di,j
χj
Definition 5.3 we have, ψi
j∈I (k) π̄i,j
(k),loc

with χj

(k),loc

, χl

(i,ρ)

which leads to (3) by lines 15, 13, 10 and 8 of Algorithm
P
(k),loc
2. For (4) we simply observe that for y ∈ J (k) , k i∈J (k) (y − y app )i χi
ka = |y −
app
y |B (k),loc .
a

= Bj,l

Let us now describe the complexity of Algorithm 2. This complexity depends on
the desired accuracy  ∈ (0, 1). Lines 1 and 2 correspond to the computation of the
(sparse) mass and stiffness matrices of (5.3). Note since A and M are sparse and
banded (of bandwidth 2d = 4 in our numerical example) this computation is of O(N )
complexity. Line 3 corresponds to the resolution of the localized linear system introduced
in Definition 5.2 using M i,ρq , the iρq × iρq sub-matrix of M . According to Theorem 5.5,
the accuracy of each solve must be |y − y app |M i,ρq ≤ C −1 H 7d/2+3 /q 2 . Since |iρq | =
O(ρdq ) and since M i,ρq is of condition number bounded by that of M , for each i the
−1 7d/2+3 /q 2 ) using
linear system of Line 3 can
 be solved efficiently (to accuracy O(C H
1
O(ρq ) = O ln max(  , q) iterations of the CG method (reminded in Subsection 4.9)
 with
1
d
d
a cost of O(ρq ) per iteration, which results in a total cost of O N ρq ln max(  , q) . Lines
(q),loc

(q)

(q)

4 and 5 are naturally of complexity O(N ρdq ). Since Ai,j
= 0 if τi and τj are at a
q
(k),loc and
distance larger than 2H ρq the complexity of Line 6 is O(N ρ2d
q ). Note that A
B (k),loc are banded and of bandwidth O(N ρdk ). It follows that Line 8 is of complexity
O(|I (k) |ρdk ). According to Theorem 5.6 the linear system of Line 9 needs to be solved
up to accuracy |y − y app |B (k),loc ≤ kgkH −1 (Ω) /2. Since B (k),loc is of uniformly bounded
condition number this can be done using O(ln 1 ) iterations of the CG method with a
cost of O(|I (k) |ρdk ) per iteration (using O(|J (k) |) = O(|I (k) |)), which results in a total
cost of O(|I (k) |ρdk ln 1 ) for Line 9. Storing the fine mesh values of u(k),loc − u(k−1),loc
in Line 11 costs O(N ρdk ) (since for each node x on the fine mesh only O(ρdk ) localized
basis functions contribute to the value of u(k),loc − u(k−1),loc ). According to Theorem
5.6, for each i ∈ I (k−1) the linear system B (i,ρ) y = b of Line 12 needs to be solved
up to accuracy |y − y app |B (i,ρ) ≤ C −1 H −k+7d/2+4 /(k − 1)2 . Since the matrix B (i,ρ)
inherits the uniformly bounded condition number from B (k),loc this can be done using
49

O(ln 1 ) iterations of the CG method with a cost of O(H −d ρdk−1 ρdk ) = O(ρdk−1 ρdk ) per
iteration. This results in a total cost of O(|I (k) |ρdk−1 ρdk ln 1 ) for Line 12. We obtain,
using the sparsity structures of D(k−1,k),loc and R(k−1,k),loc that the complexity of Line
d
13 is O(|I (k−1) |ρdk−1 H −d ) = O(|I (k−1) |ρdk−1 ) and that of Line 14 is O(|I (k−1) |ρ2d
k−1 ρk ).
The complexity of lines 15 to 16 is summarized in the display of Algorithm 2 and a simple
consequence of the sparsity structure of R(k−1,k),loc . Line 18 is complexity O(|I (1) |ρd1 ln 1 )
(using CG as in Line 9). As in Line 11, storing the values of u(1),loc costs O(N ρd1 ). Finally,
obtaining uloc in Line 20 costs O(N q) (observe that q = O(ln N )).
(k),loc

(k),loc

Compute and store ψi
, χi
, A(k),loc , B (k),loc
loc
loc
and u s.t. ku − u kH01 (Ω) ≤ kgkH −1 (Ω)
First solve
Subsequence solves

 ≤ Hq
N ln3d 1
N lnd+1 1

 ≥ Hq
N ln3d N
N lnd N ln 1

Subsequent solves to compute u(k),loc s.t.
ku − u(k),loc kH01 (Ω) ≤ CkgkL2 (Ω)

N lnd+1

1


−d lnd+1

1


−d lnd+1

1


(k)

Subsequent solves to compute the coefficients ci
P
(k) (k)
of u(k),hom = i∈I (k) ci ψi
(k),loc
s.t. ku − u
kH01 (Ω) ≤ C(kgkL2 (Ω) + kgkLip )
Subsequent solves to compute u(k),hom s.t.
ku − u(k),hom kL2 (Ω) ≤ C(kgkL2 (Ω) + kgkLip )
a periodic/ergodic with mixing length H p ≤ ,
first solve of u(k),hom s.t.
ku − u(k),hom kL2 (Ω) ≤ C(kgkL2 (Ω) + kgkLip )

(N (ln3d N )H p
+−d )
lnd+1 1

Table 1: Complexity of Algorithm 2.

Total computational complexity, first solve. Summarizing we obtain that the
(k),loc
(k),loc
complexity of Algorithm 2, i.e. the cost of computing the gamblets (ψi
), (χj
),
(k),loc
(k),loc
loc
their stiffness matrices (A
,B
), and the approximation u such that ku −
1 3d 
loc
u kH01 (Ω) ≤ kgkH −1 (Ω) is O N ln max( 1 , N d )
(with Line 14 being the correspond(1),loc

ing complexity bottleneck). The complexity of storing the gamblets (ψi
1 d 
and their stiffness matrices (A(1),loc , B (k),loc ) is O N ln max( 1 , N d ) .

(k),loc

), (χj

)

Computational complexity of subsequent solves with g ∈ H −1 (Ω). If (5.3) (i.e.
(k),loc
(k),loc
(1.1)) needs to be solved for more than one g then the gamblets ψi
, χi
and the
stiffness matrices B (k),loc do not need to be recomputed. The cost of subsequent solves

1 d
is therefore that of Line 9 i.e. O N ln max( 1 , N d ) ln 1 to achieve the approximation
accuracy ku − uloc kH01 (Ω) ≤ kgkH −1 (Ω) .
50

Computational complexity of subsequent solves with g ∈ L2 (Ω) and  ≥ H q .
If g ∈ L2 (Ω) (i.e. if kgkL2 (Ω) is used to express the accuracy of the approximation)
and  ∈ [H k , H k−1 ] then, by Theorem 5.5, u(k),loc achieves the approximation accuracy
ku − u(k),loc kH01 (Ω) ≤ CkgkL2 (Ω) (i.e. u(j+1),loc − u(j),loc does not need to be computed

for j ≥ k) and the corresponding complexity is O N (ln 1 )d+1 (if g ∈ H −1 (Ω) then the
energy of the solution can be in the fine scales and u(j+1),loc − u(j),loc do need to be
computed for j ≥ k).
Computational complexity of subsequent solves with g Lipschitz continuous
and  ≥ H q . Note that the computational complexity bottleneck for computing the
(k),loc
coefficients of u(k),loc in the basis (ψi
) when g ∈ L2 (Ω) and  ∈ [H k , H k−1 ] is in
the computation of the vectors g (j),loc for j > k. If g is Lipschitz continuous then
(k),loc
(k)
(k)
(k)
gi
be approximated g(xi ) where xi is any point in τi without loss of accuracy.
Note that this approximation requires (only) O(H −kd ) evaluations of g and leads to a
corresponding u(k),loc satisfying ku − u(k),loc ka ≤ C(kgkL2 (Ω) + kgkLip ) (with kgkLip =
supx,y∈Ω |g(x) − g(y)|/|x − y|). Therefore the computational complexity of subsequent
P
(k)
(k) (k),loc
solves to obtain the
coefficients ci in the decomposition u(k),loc = i∈I (k) ci ψi

is O −d (ln 1 )d+1 (i.e. independent from N if g is Lipschitz continuous). Of course,
obtaining an H01 (Ω)-norm approximation of u with accuracy H k requires expressing the
(k),loc
values of ψi
(and therefore u(k),loc ) on the fine mesh, which leads to a total cost of
1 d
O(N (ln  ) ). However if one is only interested expressing the values of u(k),loc on the fine
mesh in a sub-domain of diameter  then the resulting complexity is O((N d +−d )(ln 1 )d )
Computational complexity of subsequent L2 -approximations with g Lipschitz
(k)
(k)
continuous and  ≥ H q . Let (xi )i∈I (k) be points of (τi )i∈I (k) forming a regular
(k)
coarse mesh of Ω of resolution H k and write ϕi the corresponding (regular and coarse)
piecewise linear nodal basis elements. If (as in classical homogenization or HMM) one is
only interested in an L2 -norm approximation of u with accuracy H k then the coefficients
(k)
P
c
(k)
(k)
ci defined above are sufficient to obtain the approximation uhom = i∈I (k) R i (k) ϕi
φ
Ω i
R
R
(k)
(k)
that satisfies ku(k),loc − uhom kL2 (Ω) ≤ CH k ku(k),loc kH01 (Ω) ( Ω uhom φi = Ω u(k),loc φi )
and therefore ku − u(k),hom kL2 (Ω) ≤ C(kgkL2 (Ω) + kgkLip ). Note that the computational

complexity of subsequent solves to obtain uhom is O −d (ln 1 )d+1 .
Total computational complexity if a is periodic or ergodic with mixing length
H p and  ≈ H k with k ≥ p. Under the assumptions of classical homogenization
or HMM [33] (e.g. a is of period H p or a is ergodic with H p as mixing length), if
(k)
the sets τi are chosen to match the period of a and the domain is rescaled so that
1/H is an integer, then the entries of A(k) are invariant under periodic translations (or
stationary if the medium is ergodic). Therefore, under these assumptions, as in classical
homogenization, it sufficient to limit the computation of gamblets to periodicity cells (or

51

ergodicity cells with a tight control on mixing as in [47]). The resulting cost of obtaining
u(k),hom (in a first solve) such that ku(k),loc − u(k),hom kL2 (Ω) ≤ C(kgkL2 (Ω) + kgkLip ), is

O N ln3d N H p + −d ) lnd+1 1 .
Acknowledgements. The author gratefully acknowledges this work supported by the
Air Force Office of Scientific Research and the DARPA EQUiPS Program under awards
number FA9550-12-1-0389 (Scientific Computation of Optimal Statistical Estimators)
and number FA9550-16-1-0054 (Computational Information Games) and the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research,
through the Exascale Co-Design Center for Materials in Extreme Environments (ExMatEx, LANL Contract No DE-AC52-06NA25396, Caltech Subcontract Number 273448).
The author also thanks C. Scovel, L. Zhang, P. B. Bochev, P. S. Vassilevski, J.-L.
Cambier, B. Suter, G. Pavliotis and an anonymous referee for valuable comments and
suggestions.

References
[1] A. Abdulle and C. Schwab. Heterogeneous multiscale FEM for diffusion problems
on rough surfaces. Multiscale Model. Simul., 3(1):195–220 (electronic), 2004/05.
[2] G. Allaire and R. Brizzi. A multiscale finite element method for numerical homogenization. Multiscale Model. Simul., 4(3):790–812 (electronic), 2005.
[3] A. Averbuch, G. Beylkin, R. Coifman, and M. Israeli. Multiscale inversion of
elliptic operators. In Signal and image representation in combined spaces, volume 7
of Wavelet Anal. Appl., pages 341–359. Academic Press, San Diego, CA, 1998.
[4] I. Babuška, G. Caloz, and J. E. Osborn. Special finite element methods for a class
of second order elliptic problems with rough coefficients. SIAM J. Numer. Anal.,
31(4):945–981, 1994.
[5] I. Babuška and R. Lipton. Optimal local approximation spaces for generalized
finite element methods with application to multiscale problems. Multiscale Model.
Simul., 9:373–406, 2011.
[6] I. Babuška and J. E. Osborn. Generalized finite element methods: their performance and their relation to mixed methods. SIAM J. Numer. Anal., 20(3):510–536,
1983.
[7] I. Babuška and J. E. Osborn. Can a finite element method perform arbitrarily
badly? Math. Comp., 69(230):443–462, 2000.
[8] G. Bal and W. Jing. Corrector theory for MSFEM and HMM in random media.
Multiscale Model. Simul., 9(4):1549–1587, 2011.

52

[9] R. E. Bank, T. F. Dupont, and H. Yserentant. The hierarchical basis multigrid
method. Numer. Math., 52(4):427–458, 1988.
[10] M. Bebendorf. Efficient inversion of the Galerkin matrix of general second-order
elliptic operators with nonsmooth coefficients. Math. Comp., 74(251):1179–1199
(electronic), 2005.
[11] M. Bebendorf. Hierarchical matrices, volume 63 of Lecture Notes in Computational
Science and Engineering. Springer-Verlag, Berlin, 2008. A means to efficiently solve
elliptic boundary value problems.
[12] G. Ben Arous and H. Owhadi. Multiscale homogenization with bounded ratios
and anomalous slow diffusion. Comm. Pure Appl. Math., 56(1):80–113, 2003.
[13] A. Bensoussan, J. L. Lions, and G. Papanicolaou. Asymptotic analysis for periodic
structure. North Holland, Amsterdam, 1978.
[14] L. Berlyand and H. Owhadi. Flux norm approach to finite dimensional homogenization approximations with non-separated scales and high contrast. Archives for
Rational Mechanics and Analysis, 198(2):677–721, 2010.
[15] G. Beylkin and N. Coult. A multiresolution strategy for reduction of elliptic PDEs
and eigenvalue problems. Appl. Comput. Harmon. Anal., 5(2):129–155, 1998.
[16] A. Brandt. Multi-level adaptive technique (mlat) for fast numerical solutions to
boundary value problems. 1973. Lecture Notes in Physics 18.
[17] L. V. Branets, S. S. Ghai, L. L., and X.-H. Wu. Challenges and technologies in
reservoir modeling. Commun. Comput. Phys., 6(1):1–23, 2009.
[18] M. E. Brewster and G. Beylkin. A multiresolution strategy for numerical homogenization. Appl. Comput. Harmon. Anal., 2(4):327–349, 1995.
[19] F.-X. Briol, C. J. Oates, M. Girolami, M. A. Osborne, and D. Sejdinovic. Probabilistic integration: A role for statisticians in numerical analysis?
arXiv:1512.00933, 2015.
[20] E. J. Candès, J. Romberg, and T. Tao. Robust uncertainty principles: exact
signal reconstruction from highly incomplete frequency information. IEEE Trans.
Inform. Theory, 52(2):489–509, 2006.
[21] E. J. Candes and T. Tao. Near-optimal signal recovery from random projections:
universal encoding strategies? IEEE Trans. Inform. Theory, 52(12):5406–5425,
2006.
[22] V. Chandrasekaran, S. Sanghavi, P. A. Parrilo, and A. S. Willsky. Rank-sparsity
incoherence for matrix decomposition. SIAM J. Optim., 21(2):572–596, 2011.

53

[23] O. A. Chkrebtii, D. A. Campbell, B. Calderhead, and M. A. Girolami. Bayesian
solution uncertainty quantification for differential equations. Bayesian Analysis.
arXiv:1306.2365, 2016.
[24] E. Chow and P. S. Vassilevski. Multilevel block factorizations in generalized hierarchical bases. Numer. Linear Algebra Appl., 10(1-2):105–127, 2003. Dedicated to
the 60th birthday of Raytcho Lazarov.
[25] P. R. Conrad, M. Girolami, S. Srkk, A. Stuart, and K. Zygalakis. Probability measures for numerical solutions of differential equations. Statistics and Computing,
arXiv:1512.00933, pages 1–18, 2016.
[26] S. Demko, W. F. Moss, and P. W. Smith. Decay rates for inverses of band matrices.
Math. Comp., 43(168):491–499, 1984.
[27] P. Diaconis. Bayesian numerical analysis. In Statistical decision theory and related
topics, IV, Vol. 1 (West Lafayette, Ind., 1986), pages 163–175. Springer, New
York, 1988.
[28] D. L. Donoho. Compressed sensing. IEEE Trans. Inform. Theory, 52(4):1289–
1306, 2006.
[29] M. Dorobantu and B. Engquist. Wavelet-based numerical homogenization. SIAM
J. Numer. Anal., 35(2):540–559 (electronic), 1998.
[30] J. Duchon. Interpolation des fonctions de deux variables suivant le principe de la
flexion des plaques minces. Rev. Francaise Automat. Informat. Recherche Operationnelle Ser. RAIRO Analyse Numerique, 10(R-3):5–12, 1976.
[31] J. Duchon. Splines minimizing rotation-invariant semi-norms in Sobolev spaces. In
Constructive theory of functions of several variables (Proc. Conf., Math. Res. Inst.,
Oberwolfach, 1976), pages 85–100. Lecture Notes in Math., Vol. 571. Springer,
Berlin, 1977.
[32] J. Duchon. Sur l’erreur d’interpolation des fonctions de plusieurs variables par les
Dm -splines. RAIRO Anal. Numér., 12(4):325–334, vi, 1978.
[33] W. E and B. Engquist. The heterogeneous multiscale methods. Commun. Math.
Sci., 1(1):87–132, 2003.
[34] Y. Efendiev, J. Galvis, and P. S. Vassilevski. Spectral element agglomerate algebraic multigrid methods for elliptic problems with high-contrast coefficients. In
Domain decomposition methods in science and engineering XIX, volume 78 of Lect.
Notes Comput. Sci. Eng., pages 407–414. Springer, Heidelberg, 2011.
[35] Y. Efendiev, V. Ginting, T. Hou, and R. Ewing. Accurate multiscale finite element
methods for two-phase flow simulations. J. Comput. Phys., 220(1):155–174, 2006.

54

[36] Y. Efendiev and T. Hou. Multiscale finite element methods for porous media flows
and their applications. Appl. Numer. Math., 57(5-7):577–596, 2007.
[37] B. Engquist and E. Luo. Convergence of a multigrid method for elliptic equations
with highly oscillatory coefficients. SIAM J. Numer. Anal., 34(6):2254–2273, 1997.
[38] B. Engquist and O. Runborg. Wavelet-based numerical homogenization with applications. In Multiscale and multiresolution methods, volume 20 of Lect. Notes
Comput. Sci. Eng., pages 97–148. Springer, Berlin, 2002.
[39] B. Engquist and P. E. Souganidis. Asymptotic and numerical homogenization.
Acta Numerica, 17:147–190, 2008.
[40] G. E. Fasshauer. Meshfree methods. In Handbook of Theoretical and Computational
Nanotechnology. American Scientific Publishers, 2005.
[41] R. P. Fedorenko. A relaxation method of solution of elliptic difference equations.
Ž. Vyčisl. Mat. i Mat. Fiz., 1:922–927, 1961.
[42] A. Frieze, R. Kannan, and S.’ Vempala. Fast monte-carlo algorithms for finding
low-rank approximations. In In Proceedings of the 39th annual IEEE symposium
on foundations of computer science, pages 370–378, 1998.
[43] A. C. Gilbert, S. Guha, P. Indyk, S. Muthukrishnan, and M. Strauss. Near-optimal
sparse Fourier representations via sampling. In Proceedings of the Thirty-Fourth
Annual ACM Symposium on Theory of Computing, pages 152–161. ACM, New
York, 2002.
[44] A. C. Gilbert, M. J. Strauss, J. A. Tropp, and R. Vershynin. One sketch for all: fast
algorithms for compressed sensing. In STOC’07—Proceedings of the 39th Annual
ACM Symposium on Theory of Computing, pages 237–246. ACM, New York, 2007.
[45] E. De Giorgi. Sulla convergenza di alcune successioni di integrali del tipo dell’aera.
Rendi Conti di Mat., 8:277–294, 1975.
[46] A. Gloria. Analytical framework for the numerical homogenization of elliptic monotone operators and quasiconvex energies. SIAM MMS, 5(3):996–1043, 2006.
[47] A. Gloria, S. Neukamm, and F. Otto. Quantification of ergodicity in stochastic
homogenization: optimal bounds via spectral gap on Glauber dynamics. Invent.
Math., 199(2):455–515, 2015.
[48] L. Grasedyck, I. Greff, and S. Sauter. The al basis for the solution of elliptic
problems in heterogeneous media. Multiscale Modeling & Simulation, 10(1):245–
258, 2012.
[49] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. J. Comput.
Phys., 73(2):325–348, 1987.
55

[50] W. Hackbusch. A fast iterative method for solving Poisson’s equation in a general region. In Numerical treatment of differential equations (Proc. Conf., Math.
Forschungsinst., Oberwolfach, 1976), pages 51–62. Lecture Notes in Math., Vol.
631. Springer, Berlin, 1978.
[51] W. Hackbusch. Multigrid methods and applications, volume 4 of Springer Series
in Computational Mathematics. Springer-Verlag, Berlin, 1985.
[52] W. Hackbusch, L. Grasedyck, and S. Börm. An introduction to hierarchical matrices. In Proceedings of EQUADIFF, 10 (Prague, 2001), volume 127, pages 229–241,
2002.
[53] N. Halko, P. G. Martinsson, and J. A. Tropp. Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev., 53(2):217–288, 2011.
[54] R. L. Harder and R. N. Desmarais. Interpolation using surface splines. J. Aircraft,
9:189–191, 1972.
[55] P. Hennig. Probabilistic interpretation of linear solvers. SIAM Journal on Optimization, 25(1):234–260, 2015.
[56] P. Hennig, M. A. Osborne, and M. Girolami. Probabilistic numerics and uncertainty in computations. Proc. A., 471(2179):20150142, 17, 2015.
[57] M. R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear
systems. J. Research Nat. Bur. Standards, 49:409–436 (1953), 1952.
[58] T. H. Hou and P. Liu. Optimal local multi-scale basis functions for linear elliptic
equations with rough coefficient. Discrete and Continuous Dynamical Systems,
36(8):4451–4476, 2016.
[59] T. Y. Hou and X. H. Wu. A multiscale finite element method for elliptic problems
in composite materials and porous media. J. Comput. Phys., 134(1):169–189, 1997.
[60] G. S. Kimeldorf and G. Wahba. A correspondence between Bayesian estimation on
stochastic processes and smoothing by splines. Ann. Math. Statist., 41:495–502,
1970.
[61] S. M. Kozlov.
The averaging of random operators.
109(151)(2):188–202, 327, 1979.

Mat. Sb. (N.S.),

[62] F. M. Larkin. Gaussian measure in Hilbert space and applications in numerical
analysis. Rocky Mountain J. Math., 2(3):379–421, 1972.
[63] D. Le and D. S. Parker. Using randomization to make recursive matrix algorithms
practical. Journal of Functional Programming, 9:605–624, 11 1999.

56

[64] A. Malqvist and D. Peterseim. Localization of elliptic multiscale problems. Mathematics of Computation, 2014.
[65] J. Mandel, M. Brezina, and P. Vaněk. Energy optimization of algebraic multigrid
bases. Computing, 62(3):205–228, 1999.
[66] P.-G. Martinsson, V. Rokhlin, and M. Tygert. A randomized algorithm for the
decomposition of matrices. Appl. Comput. Harmon. Anal., 30(1):47–68, 2011.
[67] F. Murat. Compacité par compensation. Ann. Scuola Norm. Sup. Pisa Cl. Sci.
(4), 5(3):489–507, 1978.
[68] F. Murat and L. Tartar. H-convergence. Séminaire d’Analyse Fonctionnelle et
Numérique de l’Université d’Alger, 1978.
[69] D. E. Myers. Kriging, co-Kriging, radial basis functions and the role of positive
definiteness. Comput. Math. Appl., 24(12):139–148, 1992. Advances in the theory
and applications of radial basis functions.
[70] J. Nash. Non-cooperative games. Ann. of Math. (2), 54:286–295, 1951.
[71] A. S. Nemirovsky. Information-based complexity of linear operator equations. J.
Complexity, 8(2):153–175, 1992.
[72] J. Nolen, G. Papanicolaou, and O. Pironneau. A framework for adaptive multiscale
methods for elliptic problems. Multiscale Model. Simul., 7(1):171–196, 2008.
[73] E. Novak and H. Woźniakowski. Tractability of multivariate problems. Volume II:
Standard information for functionals, volume 12 of EMS Tracts in Mathematics.
European Mathematical Society (EMS), Zürich, 2010.
[74] A. O’Hagan. Bayes-Hermite quadrature. J. Statist. Plann. Inference, 29(3):245–
260, 1991.
[75] A. O’Hagan. Some Bayesian numerical analysis. In Bayesian statistics, 4
(Peñı́scola, 1991), pages 345–363. Oxford Univ. Press, New York, 1992.
[76] H. Owhadi. Anomalous slow diffusion from perpetual homogenization. Ann.
Probab., 31(4):1935–1969, 2003.
[77] H. Owhadi. Bayesian numerical homogenization.
13(3):812–828, 2015.

Multiscale Model. Simul.,

[78] H. Owhadi and E. J. Candès. Private communication. February 2004.
[79] H. Owhadi and C. Scovel. Brittleness of Bayesian inference and new Selberg formulas. Communications in Mathematical Sciences, 14:83–145, 2016. arXiv:1304.7046.
[80] H. Owhadi and C. Scovel. Toward Machine Wald, pages 1–35. Springer International Publishing, 2016. arXiv:1508.02449.
57

[81] H. Owhadi, C. Scovel, and T. J. Sullivan. Brittleness of Bayesian Inference under
finite information in a continuous world. Electronic Journal of Statistics, 9:1–79,
2015. arXiv:1304.6772.
[82] H. Owhadi, C. Scovel, and T. J. Sullivan. On the Brittleness of Bayesian Inference.
SIAM Review, 57(4):566–582, 2015.
[83] H. Owhadi, C. Scovel, T. J. Sullivan, M. McKerns, and M. Ortiz. Optimal uncertainty quantification. SIAM Rev., 55(2):271–345, 2013.
[84] H. Owhadi and L. Zhang. Metric-based upscaling. Comm. Pure Appl. Math.,
60(5):675–723, 2007.
[85] H. Owhadi and L. Zhang. Localized bases for finite dimensional homogenization
approximations with non-separated scales and high-contrast. SIAM Multiscale
Modeling & Simulation, 9:1373–1398, 2011. arXiv:1011.0986.
[86] H. Owhadi and L. Zhang. Gamblets for opening the complexity-bottleneck of
implicit schemes for hyperbolic and parabolic odes/pdes with rough coefficients.
arXiv:1606.07686, 2016.
[87] H. Owhadi, L. Zhang, and L. Berlyand. Polyharmonic homogenization, rough
polyharmonic splines and sparse super-localization. ESAIM Math. Model. Numer.
Anal., 48(2):517–552, 2014.
[88] E. W. Packel. The algorithm designer versus nature: a game-theoretic approach
to information-based complexity. J. Complexity, 3(3):244–257, 1987.
[89] I. Palasti and A. Renyi. On interpolation theory and the theory of games. MTA
Mat. Kat. Int. Kozl, 1:529–540, 1956.
[90] G. C. Papanicolaou and S. R. S. Varadhan. Boundary value problems with rapidly
oscillating random coefficients. In Random fields, Vol. I, II (Esztergom, 1979),
volume 27 of Colloq. Math. Soc. János Bolyai, pages 835–873. North-Holland,
Amsterdam, 1981.
[91] L. E. Payne and H. F. Weinberger. An optimal Poincaré inequality for convex
domains. Arch. Rational Mech. Anal., 5:286–292 (1960), 1960.
[92] H. Poincaré. Calcul des probabilités. Georges Carrés, Paris, 1896.
[93] K. Ritter. Average-case analysis of numerical problems, volume 1733 of Lecture
Notes in Mathematics. Springer-Verlag, Berlin, 2000.
[94] J. W. Ruge and K. Stüben. Algebraic multigrid. In Multigrid methods, volume 3
of Frontiers Appl. Math., pages 73–130. SIAM, Philadelphia, PA, 1987.
[95] A. Sard. Linear approximation. American Mathematical Society, Providence, R.I.,
1963.
58

[96] M. Schober, D. K. Duvenaud, and P. Hennig. Probabilistic ODE solvers with
Runge-Kutta means. In Z. Ghahramani, M. Welling, C. Cortes, N.D. Lawrence,
and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems
27, pages 739–747. Curran Associates, Inc., 2014.
[97] J. E. H. Shaw. A quasirandom approach to integration in Bayesian statistics. Ann.
Statist., 16(2):895–914, 1988.
[98] J. R Shewchuk. An introduction to the conjugate gradient method without the
agonizing pain. Technical report, Pittsburgh, PA, USA, 1994.
[99] J. Skilling. Maximum Entropy and Bayesian Methods: Seattle, 1991, chapter
Bayesian Solution of Ordinary Differential Equations. Springer Netherlands, Dordrecht, 1992.
[100] S. Spagnolo. Sulla convergenza di soluzioni di equazioni paraboliche ed ellittiche.
Ann. Scuola Norm. Sup. Pisa (3) 22 (1968), 571-597; errata, ibid. (3), 22:673,
1968.
[101] K. Stüben. A review of algebraic multigrid. J. Comput. Appl. Math., 128(1-2):281–
309, 2001. Numerical analysis 2000, Vol. VII, Partial differential equations.
[102] A. V. Sul0 din. Wiener measure and its applications to approximation methods. I.
Izv. Vysš. Učebn. Zaved. Matematika, 1959(6 (13)):145–158, 1959.
[103] W. Symes. Transfer of approximation and numerical homogenization of hyperbolic
boundary value problems with a continuum of scales. TR12-20 Rice Tech Report,
2012.
[104] J. F. Traub, G. W. Wasilkowski, and H. Woźniakowski. Information-based complexity. Computer Science and Scientific Computing. Academic Press, Inc., Boston,
MA, 1988. With contributions by A. G. Werschulz and T. Boult.
[105] J. A. Tropp. Recovery of short, complex linear combinations via l1 minimization.
IEEE Trans. Inform. Theory, 51(4):1568–1570, 2005.
[106] P. S. Vassilevski. Multilevel preconditioning matrices and multigrid V -cycle methods. In Robust multi-grid methods (Kiel, 1988), volume 23 of Notes Numer. Fluid
Mech., pages 200–208. Vieweg, Braunschweig, 1989.
[107] P. S. Vassilevski. On two ways of stabilizing the hierarchical basis multilevel
methods. SIAM Rev., 39(1):18–53, 1997.
[108] P. S. Vassilevski. General constrained energy minimization interpolation mappings
for AMG. SIAM J. Sci. Comput., 32(1):1–13, 2010.
[109] P. S. Vassilevski and J. Wang. Stabilizing the hierarchical basis by approximate
wavelets. I. Theory. Numer. Linear Algebra Appl., 4(2):103–126, 1997.
59

[110] P. S. Vassilevski and J. Wang. Stabilizing the hierarchical basis by approximate wavelets. II. Implementation and numerical results. SIAM J. Sci. Comput.,
20(2):490–514 (electronic), 1998.
[111] J. Von Neumann. Zur Theorie der Gesellschaftsspiele. Math. Ann., 100(1):295–320,
1928.
[112] J. von Neumann and O. Morgenstern. Theory of Games and Economic Behavior.
Princeton University Press, Princeton, New Jersey, 1944.
[113] A. Wald. Statistical decision functions which minimize the maximum risk. Ann.
of Math. (2), 46:265–280, 1945.
[114] W. L. Wan, Tony F. Chan, and Barry Smith. An energy-minimizing interpolation
for robust multigrid methods. SIAM J. Sci. Comput., 21(4):1632–1649, 1999/00.
[115] X. Wang. Transfer-of-approximation approaches for subgrid modeling. PhD thesis,
Rice University, 2012.
[116] H. Wendland. Scattered data approximation, volume 17 of Cambridge Monographs
on Applied and Computational Mathematics. Cambridge University Press, Cambridge, 2005.
[117] C. D. White and R. N. Horne. Computing absolute transmissibility in the presence
of finescale heterogeneity. SPE Symposium on Reservoir Simulation, page 16011,
1987.
[118] H. Woźniakowski. Probabilistic setting of information-based complexity. J. Complexity, 2(3):255–269, 1986.
[119] H. Woźniakowski. What is information-based complexity? In Essays on the
complexity of continuous problems, pages 89–95. Eur. Math. Soc., Zürich, 2009.
[120] Z. Min Wu and R. Schaback. Local error estimates for radial basis function interpolation of scattered data. IMA J. Numer. Anal., 13(1):13–27, 1993.
[121] J. Xu and Y. Zhu. Uniform convergent multigrid methods for elliptic problems with
strongly discontinuous coefficients. Math. Models Methods Appl. Sci., 18(1):77–105,
2008.
[122] J. Xu and L. Zikatanov. On an energy minimizing basis for algebraic multigrid
methods. Comput. Vis. Sci., 7(3-4):121–127, 2004.
[123] I. Yavneh. Why multigrid methods are so efficient. Computing in Science and
Engg., 8(6):12–22, November 2006.
[124] H. Yserentant. On the multilevel splitting of finite element spaces. Numer. Math.,
49(4):379–412, 1986.

60

