Multi-level Compressed Sensing Petrov-Galerkin discretization of high-dimensional parametric PDEs

0.63 MB·

by  Jean-Luc Bouchot

MULTI-LEVEL COMPRESSED SENSING PETROV-GALERKIN DISCRETIZATION
OF HIGH-DIMENSIONAL PARAMETRIC PDES
JEAN-LUCBOUCHOT,HOLGERRAUHUT,ANDCHRISTOPHSCHWAB
7
Abstract. Weanalyzeanovelmulti-levelversionofarecentlyintroducedcompressedsensing(CS)Petrov-
1
Galerkin(PG)methodfrom[H.RauhutandCh. Schwab: CompressivesensingPetrov-Galerkinapproxima-
0 tionofhigh-dimensionalparametricoperatorequations,Math. Comp. 304(2017)661–700]forthesolution
2
of many-parametricpartialdifferential equations. We proposetousemulti-levelPGdiscretizations, based

n onahierarchyofnestedfinitedimensionalsubspaces,andtoreconstructparametricsolutionsateachlevel
a from level-dependent random samples of the high-dimensional parameter space via CS methods such as
J weighted ℓ1-minimization. For affine parametric, linear operator equations, we prove that our approach
allows to approximate the parametric solution with (almost) optimal convergence order as specified by
6
certain summability properties of the coefficient sequence in a general polynomial chaos expansion of the

parametric solution and by the convergence order of the PG discretization in the physical variables. The
]
A computationsoftheparametersamplesofthePDEsolutionis“embarrasinglyparallel”,asinMonte-Carlo
Methods. Contrary to other recent approaches, and as already noted in [A. Doostan and H. Owhadi: A
N
non-adaptedsparseapproximationofPDEswithstochasticinputs. JCP230(2011)3015-3034]theoptimal-
. ityofthecomputedapproximationsdoesnotrequirea-prioriassumptionsonorderingandstructureofthe
h
index sets of the largest gpc coefficients (such as the “downward closed” property). We prove that under
t
a certainassumptionsworkversusaccuracyofthenewalgorithmsisasymptoticallyequaltothatofonePG
m solveforthecorrespondingnominalproblemonthefinestdiscretizationleveluptoaconstant.
[

1. Introduction
1
v Motivated in particular by uncertainty quantification, the numerical solution of parametric operator
1
equations has gained significant attention in recent years. In many cases, the underlying parameter space
7
is high dimensionaloreveninfinite dimensionalso that standardapproximationmethods are subjectto the
6
1 curse ofdimensionality, see e.g. [17, 16]. Monte Carlo (MC) sampling, however,may be usedin the context
0 thattheparametricmodelarisesfromastochasticmodelandleadstoamean-squarerateofm 1/2 interms
−
. of the number m of sample evaluations, with constants that are independent of the parameter dimension.
1
0 The (dimension-independent) rate 1/2 is not improvable in MC methods, in general, and the challenge
7 consistsindevelopingmethods thatachieveafasterconvergencerateandatthe sametime alleviateoreven
1 overcome the curse of dimensionality.
:
v A number of computationalapproacheshave emergedin recentyearstowardsthis end. Among these are
i adaptivestochasticGalerkinmethods,asdevelopedin[26,25,34],reducedbasisapproaches(see,eg.,[6,12]),
X
adaptive Smolyak discretizations [49, 50], adaptive interpolation methods [14] as well as sampling methods
r
a [52]. AdaptiveGalerkinmethods[26,25,34]areintrusiveinthesensethattheycannotsimplyreuseasolver
developed for the corresponding problem with fixed parameter. In contrast, the other above mentioned
methods and algorithms are non-intrusive, but they rely on successive numerical solutions of the operator
equations for various parameter instances that are chosen based on suitable precomputations. In contrast,
(multilevel)Monte-Carlo(MLMC)[44],orQuasi-MonteCarloapproaches(QMC)[23]compute expectations
or statistical moments of the (random parametric) solution via solutions for parameter instances chosen at
random or “quasi-random”,which allows to compute the “parameter snapshot” solutions in parallel.
Date:January9,2017.
CSsupportedinpartbyEuropeanResearchCouncil(ERC)throughthegrantAdG247277,JLBandHRsupportedinpart
bytheERCgrantStG258926. JLBandHRwouldliketothanktheHausdorffResearchInstituteforMathematics,University
ofBonn,wherepartsofthisworkhavebeenperformedduringthetrimesterprogramMathematics ofSignalProcessing. JLB
thanksthesupportoftheClayMathematicsInstituteforhisvisittotheCRMtoattendtheIRPConstructiveApproximation
andHarmonicAnalysiswherepartsofthisworkhasbeendone.
1
2 JEAN-LUCBOUCHOT,HOLGERRAUHUT,ANDCHRISTOPHSCHWAB
In this article, we build on a compressed sensing approach for numerically computing parametric solu-
tions developed and analyzed in [46, 9] (see also [24, 45] for earlier work, and [15] for recent developments)
and combine it with ideas originating from MLMC methods, see e.g. [33, 4]. For Petrov-Galerkin (PG)
discretizations on a finite hierarchy of nested subspaces, ordered with respect to discretization levels, the
presently proposed method “samples”, in a judicious fashion, the parameter space and computes corre-
sponding PG approximations for random choices of the parameter vector. As in MLMC-PG approaches,
the number of such snapshot evaluations decreases with increasing discretization level (corresponding to
increasing refinement of the discretization). In contrast to (ML)MC sampling, we employ a CS technique
based on weighted ℓ -minimization [47, 1] or Iterative Hard Thresholding Pursuit [40, 7] in order to re-
1
construct the coefficients of a generalized polynomial chaos expansion of the difference of the parametric
solution at two subsequent discretization levels. Finally, these differences are summed together to obtain
a PG approximation of the parametric solution at the finest level. One contribution of this paper is to
show that the generalized polynomial chaos (GPC) expansion of the differences of PG approximations of
the parametric solution is approximately sparse by estimating the weighted ℓ -norm for 0 < p < 1 of the
p
sequence of Chebyshev coefficients by a term that depends in a controlled way on the discretization level.
This fact makes the presently developed, multi-level version of the compressive sensing approach feasible.
We provide dimension-independent convergence rates which exceed 1/2 under certain sparsity assumptions
on the parametric solution family of the operator equation and estimate the computational complexity for
achievingsuchrates. SimilartoMLMCmethods,theworkloadforapproximatingtheparametricsolutionis
asymptotically the same as the one for computing one snapshot solution at the finest levelup to a constant
thatdependsonlyonsmoothnessparametersandp (0,1). However,incontrasttomultilevelMonteCarlo,
∈
the convergencerates affordedby our scheme arepractically independent ofthe dimensionandonly limited
by the solutions’ sparsity; in particular, they may significantly exceed (m 1/2).
−
O
In mathematical terms, we consider linear, parametric operator equations of the generic form
A(y)u(y)=f. (1)
Heretheparametervectory U liesinahigh-dimensionalspaceU makingitchallengingtocomputationally
∈
approximatethe solutionmapy u(y), due tothe mentionedcurse of dimensionality, anotiongoingback
7→
toR.E.Bellman[5],see[17,16]foritsrelevanceinthepresentcontext. Assumingthattheparametervector
y = (y )d takes values in finite intervals, we can consider, without loss of generality, U =[ 1,1]d, where
j j=1 −
the parameter set dimension d may be finite or infinite.
In our setting, the parametric family of operators A(y) : maps from a reflexive Banach space
′
X → Y
to the topological dual of, potentially, another reflexive Banach space . A canonical example is the
X Y
affine-parametric diffusion equation considered in [18, 19] and in the single-level version of the present
work [46, 9]. For a bounded Lipschitz domain D Rn (one should think of n = 1,2,3) and a parametric
⊂
diffusion coefficient a(,y) L (D) that depends affinely on a parameter vector y, i.e.,
∞
· ∈
a(x,y)=a¯(x)+ y ψ (x), x D, (2)
j j
∈
j 1
X≥
we consider the model parametric, second order divergence form elliptic Dirichlet problem
A(y)u:= (a(,y) u)=f in D, u =0. (3)
∂D
−∇· · ∇ |
The weak formulation of (3) in the Sobolev space = := H1(D) reads: Given f , for every
y U :=[ 1,1]N find u(y) such that X Y 0 ∈ Y′
∈ − ∈X
a(x,y) u(x) v(x)dx = f(x)v(x)dx, for all v . (4)
∇ ·∇ ∈Y
ZD ZD
Eq. (3) is a particular example of an affine-parametric operator equation of the form
A(y):=A + y A , y =(y ) U :=[ 1;1]N, (5)
0 j j j j 1
≥ ∈ −
j 1
X≥
with A := (ψ ),A = (a¯ ). In (5), the operator A ( , ) is traditionally referred
j j 0 0 ′
−∇· ∇ −∇· ∇ ∈ L X Y
to as nominal operator or mean field while the operators A ( , ), for j 1, are referred to as
j ′
∈ L X Y ≥
fluctuations. For the parametric problem to be well-posed uniformly with respect to the parameter y U,
∈
MULTI-LEVEL COMPRESSED SENSING DISCRETIZATION OF HIGH-DIMENSIONAL PARAMETRIC PDES 3
awnedaaspsupmliceatbhilaittyofj≥ou1rkAapjpkrLo(Xac,Yh′w) <ill∞beisnpewchifiaetdfoalhloeawds.. PFaurrathmeertraiscsuexmppatniosinosnsresquucihreadsf(o5r)tchaencboenovbertgaeinnecde
P
e.g.byaKarhunen-Lo`eveexpansionofrandominputdatafordivergence-formpartialdifferentialequations,
as explained in [51, 19].
Inordertoensurewell-posednessoftheparametricdiffusionproblem(3)asin[19]werequiretheuniform
ellipticity assumption: there exist constants 0<r R< such that
≤ ∞
r a(x,y) R, for almost all x D, for all y U . (6)
≤ ≤ ∈ ∈
The Lax-Milgram Lemma ensures that for every y U, the weak formulation (4) admits a unique solution
∈
u(,y) which satisfies the uniform a priori estimate
· ∈X
sup u(y) r−1 f ′.
y Uk kX ≤ k kY
∈
Here and throughout the remainder, the term “uniform” refers to uniform with respect to the parameter
sequence y U.
∈
For the sake of simplicity we detail here only the approximationof functionals of solutions to the
′
G ∈X
parametric operator equation (1), i.e., we are interested in the numerical approximation of
F(y):= (u(y)), y U =[ 1,1]d,
G ∈ −
pointwise with respect to y. We expect that our approachcan be generalized to the recoveryof the vector-
valued solution map y u(y), but we postpone this generalization to later contributions. We are aiming
7→
at numerical schemes that are:
Reliable: the convergence and accuracy should be verified and customizable;
•
Parallelizable: parallel sampling as in Monte-Carlo methods should be allowed, with a convergence
•
rate in terms of the number of samples which (up to possibly logarithmic terms) equal the best
possiblerateensuredbythecompressibilityofF(y),i.e.,byweightedℓ -estimatesoftheChebyshev
p
coefficients of F;
Non-intrusive: the approximation should use existing numerical solvers of the problem with fixed
•
parameters, without any re-implementation of PDE solvers.
It is important to notice the difference to usual MC methods where the results obtained from random
samplingusuallyholdinexpectation. Incontrast,ourapproachprovidesapproximationsthatholdpointwise
with respect to y. We estimate the coefficients of a tensorized Chebyshev expansion; whence only matrix-
vectormultiplicationsarerequiredinordertocomputethesolutionF(y)= (u(y))foranygivenparameter
G
vector y = (y )d up to a prescribed accuracy. The computation scheme analyzed here differs from the
j j=1
single-levelone introduced in [46] in the sense that computing the approximationis done in a more efficient
and computationally tractable manner. To this end, an unknown function u(y) is approximated by a
telescopic sequence of so-called “details” at successively finer spatial resolutions: u(y) L (ul(y)
≈ l=1 −
ul 1(y)) where ul corresponds to a PG approximation on a discretization level l. This is analogous to
−
P
MLMC methods, but is achieved here by compressive sensing of the parameters with a level-dependent
number of parameter samples y(i) on each discretization level in the physical domain.
We outline key ideas of the compressive sensing approach. We assume at our disposal a countable
orthonormalbasis (ϕ ) of L2(U,η) with η denoting a probability measure on the parameter set U to be
ν ν Λ
specified,anddenoteby∈L2(U,η; )theBochnerspaceofstronglymeasurablemapsfromU tothe(separable
X
Hilbert) space containing solution instances, which are square integrable w.r.t. η. We represent any
X
function u(y) with values in as u(y)= α ϕ (y), where α= (α ) denotes the unique sequence
X ν Λ ν ν ν j Λ
ofcoefficients α . Hence, inorderto comp∈ute anapproximationofthe∈parametricsolutionfor anyy it
ν
∈X P
sufficestocalculateanapproximationofthecoefficientsα . Foranewinputparametery,oneevaluatesthe
ν
basis functions ϕ at y and forms a linear combination to recovera direct estimation of the solution. Later
ν
on, we analyze the use of tensorized Chebyshev polynomials as orthonormal basis. The approximation is
computedbyevaluatingthefunctionatafewparameterpointsy(i),1 i m,andsolvingthelinearsystem
g =Φα, where g =(g )m = u(y(i)) m and where Φ corresponds t≤o th≤e sensing matrix Φ Rm N with
i i=1 i=1 ∈ ×
entries Φ =ϕ (y(i)), where N correspondsto the number ofbasis functions taken for the approximation.
i,ν ν (cid:0) (cid:1)
However at this stage the coefficients α and the components g are elements in , and therefore, we first
ν i
X
4 JEAN-LUCBOUCHOT,HOLGERRAUHUT,ANDCHRISTOPHSCHWAB
deal with the simpler case where a functional (also known as the Quantity of Interest (QoI for short) in
G
the uncertainty quantification literature) is applied to the solution, resulting in
b = (g )= u(y(i)) = z ϕ (y(i)), i=1,...,m, z = (α ), ν Λ.
i i ν ν ν ν
G G G ∈
(cid:16) (cid:17) νX∈Λ
We are particularly interested in the situation that the number m of evaluations is smaller than the cardi-
nality of Λ, so that the linear system b = Φz is underdetermined. Approximate sparsity of the coefficient
sequence (α ), and of (z ), allows to apply techniques from compressive sensing such as (weighted) ℓ -
ν ν 1
minimization or iterative hard thresholding (pursuit) in order to recover z accurately. In fact, approximate
sparsity follows from the fact that ( α ) and (z ) are contained in weighted ℓ ( )-spaces, as shown in
ν ν p
k kX F
[3, 16, 18, 19, 46] and, for a related coefficient sequence, in this paper.
We expect that anapproximationof the full solutionu(y), y U,taking values in the function space ,
∈ X
canbecomputedbyavariantofourcompressivesensingscheme. Onemayuseideasfromjoint/blocksparsity
inordertorecoverthesequence(α )withα viamixedℓ /ℓ -minimization,seee.g.[27,28,30](atleast
ν ν 1 2
∈X
inthecasethat and areHilbertspaces). However,wepostponeadetailedanalysistoalatercontribution
X Y
and restrict ourselves here to the simpler case of recovering the real-valued function y (u(y)).
7→G
The multi-level approximation scheme uses discretization levels l = 1,...,L, where the meshwidth at
discretization level l is 2 lh , so that the finest discretization is h = 2 Lh . With n being the dimension
− 0 L − 0
of the domain D, we assume that the number of degrees of freedom at level l scales like (2nl), and we
O
further assume available linear complexity, multilevel solvers for the approximate solution of the discretized
linear system of equations (uniformly with respect to the parameter y) resulting in computational costs
per PG solution ul(y(i)) that scales linearly in the number of degrees of freedom: (2nl).
O
The presently proposed multi-level extension of the CS PG approach from [46] proceeds analogous to
MLMC(see,e.g.,[39]or[33]andthe referencestherein): forparameterchoices y(i) ondiscretiza-
{ l }i=1,...,ml
tion level l, compute PG solutions ul(y(i)), ul 1(y(i)) at two consecutive discretization levels l and l 1
l − l −
(setting u0 0). From the differences dul(y(i)) = ul(y(i)) ul 1(y(i)), we compute an approximation
≡ l l − − l
dul(y) via the single level compressive sensing approach of [46] for each l = 1,...,L. Finally, we combine
the approximations at all levels similarly as in MLMC methods, i.e., u(y) = L dul(y), to obtain an
l=1
apeproximation of the full parametric solution. The main result of this paper consists of an analysis of this
P
method and provides, in its proof, a strategy on how to choose the numbeer ml of parameeter points at each
level l. Its precise statement, Theorem 9, is postponed to later in the exposition. To illustrate the type of
resultsobtainedhere,westatenowaversionofTheorem9intheparticularcaseofalinear,divergenceform
diffusion operator with affine dependence on the parameters (see Eqs. (2) and (3)). Ahead, we say that the
weight sequence v is constant, if it is of the form v = β for j = 1,...,d for some β > 1 and v = for
j j
∞
j > d, which corresponds to the case that the expansion (2) is finite (with d terms). We say that v has
polynomial growth if v = cjα, j N, for some c > 1, α > 0. We refer to Section 5.3 for details on the
j
∈
weight sequences.
Theorem 1. Let L N and γ (0,1). Consider the diffusion equation (3) with affine parametric co-
efficient (2), forcing ∈term f H∈1+t(D) and functional H 1+t′(D), with the respective smoothness
− −
∈ G ∈
parameters t,t 0. Assume that in (2) holds a¯ Wt, (D) and that the fluctuations fulfill the weighted
′ ∞
≥ ∈
p-summability1
kψjkpWt,∞(D)vj2−p <∞,
j 1
X≥
for a sequence v = (v ) of weights as well as the following stronger, weighted version of the Uniform
j j 1
≥
Ellipticity Assumption (6): there exists 0<r R< such that
≤ ∞
vj(2−p)/p|ψj(x)|≤min{a¯(x)−r,R−a¯(x)}, for all x∈D. (7)
j 1
X≥
1To ease the presentation, here and throughout the paper, we have not highlighted the dependence of the summability
parameterpontheregularitytoftheright-hand-sidef. Itshouldbenotedthatthecompressibilityofthegpcexpansion,the
choiceoftheweightsequence, thenumberofsamplesperlevelalldependontheregularityofthedataa¯,ψj,D andf.
MULTI-LEVEL COMPRESSED SENSING DISCRETIZATION OF HIGH-DIMENSIONAL PARAMETRIC PDES 5
With probability at least 1 γ, the function F(y):= (u(y)), y U, can be approximated by L (weigthed)
− G ∈
compressed sensing approximations based on a sequence of Galerkin projections into spaces of piecewise
polynomials on regular, simplicial triangulations of meshwidth h =2 ℓh from
ℓ − 0
m &max s log3(s )log(N ),log(L/γ)
l l l l
{ }
solution evaluations at discretization level l for l = 1,...,L, where s 2(L l)(t+t′)p/(1 p), and N is the
l − − l
≍
size of the (level-dependent) active set of tensorized Chebyshev polynomials. The resulting approximation
F# satisfies
kF −F#k∞ ≤CpkfkH−1+t(D)kGkH−1+t′(D)L2−(t+t′)Lht0+t′
kF −F#k2 ≤Cp′kfkH−1+t(D)kGkH−1+t′(D)2−(t+t′)Lht0+t′
Under the assumption that the computational cost of a single solve at level l scales linearly with respect to
the number of degrees of freedom, i.e., is (2nl) (for an n-dimensional domain D), this result is achieved
O
witha totalwork for thecomputation of snapshot solutions that scales as max 2nL,Lβ2L(t+t′)p/(1 p) ,
−
O
where β =4 for constant weights v and β =5 for polynomially growing wei(cid:16)ghts. Tnhe constant hidden inoth(cid:17)e
-notation includes a factor of log(d) in the case of constant weights.
O
We note in passing that in what follows, the estimates of the overallcomputational work do not account
forthe numericalsolutionofthe convexprogramsrequiredfor the compressedsensingapproximationofthe
mapping F. We justify this convention by the observation that the computational cost of ℓ -minimization
1
is often of lower order compared to the total cost of evaluating the PDE samples.
Our theorem shows that in the case of sufficiently strong summability, i.e., (t+t′)p < n, at a total
1 p
cost that scales as a constant times a single PDE solve at the finest discretization le−vel L, the multilevel
CSPG (MLCSPG) strategy can approximate a fixed function F for any parameter vector y U. This is
∈
analogousto what is affordedby MLMC methods, but the presentMLCSPG strategyallowsto achieveany
convergencerateaffordedbythegpcsummability,andallowstoapproximatethefullparametricdependence,
while MLMC only yields expectations (or moments). Moreover, in our case the computational work scales
favorablywithdecreasingp,whichcorrespondstobettersparseapproximationratesimpliedbytheweighted
p-summability of (norms of) polynomial chaos coefficients of the parametric solution. To be more precise,
in the case of higher smoothness t+t >0, we obtain an approximation error that scales with ht+t′. With
′ L
a small enough value of p, we may exploit smoothness in the physical domain (allowing t+t such that
′
(t+t′)p < n) and balance approximation error for the PDE solves. In contrast, the computational work
1 p
req−uired by MLMC to achieve an expected approximation error scaling as ht grows proportionally to 22tL
L
when2t n(wheretcorrespondstothesmoothnessofthesolutioninthephysicaldomain),see[4,Theorem
≥
5.7], and there is no parameter p in MLMC whose tuning allows to avoid such growth.
Nevertheless,wenotethattandpmaynotbetunedindependently: inmanyinstancesincreasedsmooth-
ness t leads to a larger value of the summability parameter p.
We emphasize that the tools and results developed here do not require a particular structure on the
support set of the best approximation. It is often the case (see e.g. [14, 43]), that proofs and/or methods
require the sets of active indices in N-term gpc approximations be downward closed, their approximation
properties then being, in particular, independent of the polynomial system adopted for implementation. In
constrast, the presently proposed, compressed sensing based approach can recover (with high probability)
any support set of active multi-degrees of tensorized Chebyshev polynomial approximations (only assuming
veryroughknowledgeof its locationas providedby weightedℓ -estimates of polynomialchaoscoefficients),
p
yetstill providingquasi-optimalratesofconvergence. Moreover,apartfromthe ℓ -minimizationpartofthe
1
algorithm, all function evaluations can be done in parallel.
Theorem 1 is a particular case of our main Theorem 9 which we prove in Section 4 after recalling some
basics about Petrov-Galerkinapproximations in Section 2 combined with compressed sensing techniques in
Section 3. Section 5 deals with pratical aspects such as truncating the dimension of the parameter space.
The paper is finally concluded by numerical experiments to illustrate the theory in Section 6.
6 JEAN-LUCBOUCHOT,HOLGERRAUHUT,ANDCHRISTOPHSCHWAB
2. Petrov-Galerkin approximations
We deal with the pointwise numerical approximations of the countably-parametric operator equation
Eq. (1). Numerically accessing the parametric solution map U y u(y) at a fixed parameter instance
∋ 7→
y U requires discretization of Eq. (1) also in “physical space”. To this end, we introduce two dense,
∈
one-parameter families of discretization spaces h and h of equal finite dimensions
h>0 h>0
{X } ⊂ X {Y } ⊂ Y
Nh :=dim( h)=dim( h) and assume that the parametric operator A(y) fulfills the discrete and uniform
X Y
inf sup conditions: there exists a µ>0 such that for any h>0 and y U
− ∈
 iinnff006=6=wvhh∈∈XYhhssuupp006=6=wvhh∈∈XYhh kkhhvvAAhh((kkyyXX))vvkkhhww,,wwhhkkhhYYii ≥≥µµ>>00. (8)
The PG projections are defined as the solution to the following weak variational problems:

Find uh(y):=Gh(y)(u(y)), such that A(y)uh(y),vh = f,vh for all vh h . (9)
h i h i ∈Y
We recall the following classical result (see for example [8, Chapter 6]).
Proposition 1. Let h and h be discretization spaces for the PG method, such that the uniform discrete
X Y
inf sup conditions (8) are fulfilled and assume that the bilinear operator (u,w) A(y)u,w is
− X ×Y ∋ 7→h i
continuous, uniformly with respect to y U.
∈
Then the PG projections Gh(y) : h are well-defined linear operators, whose norms are uniformly
X → X
bounded with respect to the parameters y and h, i.e.,
1
supsup uh(y) f ′, (10)
y Uh>0k kX ≤ µk kY
∈
C
supsup Gh(y) (11)
y Uh>0k kL(X) ≤ µ
∈
The Galerkin projections are uniformly quasi-optimal: for every y U we have the a-priori error bound
∈
C
u(y) uh(y) 1+ inf u(y) vh . (12)
k − kX ≤ µ vh∈Xhk − kX
(cid:18) (cid:19)
As is classical in the theory of polynomial approximation (see, e.g. [20, 48]), we use a holomorphic
extension of the parametric operator family A(y) to complex parameter sequences z U, where is
∈O ⊃ O
some suitable subset of the complex plane. Here, when dealing with extensions of operators and solutions
to parameters taking values in the complex domain, we identify the function spaces and with their
X Y
complexifications 1,i and 1,i for the sake of simplicity. We require the parametric operator
X ⊗{ } Y ⊗{ }
z A(z)to be holormorphicwithrespecttoanyfinite setofvariablesandto be boundedly invertible.
O∋ 7→
Hereby,aBanach-spacevaluedmappingz R(z) E ofasinglecomplexvariableissaidtobeholomorphic
7→ ∈
(in some open domain ) if
O
R(z +h) R(z )
0 0
lim −
h 0 h
exists in E for any z , with h 0 un→derstood in C. Note that our assumption on A requires it to
0
∈ O →
be holomorphic with respect to any component z of z independently. Joint holomorphy with respect to
j
an arbitrary, finite subset of variables z = (z ) with Λ < of z then follows from Hartogs’ theorem.
′ j j Λ
In the sequel, we will often assume that the op∈en set | ,|on∞which z A(z) is holomorphic, contains
the product of Bernstein ellipses = with O:= (z +z 1)7→/2,z C : z = σ . In the case
Eρ j 1Eρj Eσ { − ∈ | | }
of complex-parametric operators, the boun≥ded invertibility of A(z) is equivalent to the complex discrete
N
inf sup conditions: there exists a constant µC >0 such that for any h>0 and z
− ∈O
 iinnff006=6=vwhh∈∈XYhhssuupp006=6=wvhh∈∈XYhh RReekkhhvvAAhh((kkzzXX))vvkkhhww,,wwhhkkhhYYii ≥≥µµCC >>00,. (13)
Approximationresultsondiscretizationspacesareusuallycombinedwithpriorknowledgeoftheregularity

ofthedata. Forthis,weassumethatthereexistsa0<t tsuchthattheparametricfamilyA(y) ( , )
′
≤ ∈L X Y
MULTI-LEVEL COMPRESSED SENSING DISCRETIZATION OF HIGH-DIMENSIONAL PARAMETRIC PDES 7
is regular in given smoothness scales , resp. , satisfying:
t t 0 t t 0
{X } ≥ {Y } ≥
= ... , = ... . (14)
0 1 t 0 1 t
X X ⊃X ⊃ ⊃X Y Y ⊃Y ⊃ ⊃Y
Here, the smoothness index t 0 denotes, for example, a differentiation order in a scale of Sobolev or
≥
Besov spaces. These spaces are defined by interpolation for non integer indices. We shall also require a
corresponding scale on the dual side, with :=( ) , and :=( ) :
Xt′ X′ t Yt′ Y′ t
= ... , = ... . (15)
X′ X0′ ⊃X1′ ⊃ ⊃Xt′ Y′ Y0′ ⊃Y1′ ⊃ ⊃Yt′
Note carefully that with this notation, ( ) generally differs from . For example, in the case of the
Xt ′ Xt′
diffusion equation (3), one may choose = H1(D) and = H1+t(D). In this case, = H 1(D) =
X 0 Xt 0 X′ −
(H01(D))′ and (Xt)′ =H−1−t(D)6=H−1+t(D)=(X′)t =:Xt′.
A first statement of solutionregularityin the scales (14), (15) takes the formof uniform bounded invert-
ibility of the family of parametric operators A(y):
A(y)∈L(Xt,Yt′), for all y ∈U, and ysupUkA(y)−1kL(Yt′,Xt) <∞. (16)
∈
ForthePGdiscretization,weassumeathandtwoone-parameterfamilies h and h of and
h>0 h>0
{X } {Y } X
of , respectively, with finite, equal dimension: N = dim( h) = dim( h) < . We assume furthermore
h
Y X Y ∞
that h and h are dense in and in , respectively. Here the discretization parameter
h>0 h>0
{X } {Y } X Y
h > 0 usually stands for the meshwidth in finite element discretizations of fixed polynomial degree, on a
quasiuniform triangulation of the physical bounded, polyhedral domain D. We assume that these spaces
admit approximation properties in the smoothness scales (14), (15),
inf w wh C ht w , for all w ,
infwvhh∈∈YXhhkkv−−vhkYkX≤≤Ct′htt′kkvkkYXt′t, for all v ∈∈YXt′t. (17)
Together with the bounded invertibility of the family of operators A(y), it holds:
Eq. (12) Eq. (17) Eq.(16)
ku(y)−uh(y)kX ≤Cvhinfhku(y)−vhkX ≤cthtku(y)kXt ≤CthtkfkYt′. (18)
∈X
Here, the constant C depends on a uniform bound on the inverse of the parametric operator in the ap-
t
pdrisocprreitaizteatsimonooptahrnaemssetseprahc.e: supy∈UkA(y)−1kL(Yt′,Xt), and on the smoothness parameter t, but not on the
Moreover,aswe confine the expositionto functionals ofsolutionsF(y)= (u(y)) forsome () ,we
′
G G · ∈X
assume adjoint regularity, i.e., there exists t 0, such that , and such that the parametric adjoint
′ ≥ G ∈ Xt′′
solution w (y) of the problem
′
G ∈Y
A(y)∗w (y)= (19)
G G
satisfies w (y) uniformly with respect to y:
G ∈Yt′′
sup w (y) ′ C ′ . (20)
y Uk G kY ≤ kGkXt′
∈
Under the adjointregularity(20), the uniformparametricdiscreteinf-supcondition(8)andthe approxima-
tion property (17), an Aubin-Nitsche duality argument as, e.g., in [41], implies superconvergence: for any
y U, with Fh the functional applied to the parametric PG solution uh(y) defined in (9) on discretization
∈
spaces of parameter h,
F(y) Fh(y) Ct+t′ht+t′ f ′ ′ . (21)
| − |≤ k kYtkGkXt′
3. Single-Level Compressed Sensing Petrov-Galerkin approximations
The multi-level compressed sensing PG (MLCSPG) discretization is a generalization of the single-level
algorithms and results developed in [46]. Analogous to MLMC path simulations (see e.g. [33] and the
referencesthere)orMLMCFiniteElementdiscretizations(seee.g.[4])theMLCSPGmethoddescribedhere
considers a sampling scheme from [46] with a number of sampling points depending on the discretization
level.
Such compressedsensing reconstructiontechniques have alreadyshownpromise in the contextof numer-
ical solutions of PDEs on high-dimensional parameter spaces: we refer, for example, to [55, 24, 45, 46, 9].
8 JEAN-LUCBOUCHOT,HOLGERRAUHUT,ANDCHRISTOPHSCHWAB
The key idea in these works is to decompose the solution of Eq. (1) via its (tensorized Chebyshev or Le-
gendre)polynomialchaosexpansionwithrespecttotheparametervectory. Astronglymeasurablemapping
u : U : y u(y) which is square (Bochner-) integrable with respect to the Chebyshev measure dη
→ X 7→
over U can be represented as a gpc expansion, i.e.,
u(y)= u T (y), (22)
ν ν
ν
X∈F
where in this case the coefficients in this expansion are functions u . Here := ν NN : supp(ν) <
ν ∈X F { ∈ 0 | |
is the set of multi-indices with finite support. The tensorized Chebyshev polynomials are defined as
∞}
∞
T (y)= T (y )= T (y ), y U, ν , (23)
ν νj j νj j ∈ ∈F
jY=1 j∈sYupp(ν)
with the univariate Chebyshev polynomials defined by
T (t)=√2cos(jarccos(t)), and T (t) 1. (24)
j 0
≡
Defining the probability measure σ on [ 1;1]as dσ(t):= dt , the univariate Chebyshev polynomials T
− π√1 t2 j
defined in (24) form an orthonormal system in L2([ 1,1];σ) i−n the sense that
−
1
T (t)T (t)dσ(t)=δ , k,l N .
k l k,l 0
∈
Z−1
Similarly, with the product measure
dy
dη(y):= dσ(y )= j ,
j
π 1 y2
Oj≥1 Oj≥1 − j
q
the tensorized Chebyshev polynomials (23) are orthonormalwith respect to η in the sense that
T (y)T (y)dη(y)=δ , µ,ν .
µ ν µ,ν
∈F
Z
y U
∈
Aresultprovenin[37]ensurestheℓ summability,forsome0<p 1,ofthepolynomialchaosexpansion(22)
p
≤
for the diffusion case, Eq. (3):
( u ) p = u p <
k νkX ν∈F p k νkX ∞
ν
(cid:13) (cid:13) X∈F
under the uniform ellipticity assump(cid:13)tion(6) andt(cid:13)he condition thatthe sequence ofinfinity normsofthe ψ
j
is itself ℓ summable:
p
p
( ψ ) = ψ p < .
k jk∞ j≥1 p k jk∞ ∞
(cid:13)(cid:13) (cid:13)(cid:13) Xj≥1
Recentresultsby[3]showthatthese(cid:13)conditionsca(cid:13)nbeimprovedbyconsideringpointwiseconvergenceofthe
series ψ insteadofinfinitynormsinthewholedomainD. Thistakesadvantageofthelocalstructure
j 1| j|
ofthebas≥iselementsψ ,e.g. whenonlyfewofthemareoverlapping,asisthecaseforwavelets. Inparticular,
j
P
ℓ summability of Legendre coefficients can be obtained when ( ψ ) ℓ for q := 2p/(2 p) provided
p k jk∞ j ∈ q −
that the interiors of the supports of the ψ do not overlap. The summability results from [37] concerning
j
Chebyshev expansions were extended to weighted ℓ estimates for the generalparametric operator problem
p
(1) with affine dependence as in (5), in [46] under slightly stronger assumptions. This result is particularly
importantforusasitensurestherecoveryofthe coefficientsu (oranyfunctionalthereof)viaCSmethods.
ν
The results on the approximation via an MLCSPG framework rely on the single-level results developed
in [47, 46], where functions are approximated via a weighted-sparse expansion in an appropriate basis. We
review here the main ideas. Given a (finite) orthonormal system (φ ) , with Λ = N < for L2(U,η)
where η is a probability measure, for any fixed function f : U νRν,∈Λthere ex|ist|s a uniqu∞e sequence of
→
coefficients f =(f ) such that
ν ν Λ
∈
f(y)= f φ (y), y U . (25)
ν ν
∀ ∈
ν Λ
X∈
MULTI-LEVEL COMPRESSED SENSING DISCRETIZATION OF HIGH-DIMENSIONAL PARAMETRIC PDES 9
We define an ℓ norm associated with this expansion as f := f .
p ||| |||p k kp
In particular, a function f is said to be sparse (or compressible) if its sequence of coefficients in expan-
sion (25) is sparse (or compressible) itself. Our goal is to recover this said sequence from seemingly few
evaluations of the function f at certain (here random) sample points y(i), for 1 i m. This can be done
≤ ≤
by CSmethods: after introducingthe sensingmatrix Φ as Φ :=φ (y(i)) andletting g :=f(y(i)), itholds
i,j j i
g =Φf .
Hence, assuming that the expansion is sparse,and the number of samples rather small, we are dealing with
the by-now classical problem of recovering a sparse vector from few linear measurements, by solving, for
instance, the convex program
min z , subject to Φz =g. (26)
1
z RNk k
∈
In our context, it is beneficial to use a weighted framework which has been developed recently in [47].
Introducing a sequence of positive weights (ω ) with ω 1 for all ν, a weighted ℓ (quasi-)norm
ν ν Λ | ν| ≥ p
(henceforth indexed ℓ when appropriate) can be∈defined as
ω,p
f p := ω2 p f p, 0<p 2.
k kω,p ν− | ν| ≤
ν Λ
X∈
In particular, it holds f = f and f = f ω , where defines the pointwise multiplication.
ω,2 2 ω,1 1
k k k k k k k ⊙ k ⊙
Moreover,choosing the constant weight ω =1 yields the original definitions of ℓ norms. Formally letting
ν p
p 0 motivates the introduction of the weighted sparsity measure
↓
f := ω2 .
k kω,0 ν
ν∈ΛX,fν6=0
A vector x is therefore called weighted s-sparse (with respect to a weight sequence ω) if x s. We
ω,0
k k ≤
may therefore define the error of best weighted s-term approximation as
σ (f)=σ (f):= inf f z .
ω,s ω,s ω,p
z:kzkω,0≤sk − k
Withtheseweightederrormeasuresathand,theBasisPursuitproblem(26)canbegeneralizedtoinclude
a-prioriinformation encoded in the sequences ω of weights, as
min z , subject to Φz =g. (27)
ω,1
z RNk k
∈
More details on such weighted spaces and weighted sparse approximations can be found in [47] where the
following fundamental result is also proved.
Theorem 2. Suppose (φ ) is a finite orthonormal system with Λ = N < and that weights ω
ν ν Λ | | ∞ ν ≥
φ are given. For a (wei∈ghted) sparsity s 2 ω 2 , draw
ν
k k∞ ≥ k k∞
m Cslog3(s)log(N) (28)
≥
sample points y(i) at random, 1 i m, according to the orthonomalization measure η. The constant
≤ ≤
C >0 in (28) is universal, i.e., independent of all other quantities including s, m and N.
Then, withprobabilityatleast1 N log(s)3,anyfunctionf = f φ canbeapproximatedbythefunction
− ν ν
−
f := f φ , where f is the solution to the weighted basis pursuit problem (27). The approximation holds
ν ν P
in the following sense:
P
b b b
f f f f c σ (f) , and f f d σ (f) /√s.
1 s ω,1 2 1 s ω,1
k − k∞ ≤ − ω,1 ≤ k − k ≤
(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)
In particular, usingbthe w(cid:12)e(cid:12)(cid:12)ightedb(cid:12)(cid:12)S(cid:12)techkin inequality from [47] b
(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)
σs(f)ω,q s ω 2 1/q−1/p f ω,p, p<q 2, ω 2 <s, (29)
≤ −k k∞ k k ≤ k k∞
we obtain the estimates, for p<(cid:0)1, in terms(cid:1)of the sparsity
f f cs1 1/p f , f f ds1/2 1/p f . (30)
− ω,p 2 − ω,p
k − k∞ ≤ k k k − k ≤ k k
b b
10 JEAN-LUCBOUCHOT,HOLGERRAUHUT,ANDCHRISTOPHSCHWAB
Choosing m sln3(s)log(N) relates the reconstruction error and the number m of samples as
≍
log3(m)log(N) 1/p−1 log3(m)log(N) 1/p−1/2
f f c f , f f d f .
ω,p 2 ω,p
k − k∞ ≤ m k k k − k ≤ m k k
(cid:18) (cid:19) (cid:18) (cid:19)
Remark 1.bRecent works [11, 38] on restricted isometry conbstants for subsampled Fourier matrices suggest
that the factor log3(s) in (28) can be reduced to log2(s), but details remain to be worked out.
4. Multi-level Compressed Sensing Petrov-Galerkin approximations
4.1. A multi-level framework. We extend the foregoing CS methods to sweeping the parameter domain
to multi-level (“ML” for short) discretizations of the parametric problems, in the spirit of the ML MC
methods for numerical treatment of operator equations with random inputs as developed in [39, 33, 4].
There, the solution of the parametric operator equation (1) is approximated on a sequence of partitions of
the physical domain D of widths h L for a prescribed, maximal refinement level L N. To simplify
{ l}l=1 ∈
the exposition, we assume dyadic refinement, i.e. h = h /2 = 2 lh for a given, small enough, initial
l+1 l − 0
resolution h >0.
0
For a given parameter sequence y, we may write the Galerkin projection uL(y) hL of u(y) as
∈X
L
uL(y)= ul(y) ul−1(y), (31)
−
l=1
X
where we define u0(y) 0 (note that we will equivalently parametrize the approximations and spaces by l
≡
or h ). The idea behind our MLCSPG approach is to estimate every difference between consecutive levels
l
of approximation (the details) via a single level CSPG as presented above. For the remaining, we let
dul(y):=ul(y) ul 1(y), 1 l L, (32)
−
− ≤ ≤
denote the difference between two scales of approximation.
Asalreadyoutlinedintheintroduction,ourmethodproducespointwisenumericalapproximationsdul(y)
of dul(y) via a (single level) CSPG method. For each level l, we choose a number m of parameter vectors
l
y(1),...,y(ml), compute the PG approximations ul(y(i)) and ul 1(y(i)) by solving the correspondingcfinite
l l l − l
dimensional linear systems, and form the samples dul(y(i)) = ul(y(i)) ul 1(y(i)), i = 1,...,m . From
l l − − l l
these samples, one approximates the coefficients in the tensorizedChebyshev expansionof dul via weighted
ℓ -minimization (or any sparse recovery method). This yields approximations dul(y), ℓ=1,...,L, and
1
L
uL (y):= dul(y) c
MLCS
l=1
X
then provides an approximation of the targeted parametricc solution u = u(y). The convergence of our
MLCSPG framework can be estimated via the triangle inequality,
L
u(y) uL (y) u(y) uL(y) + dul(y) dul(y) . (33)
MLCS
k − kX ≤k − kX −
Xl=1(cid:13) (cid:13)X
For simplicity, we constrain our considerations to a functional (cid:13)(cid:13) ′ appliedcto th(cid:13)(cid:13)e solution, leading to
G ∈ X
the real-valued QoI F(y) = (u(y)) to be approximated. The above considerations apply verbatim when
G
replacing u(y) by F(y), and dul by dFl, the levelwise PG approximation, and dul(y) by dFl(y). The
triangle inequality leads to the error estimate
c d
L
F(y) FL (y) F(y) FL(y) + dFl(y) dFl(y) . (34)
MLCS
− ≤| − | −
(cid:12) (cid:12) Xl=1(cid:12) (cid:12)
The first term on the(cid:12)right hand side of E(cid:12)q. (33) can be estimated(cid:12)(cid:12)using the udniform(cid:12)(cid:12)parametric regularity
(16), the uniform parametric inf-sup condition (8) and the approximation property (17): for a regularity
parameter 0<t t of the data f,
≤
ku(y)−uL(y)kX ≤CthtLkfkYt′. (35)