Journal of Scientific Computing manuscript No.
(will be inserted by the editor)

Low-Dimensional Spatial Embedding Method for
Shape Uncertainty Quantification in Acoustic
Scattering
Yuval Harness

arXiv:1704.07727v1 [cs.CE] 25 Apr 2017

Received: date / Accepted: date

Abstract This paper introduces a novel boundary integral approach of shape
uncertainty quantification for the Helmholtz scattering problem in the framework
of the so-called parametric method. The key idea is to form a low-dimensional
spatial embedding within the family of uncertain boundary deformations via the
Coarea formula. The embedding, essentially, encompasses any irregular behavior
of the boundary deformations and facilitates a low-dimensional integration rule
capturing the bulk variation of output functionals defined on the boundary. In a
second phase a matching parametric grid is imposed. For the ease of presentation
the theory is restricted to 2D star-shaped obstacles in low-dimensional setting.
We employ the null-field reconstruction technique which is capable of handling
large shape deformations. Higher spatial and parametric dimensional cases are
discussed, though, not extensively explored in the current study.
Keywords Uncertainty Quantification · Shape Uncertainty · Helmholtz ·
Parametric Method · Low-dimensional Embedding · Coarea Formula
Mathematics Subject Classification (2010) MSC 49Q15 · MSC 65C99 ·
MSC 65N35 · MSC 65R20 · MSC 65Z05

1 Introduction
Considerable effort has been devoted in recent years to develop robust and efficient
computational strategies for the simulation of physical phenomena, that take into
account shape uncertainty. Often, the problem is formulated as an elliptic partial
differential equation (PDE) whose domain boundaries are uncertain. Such problems arise due to imperfections in manufacturing processes, e.g., in nano-optics
where the production of nano particles is, often, inaccurate relatively to nanoscale electromagnetic wave lengths [2]. Other examples arise in the context of
Y. Harness
Inria Bordeaux
Tel.: +33-5-35-00-26-19
Fax: +33-5-61-19-30-00
E-mail: yuval.harness@inria.fr

2

Yuval Harness

inverse problems, such as tomography where the visual representation of some
hidden object is constructed by partial, and possibly noisy, measurements [25].
The common practice for quantifying uncertainty in computational models,
is employ to the parametric method. The theoretical basis for the method was
laid down by Wiener [32]. The method itself was initially developed by Ghanem
and Spanos [9], and later generalized by Xiu and Karniadakis [28, 35, 36]. In this
approach the uncertain parameters are replaced by random quantities, and the
problem is recast as a system with random input. The solution is estimated via
global expansion of the random variables into a basis of uncorrelated functions.
Thus, the stochastic problem is transformed into a deterministic system in higher
dimension. The most popular expansions employed are the Karhunen-Loève expansion, and the generalized Polynomial Chaos (gPC) expansion.
The parametric approach, often, demonstrates superior performance in terms of
computational effort over other traditional methods, see [34] for a detailed review.
However, when the physical domain of the problem is uncertain the quantification
by the parametric approach becomes much more challenging. The main difficulty
stems from the need to employ an accurate discretization which captures desired
features of the solution for any realization of the random shape. Producing a
family of physical discretization schemes even for a selected finite set of random
realizations, generally, requires a considerable effort. The latter is especially true
when the geometry undergoes complicated changes as a function of the random
parameters. Note, that a single deterministic computation for a specific realization,
typically, requires a large computational effort as well.
To overcome the difficulties associated with the quantification of random shape
or domain, various techniques have been proposed. Typically, these are classified as
one of the following: perturbation, fictitious domain, level-set or random domain
mapping. Perturbation techniques [29] are straightforward and simple to apply,
however, their applicability is limited to small shape deformations. The fictitious
domain [3] and level-set methods [22, 23] are based on embedding the random
domain in a larger, deterministic domain containing all possible realizations. These
methods are capable of handling very irregular non-smooth geometries. However,
the embedding introduces non-smoothness in the spatial region, that intersects
with the random boundary. Thus, high-order convergence is only partially ensured
in the entire computational domain.
The random domain mapping method [26, 27] is the most common tool used
for solving PDEs on uncertain domains. The method is based on a realizationdependent coordinate transformation uniformly mapping all the realizations of
the domain to a fixed, reference configuration. The variational formulation of the
PDE on the random domain can then be posed on the reference domain, reducing
the problem to a PDE on a fixed domain with stochastic coefficients. The transformed PDE whose domain is fixed is solved using standard techniques. However,
the method is highly sensitive to the non-linear dependence of the problem on
the random boundary. In case of complex evolution of the shape, the random
coefficients can exhibit highly varying behavior. Even worse, non-smoothness of
the boundaries can produce singularities which require non-trivial treatment. To
overcome these difficulties a very accurate discretization grid is typically imposed,
often combined with dimensionality reduction techniques, e.g., sparse grids, to
ensure reasonable computational effort. See [4, 11, 14] for further details.

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

3

In this work an alternative method that attempts to mitigate the difficulties
associated with the more standard techniques for PDEs on uncertain domains is
proposed. The method is of boundary integral type [5] and, thus, can handle large
shape deformations with great ease. Our analysis is based on the observation that
as a function of the random interface the solution is piece-wise smooth in the
spatial domain, i.e., R2 or R3 . The proposed method constructs a discretization
grid of the random surface for all possible realizations in two stages. In the first
stage a spatial low-dimensional embedding of the family of random surfaces is
constructed via the Coarea formula [8]. The embedding, essentially, captures any
irregular behavior of the random surface and a discretization is applied only on a
compact region in R2 or R3 . In the second stage a parametric grid corresponding to
the low-dimensional spatial grid is imposed. A sparse or hierarchical parametric
grid can be applied for dealing with high dimensionality, while the spatial grid
effectively ensures that the bulk variation of output functionals defined on the
boundary is captured. In general, the method allows the handling of non-trivial
geometries without the loss of accuracy in the region intersecting with the random
interface.
Since this is a first case study and for the ease of presentation, the discussion has been limited to time-harmonic wave scattering by star-shaped obstacles.
In the analysis and numerical study a 2D scattering object and low-dimensional
parametric space are assumed. More complicated examples in higher spatial and
parametric dimensions are discussed. However, in-depth study of this topic is deferred to future work. For its simplicity, acoustic fluid-structure interaction has
been chosen as the physical application. In that case, the solution represents small
oscillations of pressure in a compressible ideal fluid. The method and ideas presented in this work can also be applied to electrodynamics and elastodynamics.
This work employs the null-field approach [20, 30, 31, 33], which in contrary to
the better known boundary element method (BEM) [24] and the Nyström method
[16], does not involve singular integrals. Null-field methods are fast and much easier to implement compared to BEM and the Nyström method. Their applicability
range, however, is typically limited to relatively low wavenumbers. The null-field
reconstruction technique [1, 13] is inherently stable, admits a-priori error evaluation, and facilitates the extraction of features of interest without prior estimation
of the entire solution. The method enables us to perform analysis from a purely
geometric point of view, which avoids the additional complications associated with
integration of weakly singular kernels. Combining low-dimensional surface embedding with BEM and Nyström method can be foreseen in a future study.
The paper is organized as follows. The fundamentals of the null-field reconstruction method for the time-harmonic wave scattering problem is presented in
Section 2. Section 3 reviews the procedure of optimal reconstruction from a numerical linear algebra point of view. Section 4 consists of the main theoretical
results of this work and includes the formulation of the problem. In Section 5 the
proposed method is applied to the general class of randomly shaped polygonal
cylinders, as a proof of concept that the suggested method can, indeed, handle
complex non-smooth shapes. Summary of the results as well as suggestions for
applying the method in more complicated scenarios are given in Section 6.

4

Yuval Harness

2 Null-Field Reconstruction for Time Harmonic Wave Scattering
In this section a brief review on the null-field reconstruction method for the timeharmonic wave scattering problem is given. First the time-harmonic acoustic scattering problem is presented, after which follows a review of the fundamental theory
of null-field methods. In the final part the main idea of the null-field reconstruction
technique for time-harmonic wave scattering is presented.

2.1 Acoustic Scattering by Impenetrable Obstacles
Let B denote a bounded domain in Rd representing an impenetrable obstacle with
boundary S, where d = 2 or 3. We denote by B = B ∪ S the closure of B. Let
Rd \ Bd be the unbounded exterior region occupied by a uniform medium. Let
r ∈ R denote a general spatial point,

(r cos θ, r sin θ) ,
if d = 2 ,
r=
(r sin ϕ cos θ, r sin ϕ sin θ, r cos ϕ) , if d = 3 .
For an incident time-harmonic field uinc (r) ’illuminating’ the obstacle, the scattered field usca (r) satisfies the following exterior boundary value problem:
∆usca + κ2 usca = 0 ∀ r ∈ Rd \ B ,

(1)

usca = −uinc or ∂ν usca = −∂ν uinc on S ,

(2)

d−1
2

(∂r usca − iκusca ) → 0 as r → ∞ ,
(3)
√
where i = −1. Equation (1) is known as the Helmholtz equation, where ∆
is the Laplacian and κ is the wavenumber. Equation (2) specifies the boundary
condition, depending on the physical problem: adopting the acoustic terminology,
it is sound-soft for Dirichlet problems and sound-hard for Neumann problems.
Here ν is the unit outward normal to S and ∂ν u is the normal derivative of u.
The last condition (3), known as the Sommerfeld radiation condition, ensures that
the scattered field propagates from the obstacle to infinity. The solution of the
exterior scattering problem is unique. A solution to Helmholtz equation is called
a wavefunction. A wavefunction satisfying the Sommerfeld condition (3) is called
an outgoing wavefunction.
r

2.2 Null-Field Theory Fundamentals
Null-field methods for the acoustic scattering problem (1) are based on Green’s
second theorem

ZZ
Z 
∂v
∂w
(w∆v − v∆w) dV =
w
−v
dS ,
(4)
∂ν
∂ν
R
S
which holds for any bounded domain R with a Lipschitz piecewise smooth boundary S, where v and w are scalar fields, and ∂v/∂ν and ∂w/∂ν denote corresponding
normal derivatives.

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

5

Let ψ be an outgoing wavefunction and let u denote the total field, uinc + usca .
Assuming ψ is analytic in Rd \ B, it can be shown by (4) that


Z 
Z 
∂ψ
∂u
∂ψ
∂uinc
u
−ψ
dS =
uinc
−ψ
dS .
∂ν
∂ν
∂ν
∂ν
S
S
Thus, for a sound-soft obstacle (u = 0 on S)

Z
Z 
∂ψ
∂uinc
∂u
dS = −
−ψ
dS ,
ψ
uinc
∂ν
∂ν
∂ν
S
S

(5)

while for a sound-hard obstacle (∂ν u = 0 on S)

Z
Z 
∂ψ
∂ψ
∂uinc
u
dS =
uinc
−ψ
dS .
∂ν
∂ν
∂ν
S
S

(6)

Using (5) or (6) an infinite set of equations can be produced, from which u or ∂ν u
on S are approximated. In practice, one chooses a finite subset of equations of the
form of (5) and (6) which are employed to optimally reconstruct the scattered field
without an explicit estimation of u or ∂ν u on S. The core idea of reconstruction
by functionals is presented in the next subsection, while the numerical procedure
for its practical implementation is covered in Section (3).

2.3 Reconstruction of Surface Functionals
Typically we are interested in estimating features of interest which are expressed by
the unknown surface density, u or ∂ν u on S. Often such features are the outcomes
of functionals in an appropriate Hilbert space. Indeed, let h denote the complex
conjugate of u or ∂ν u on S. Then for a general non-smooth surface S, the surface
density h belongs to the complex Hilbert space L2dS (S) whose inner-product and
norm are defined by
Z
q
hg , hiL2 =
h∗ (r)g(r) dS , khkL2 = hh , hiL2 ,
dS

S

dS

dS

where h∗ denotes the complex conjugate of h and dS is the induced volume form
on the surface. Recall that by Riesz representation theorem any bounded linear
functional L2dS → C operating on surface densities, is of the form hf , hiL2 . We
dS
call such functionals surface functionals.
In this work we focus on the estimation of the scattering coefficients of the
expansion of usca to cylinder or spherical harmonics in the 2D or 3D cases, respectively. These coefficients, denoted by bm or bm
n , satisfy
(P
∞
bm Hm (κr)eimθ
if d = 2 ,
usca (r) = Pm=−∞
∞ Pn
m (1)
m
b
h
(κr)Y
(r̂)
if d = 3 ,
n
n
n=0
m=−n n
(1)

where Hm (z) denotes the mth-order Hankel function of the first kind, hn (z) is
the nth-order spherical Hankel function of the first kind, and
Ynm = Pnm (cos ϕ) eimθ ,

6

Yuval Harness

where Pnm are the associated Legendre functions. See [5] for further details. The
scattering coefficients are very useful features of the surface density, since they
can easily express other important quantities such as the far-field pattern and the
radar cross section [12].
Consider the sound-soft case (5). Using the Hilbert space notation, it follows
that the scattering coefficients for the 2D and 3D cases satisfy
E
iD
Jm (κr)e−imθ , ∂ν u∗ (r) 2 = bm ,
∀ m ∈ Z,
(7)
4
LdS
D
E
+
− jn (κr)Pn|m| (cos ϕ)eimθ , ∂ν u∗ (r) 2 = bm
(8)
n , ∀ 0 ≤ |m| ≤ n ∈ Z ,
LdS

respectively, where Jm (z) denotes the mth-order Bessel function of the first kind,
jn (z) is the nth-order spherical Bessel function of the first kind, and ∂ν u∗ (r)
denotes the complex conjugate of the surface density ∂ν u(r). In practice only
{bm }|m|≤µ satisfying
o
n
(9)
µ = d3κrmax e , ρmax = inf |r| r ∈ Rd \ B ,
are required for an accurate description of usca , see [13] for more details. Similar
expressions can be derived for the sound-hard case.
The core idea of the reconstruction procedure is to approximate the outcome
of target functionals (7) and (8) without producing an explicit approximation of
the surface density ∂ν u on S. This, generally, allows us to handle much more
accurately complex geometries as well as irregular or singular surface densities.
Explicitly, we approximate the outcomes (7) and (8) by linear combinations of the
following form
E
DX
E
iD
Jm (κr)e−imθ , ∂ν u 2 ≈
cm,l ψ`2D , ∂ν u 2 ,
4
LdS
LdS
D
E
E
DX
− jn (κr)Pn|m| (cos ϕ)eimθ , ∂ν u 2 ≈
cm,l ψ`3D , ∂ν u 2 ,
LdS

LdS

where {ψ`2D } and {ψ`3D } are predetermined sets of functionals, in 2D and 3D
respectively, whose outputs are either known or can be calculated directly. We call
such functionals the information functionals.
As shown in (5) the outcome of the information functionals are readily available
if ψ`2D and ψ`3D are outgoing wavefunctions whose singularities are located in B.
Hence, given the information
D
D
E
E
ψ`3D , ∂ν u 2 = a3D
ψ`2D , ∂ν u 2 = a2D
` ,
` ,
LdS

LdS

the outcome of the target functional can be approximate by
E
X
iD
Jm (κr)e−imθ , ∂ν u 2 ≈
cm,l a2D
` , ∀ m ∈ Z,
4
LdS
and
D
E
− jn (κr)Pn|m| (cos ϕ)eimθ , ∂ν u

L2dS

≈

X

cm,l a3D
` ,

∀ 0 ≤ |m| ≤ n ∈ Z+ .

Obtaining the coefficients {cm,` } while ensuring measurable error bounds of the
estimations of the scattering coefficients can be achieved by reconstruction kernel
approximation which is the main topic of Section (3).

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

7

3 Optimal Reconstruction with A-priori Error Estimate
In this section the procedure of recovery of target functionals by information functionals in general Hilbert space setting is reviewed. Reconstruction problems often
involve regularization parameters which govern the stability and accuracy of the
procedure. The result of the optimization of the reconstruction with respect to the
regularization parameters is referred to as optimal reconstruction. Optimal reconstruction can be traced back to the notion of optimal recovery [10, 21]. A more
modern analysis from an inverse problem point of view can be found in [17, 18, 19].
We begin with the definition of the reconstruction problem and the notion
of reconstruction kernel. This is followed by a brief description of the numerical
procedure including error analysis. The final part elaborates on proper numerical
integration rules, that are needed for the error estimates. The method and error
analysis presented here, as well as further technical details have been initially
introduced in [13].

3.1 The Reconstruction Problem and Reconstruction Kernels
Let H be a complex Hilbert space, whose inner-product is denoted by h , iH . The
reconstruction problem is to approximate a finite set of target functionals
hgm , hiH = bm ,

gm ∈ H ,

m = 1, 2, . . ., M ,

(10)

l = 1, 2, . . ., L ,

(11)

by a given finite set of information functionals,
hf` , hiH = a` ,

f` ∈ H ,

where the element h ∈ H is unknown.
L
denote
Definition 1 Let k kH denote the norm induced by h , iH in H, and let Cm
P
L
L
a closed convex subset of C . A linear combination `=1 cbm,` f` whose coefficients
L
satisfy the minimality condition
(b
cm,1 , . . ., cbm,L ) ⊂ Cm

gm −

L
X

≤ gm −

cbm,` f`

`=1

L
X
`=1

H

L
∀ (cm,1 , . . ., cm,L ) ∈ Cm
,

cm,` f`

(12)

H

is called an optimal reconstruction kernel of the target functional gm by the inforL
mation functionals {f` } over Cm
.
Remark 1 Clearly, (12) is a projection on a convex set. The key point which is
L
. Note that almost no prior
addressed later, is how to determine the convex set Cm
knowledge on the element h is assumed.
We will show in the next subsection, that obtaining (12) vastly exceeds our
needs. In practice, it is sufficient to obtain an approximation satisfying
gm −

L
X
`=1

≤  kgm kH

cbm,` f`

L
(b
cm,1 , . . ., cbm,L ) ∈ Cm
,

(13)

H

with respect
P to some predetermined threshold,  > 0. In that case the linear combination L
bm,` f` is simply called a reconstruction kernel (i.e., not optimal).
`=1 c

8

Yuval Harness

3.2 The Discrete Reconstruction Procedure with Error Analysis
For the evaluation of the reconstruction kernel, we assume a finite dimensional
discretization satisfying the following definition.
Definition 2 Let HB be a bounded subset of a Hilbert space H whose innerproduct is denoted by h , iH . A mapping,
THB : HB → C1×P ,

P ∈ N,

is called an inner-product preserving discretization of HB of accuracy dis > 0 if
−
→−
f→
g ∗ − hf, giH < dis ,

∀ f, g ∈ HB ,

(14)

−
→
−
→
→
→
where f = THB (f ) and −
g = THB (g). The vectors f and −
g are called the
corresponding inner-product preserving discretizations of f and g on HB .
−
→ →
Let f` , −
gm ∈ C1×P denote inner-product preserving discretizations of some
f` , g` ∈ HB , respectively. Let η ∈ H denote the orthogonal projection of h on the
subspace spanned by HB . By definition (2) we obtain
−
→−
−
→∗
η ∗ − hf` , hiH , −
g→
f` →
m η − hgm , hiH ≤ kηkdis ,
for all ` ∈ {1, 2, . . ., L} and m ∈ {1, 2, . . ., M }. Hence, given the information (11)
and an approximation of −
g→
m,
gbm =

L
X

−
→
cbm,` f` ,

(15)

`=1

we can reconstruct the unknown target coefficients (10) via
→
bm = hgm , hiH ≈ gbm −
η∗=

L
X

L

−
→→∗ X
cbm,` f` −
η ≈
cbm,` a` .

`=1

(16)

`=1

To evaluate the error of the reconstruction (16) we denote for each ` = 1, 2, . . . , L
and m = 1, 2, . . . , M the discretization errors
−
→→∗
−
→∗
` = hf` , hiH − f` −
η , δm = hgm , hiH − −
g→
m η ,
and obtain the following estimate
bm −

L
X

→
cbm,` a` = δm + (−
g→
bm )−
η∗+
m −g

`=1

L
X

cbm,` ` ,

`=1

→
where our assumption ensure that |δm |, |` | ≤ kηkdis . Note that (−
g→
bm )−
η ∗ is
m −g
the projection error which can not be reduced if the set of information functionals,
{f` }, is predetermined.
To control the error we impose the following regularization constraint
|b
cm,` | ≤ evl /dis ,

evl ≥ dis ,

(17)

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

9

where evl is a chosen or given evaluation error bound. Thus, we obtain
bm −

L
X

→
g→
bk )−
η ∗ + L · kηk · evl .
cbm,` a` ≤ kηk · dis + (−
m −g

(18)

`=1
L
The regularization constraint (17) explicitly defines the convex set Cm
in (12) as


evl
L
= (c1 , c2 , . . . , cL ) |c` | ≤
Cm
.
dis

The error estimate (18) implies that it is sufficient to obtain an approximation
(15) satisfying k(−
g→
bk )k ≤ evl . Indeed, in that case (18) reduces to
m −g
bm −

L
X

cbm,` a` ≤ kηk · dis + (L + 1) · kηk · evl .

`=1

P
Often, the summation of evaluation errors L
bm,` ` is not cumulative. Thus,
`=1 c
the overall error is typically O (evl ) assuming kηk = O(1). This facilitates an
efficient approximation technique which is based performing successive singular
value decompositions on subsets of information functionals. The technique was
presented in [13] and demonstrated high stability and good convergence properties.
Further details including different variants of the technique can be found in [12].
3.3 Inner-Product Preserving Discretization and Numerical Integration
Obtaining inner-product preserving discretizations of surface functionals is a fundamental issue. Let us focus on the case, as in this work, where {f` } and {gm } are
smooth functions in some H = L2 space with an inner-product,
Z
hf , giH =
f g ∗ ω ds ,
D

where D ⊂ R is compact and Jordan measurable, g ∗ is the complex conjugate of
g and ω is a proper weight function.
We compute the integrals by a numerical integration, typically by some standard rule of the form
d

hf, giH ≈

N
X

f (σ (i) )g ∗ (σ (i) )ω (i) ,

i=1
(1)

with integration nodes σ , . . ., σ (N ) contained in D and real positive weights
ω (1) , . . ., ω (N ) . Discretizing an element f ∈ H as a weighted gridfunction

√
√
−
→ 
f = f (σ (1) ) ω (1) , . . ., f (σ (M ) ) ω (N ) ,
essentially, satisfies the inner-product preserving assumption (14) if N ∈ N is sufficiently large. The number of elements N required for an effective inner-product
preserving discretization depends on the convergence rate of the numerical integration formula and, typically, under some smoothness assumption of the integrands.
Indeed, if the weight function ω encompasses all the singularities while f (s) and
g(s) are analytic, a Gaussian numerical integration rule with respect ω, essentially, ensures exponential convergence. Note that the weights of Gaussian rules
are always positive and uniformly bounded. See [6] for more details.

10

Yuval Harness

4 Surface Embedding of 2D Random Star-Shaped Obstacles
In this section the main theoretical contribution of this paper is presented. The
first two subsections cover the setting of the problem, where Subsection (4.1)
defines the random shape properties, and Subsection (4.2) covers the relevant
components generalized Polynomial Chaos (gPC) expansion theory. The chosen
framework leads to a reconstruction problem in a Hilbert space. A concise discussion on the disadvantages of naive discretization of the reconstruction problem
concludes Subsection (4.2).
In Subsection (4.3) we present an analytic approach for overcoming the difficulties associated with the naive discretization approach. Using the Coarea formula
we construct a low-dimensional spatial embedding within the family of random
surfaces, which facilitates a natural choice for setting a cubature rule in a compact
region of R2 . The chosen integration weight function is a strictly positive minimal
variance quantity encompassing the irregularities of the family of random surfaces.
In Subsection (4.4) we focus on the case of a single random variable describing the randomness of the object. Using the implicit function theorem we obtain
explicit formulas including full characterization of the singular behavior of the integration weight function. The usage of the single random variable formulation as a
building block for the more general case of multiple random variables is considered
and discussed in Section (6).
Subsections (4.5) and (4.6) are devoted to the demonstration of the preceding
theoretical parts on a model problem of a randomly oriented elliptic cylinder. The
random orientation problem is a very simple ’toy’ problem. However, it allows us
us to demonstrate in a lucid and affable fashion the theoretical ingredients.

4.1 The 2D Random Shape Setting
For brevity, we focus on the sound-soft case and assume that B represents a starshaped obstacle in R2 whose boundary, S, depends smoothly on a real valued
vector of mutually independent and continuous random variables
Z = (Z1 , . . ., ZP ) ,

P ∈ N.

The boundary S is, however, not assumed to be uniformly smooth in the spatial
domain. We assume that each random variable Zp has finite even moments

E

h

i Z
Zp 2n =
IZp

zp2n

dFZp
dzp < ∞ ,
dzp

n ∈ {0, 1, . . ., N } ,

(19)

dFZ (zp )

p
where IZp is the support of Zp and
is the probability density function
dzp
of Zp . The property (19) effectively ensures the existence of surface functionals
suitable for the reconstruction of the scattering coefficients.
Our assumption that the obstacle is star-shaped for any realization of the
random vector Z, ensures that its boundary possesses a polar representation,

S(Z) = { (ρ(θ; Z) · cos θ, ρ(θ; Z) · sin θ)| θ ∈ [0, 2π]} ,

(20)

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

11

and the existence of two positive radial bounds, ρmax and ρmin , satisfying
0 < ρmin = inf ρ(θ; z) < sup ρ(θ; z) = ρmax < ∞ .
θ,z

(21)

θ,z

Thus, as illustrated in Figure (1), S(Z) is confined to the transition region,
o
n
(22)
Rtra = r ∈ R2 ρmin < |r| < ρmax ⊂ R2 .

Fig. 1: The Transition Region Rtra .

4.2 Random Shape and Generalized Polynomial Chaos Expansion
Given our assumptions we observe that the scattering coefficients (7) are finite
dimensional random fields,
bm = bm (Z) .
A common method to approximate these fields is to obtain their generalized Polynomial Chaos (gPC) expansions,
X
bm (Z) ≈ bN
bm,n Pn (Z) ,
(23)
m (Z) =
|n|≤N

where n = (n1 , . . . , nP ) is a multi-index and {Pn } is an orthogonal basis of the
inner-product space induced by the probability density function of Z,
Z
hφ(Z) , ψ(Z)iL2

dFZ

=

φ(z)ψ(z)dFZ (z) ,
IZ

dFZ (z) =

P
Y

dFZp (zp ) ,

k=1

whose support is IZ = IZ1 ×· · ·×IZP . Hence, the expansion coefficients are readily
available by the orthogonality via
bm,n =

1
hbm , Pn iL2 ,
dFZ
γn

γn = hPn , Pn iL2

dFZ

.

12

Yuval Harness

For a randomly shaped obstacle each coefficient in the gPC expansion (23) is
a target functional of the following form
Z Z
bm,n =
gm (r)h(z, ρ) dS(z)Pn (z) dFZ (z) .
(24)
IZ

S(z)

Using the polar form representation
(20) whose associated induced volume form
p
on the surface is dS(θ; z) = ρ2 + ρ2θ dθ, the general representation (24) can be
explicitly written as
Z Z
bm,n =
gm (ρ)h(z, θ)s(θ; z)ρ dθPn (z) dFZ (z) ,
(25)
IZ

[0,2π]

where the normalized polar metric tensor is given by
s
2

1
1 ∂ρ
∂S(θ; z)
·
= 1+
·
.
s(θ; z) =
ρ(θ; z)
∂θ
ρ ∂θ
In principle, we need to devise a discretization scheme for (25) and apply the
optimal reconstruction procedure of Section (3). However, ρ ∈ S(Z) inherits any
irregularity of family of surfaces; e.g., lack of smoothness and oscillatory behaviour,
which often necessitates specialized high-order discretization of the surface S(Z).
Additionally, discretizing the random surface integral with a grid of numerical
integration nodes has to be realized for every grid point in the parameters domain, IZ . Hence, in general, the practical implementation of an inner-product
preserving discretization satisfying (14) is a difficult task. An analytic approach
for overcoming this fundamental difficulty is presented in the next subsection.

4.3 Random Surface Embedding and the Coarea Formula
In this subsection we present an analytic approach for producing inner-product
preserving discretizations of functionals of the form of (25). The key idea is to apply
a change of variables transforming (25) to the following equivalent representation
ZZ
gm (r)φn (h)ω(r)r dθ dr ,
(26)
bm,n =
Rtra

where the weight function ω(r) > 0, is proportional to the conditional expectation of s(θ; z) given the information ρ(θ; z) = r ∈ Rtra . Thus, ω(r) has minimal
variance while, essentially, encompassing the irregularities of the family of random
surfaces, S(Z). The term φn (h) is a linear functional uniformly bounded in Rtra
operating on h. A high-order numerical integration rule with respect to ω(r) would
serve as a discretization satisfying (14). The transformed representation (26) is obtained by the so-called Coarea Formulae [8] which allows us to express the surface
integral in terms of the integral of the level sets of another function.

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

13

Theorem 1 (The Coarea Formula)
Let Φ : Rd+δd → Rd be a piecewise smooth Lipschitz function, where d ≥ 0. Let D
be an open Jordan measurable subset of Rd+δd , then for any integrable function
g : D → R,
!
Z
Z
Z
g(x)
g(x) =
dSΦ−1 (y) (x) dy ,
(27)
D
Rd Φ−1 (y) JΦ (x)
r


where JΦ (x) = det [DΦ (x)] · [DΦ (x)]T is the Jacobian of Φ, and dSΦ−1 (y)
denotes volume density of Φ−1 (y).
Remark 2 The coarea formula expresses the integral in terms of the integral of the
level sets of the function Φ. The level sets, Φ−1 (y), are called fibers of the domain
D. The formula is, in fact, a generalization of Fubini’s theorem.
Let us consider the function
Φ(z, θ) = ρ(θ; z) = (ρ(θ; z) · cos θ, ρ(θ; z) · sin θ) ,

Φ : IZ × [0, 2π] → Rtra .

By direct calculations we obtain

DΦ =

∇z ρ cos θ ∇z ρ sin θ
∂θ (ρ cos θ) ∂θ (ρ sin θ)

T

∈ R2×(P +1) ,

JΦ = ρ · |∇z ρ| ,

where
∇z ρ = (∂z1 ρ, . . . , ∂zP ρ)T ,

|∇z ρ|2 =

P
X

∂zp ρ

2

.

p=1

Hence, applying (27) on (25) with the chosen implicit function, satisfies
ZZ
bm,n =

Z
gm (r)

Rtra

h(z, θ)Pn (z)
IZr

P
s(θ; z) Y dFZp
dSr r dθ dr ,
|∇z ρ| p=1 dzp

(28)

where the domain of integration of the inner integral is the fiber set
IZr = { z ∈ IZ | ρ(θ; z) = r} ⊂ IZ .

(29)

The advantage of the representation (28) is that the wave function, gm , is no longer
composed with the boundary and does not inherit its irregular properties.
The main challenge is to efficiently evaluate the inner integral in (28),
Z
h(z, θ)Pn (z)
I Zr

K
s(θ; z) Y dFZp
dSr .
|∇z ρ|
dzp

(30)

k=1

Q
This integral can become infinite since |∇z ρ|−1 and K
k=1
singular. Accordingly, using the following weight function
Z
ω(r) =
I Zr

K
s(θ; z) Y dFZp
dSr ,
|∇z ρ|
dzp
k=1

dFZp
dzp

are, essentially,

(31)

14

Yuval Harness

we have that
ZZ
bm,n =

gm (r)φn (h)ω(r)r dθ dr ,

(32)

Rtra

where
φn (h) = ω −1 (r) ·

Z
h(z, θ)Pn (z)
IZr

K
s(θ; z) Y dFZp
dSr ,
|∇z ρ|
dzp

(33)

k=1

is a linear functional operating on the surface density h. If φn (h) is bounded then
the choice (31) ensures that (30) is effectively desingularized.
To show that φn (h) is a bounded linear functional, we observe that (31) is, in
fact, proportional to the conditional expectation of s(θ; z) given ρ(θ; z) = r,
ω(r) = E [ s(θ; z)| ρ(θ, z) = r] · fρ (z) ,
where fρ (z) is the probability density function of ρ(θ; Z). It is well known that
conditional expectation is a minimum variance predictor as a function of the given
information. Hence, the choice (31) implies minimization of oscillatory behaviour
of s(θ; z) as a function of ρ(θ; z). Employing a similar argument we obtain that
the linear functional (33) is, in fact,
φn (h) =

E [ h(z, θ)Pn (z)s(θ; z)| ρ(θ; z) = r] .
E [ s(θ; z)| ρ(θ; z) = r]

Hence, the Cauchy-Schwarz inequality implies that


|φn (h)| ≤ max |Pn (z)| · E [ |h(z, θ)| | ρ(θ; z) = r] ,
z∈IZr

which shows that the linear functional, φn , is, indeed, bounded.
Setting the integration grid for (32) can be done in two stages. First we obtain
a cubature rule with nodes {r(i,j) } ⊂ Rtra and corresponding cubature weights
{ω (i,j) } with respect to the weight function ω(r), regardless of z. In the second
stage we identify the integration grid in IZr which corresponds {r(i,j) } for the
evaluation of φn (h).

4.4 Explicit Representation for a Single Random Variable
Let us consider a simplified case of one-dimensional random vector,
Z = Z ∈ IZ ⊂ R .
Then dFZ = dFZ and the target functional (25) reduces to
Z Z
bm,n =
gm (ρ)h(z, θ)s(θ; z)ρ dθPn (z) dFZ (z) .
IZ

(34)

[0,2π]

We will show that in this case the functional φn (h) (33) can be explicitly represented. This approach can serve as a building block for the case of a general
random vector, which is discussed in Section (6).

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

15

For applying formula (28) we observe that the fiber set IZr (29) is, essentially,
a finite set of discrete points; i.e., a zero-dimensional sub-surface. Accordingly, we
obtain the following explicit representation of (28) for a single random variable,

ZZ
Rtra

X

gm (r) 

zr ∈IZr


s(θ; zr ) dFZ
(zr ) r dθ dr ,
h(zr , θ)Pn (zr )
|∂z ρ|z=zr dz

(35)

where IZr is the reduction of IZr to a zero dimensional subset of IZ ,
IZr (r) = {z ∈ IZ | ρ(θ; z) = r } .

(36)

To efficiently evaluate (35), we employ the implicit function theorem [7] which
ensures the following identities
1 = ∂r ρ(θ; zr ) = [∂z ρ(θ; z)]z=zr ∂r zr ,
∂θ ρ(θ; zr ) = [∂z ρ(θ; z)]z=zr ∂θ zr + [∂θ ρ(θ; z)]z=zr .
Thus,
[∂z ρ(θ; z)]z=zr =

1
,
∂r zr

[∂θ ρ(θ; z)]z=zr = −

∂θ zr
,
∂r zr

(37)

and (35) can be equivalently represented by



dF
Z
=
gm (r) 
h(zr , θ)Pn (zr )s(θ; zr )|∂r zr |
(zr ) r dθ dr ,
dz
Rtra
zr ∈IZr


ZZ
X
dF
Z
(zr ) r dθ dr ,
(38)
=
gm (r) 
h(zr , θ)Pn (zr )|∇zr |
dz
Rtra
ZZ

bm,n

X

zr ∈IZr

since the gradient of zr ,
1
∇zr = x̂∂x zr + ŷ∂y zr = r̂∂r zr + θ̂ ∂θ zr ,
r
satisfies |∇zr | = s(θ; zr )|∂r zr |.
By our definitions s(θ; z) is strictly positive. Hence, the zeros of ∂z ρ and the
Z
singularities of dF
dz define the integration rule in the sense, that we can apply a
cubature rule whose weight function captures the singular behaviour of |∂r zr | =
dFZ
|∂z ρ|−1
z=zr and dz . Now, employing (31) yields the following weight function
ω(r) =

X

s(θ; zr )|∂r zr |

zr ∈IZr

X
dFZ
dFZ
(zr ) =
|∇zr |
(zr ) ,
dz
dz
zr ∈IZr

which reduces (38) to
ZZ
bm,n =

gm (r)φn (h)ω(r)r dθ dr ,
Rtra

(39)

16

Yuval Harness

where the bounded linear functional (33) reduces to
φn (h) =

X

h(θ; zr )Pn (zr ) P

zr ∈IZr

Z
|∇zr | dF
dz (zr )

yr ∈IZr

Z
|∇yr | dF
dz (yr )

,
−1


=

X

h(θ; zr )Pn (zr ) 

zr ∈IZr

X

|∇yr |

yr ∈IZr \{zr }

dFZ
(yr )
dz

.

(40)

Once the cubature nodes {r(i,j) } and corresponding weights {ω (i,j) } have been
chosen, we must also identify the value of zr ∈ IZr at these nodes for the evaluation of φn (h) (40). This, essentially, requires the solution of the following convex
minimization problem,
2
o

n
z (i,j,k) = arg min ρ(θ(i,j) ; z) − r(i,j) .
|z|≤1

Thus, in practice we obtain a cubature grid and corresponding weights
n
o
n
o
θ(i,j) , z (i,j,k)
,
ω (i,j) ,
respectively, which apply to the original form (34) in the sense that
X
r(i,j) ω (i,j) ,
gm (r(i,j) )φ(i,j)
bm,n ≈
n
i,j

where r(i,j) = ρ(θ(i,j) , z (i,j,k) ) independently of k, and

−1
X
X
dFZ 
(i,j)
(i,j) (i,j,k)
(i,j,k) 
|∇zr |
φn =
h(θ
;z
)Pn (z
)
dz
q6=k

k

,

z=z (i,j,q)

where |∇zr |z=z(i,j,q) should, generally, be obtained numerically.
4.5 Example: Randomly Oriented Elliptic Cylinder
Let B be a sound-soft elliptic cylinder with major radius a and minor radius b, i.e.,
a > b > 0. The symmetry axis of the cylinder is located at the origin O = (0, 0).
The major and minor axes of the elliptic cross-section are assumed to be rotated
counter-clockwise by z ∈ [0, 2π], see Figure (2). Note, that the radial bounds (21)
are b = ρmin < ρmax = a.
The polar form of the obstacle’s boundary over all random orientation states
is given by
S(z) = { (ρ(θ − z) · cos θ, ρ(θ − z) · sin θ)| θ ∈ [0, 2π]} ,
where ρ(t) ∈ (b, a) satisfies
ρ(t) = p

ab
b2 cos2 t + a2 sin2 t

,

ρ0 (t)
(a2 − b2 ) sin(2t)/2
= 2
.
ρ(t)
b cos2 t + a2 sin2 t

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

17

Fig. 2: Randomly Oriented Elliptic Cylinder.
Thus, the zeros of ∂z ρ(θ − z) are attained at θ − z = 0, π/2, π, 3π/2. Note, that
[ρ(θ − z)]θ−z=0,π = a ,

[ρ(θ − z)]θ−z=π/2,3π/2 = b ,

which are, indeed, the points of the random surface that do not vary in the radial
direction as a function of the parameter z.
The equality ρ(t) = r can be solved analytically which yields
{ρ(t) = r | r ∈ (b, a)} = {−ζ , ζ , π − ζ , −π + ζ} ,
where
ζ(r) = arccos

a
r

(1)

r

r 2 − b2
a2 − b2

(2)

(3)

(41)

!
∈ (0, π/2) ,

∀ r ∈ (b, a) .

(42)

(4)

Hence, we obtain IZr = {zr , zr , zr , zr } (36) where
zr(1) = θ − ζ ,
and

zr(2) = θ + ζ ,

(k)

∂zr
∂r

= (−1)p

zr(3) = θ − ζ + π ,

dζ
dζ
,
= (−1)k+1
dr
dr

zr(4) = θ + ζ − π ,

k = 1, 2, 3, 4 .

Thus, the target functional (38) takes the following form


s
 2
ZZ
4
X
dζ
dF
Z
h(θ, zr(k) )Pn (zr(k) ) 1 +
bm,n =
gm (r) 
(zr(k) ) r dθ dr ,
dr
dz
Rtra
k=1

where
Rtra =

n

r ∈ R2

b < |r| < a

o

.

Let us now assume that Z is a uniformly distributed random variable in [0, 2π],
Z z0
dz
z0
P(Z < z0 ) =
=
.
2π
2π
0
The problem is 2π-periodically smooth in z, hence, it is natural to employ

cos(bn/2cθ) if n = 2bn/2c
P0 (z) = 1 , Pn (z) =
n = 1, 2, . . . ,
sin(bn/2cθ) if n 6= 2bn/2c

18

Yuval Harness

in the gPC expansion (23), which leads to the following representations of the
target functional (38)
a

Z

2π

Z

bm,n =

gm (r)φn (h)ω(r)rdθdr ,
b

(43)

0

where the linear functional is
φn (h) =

4
X

h(θ, zr(k) )Pn (zr(k) ) .

k=1

The weight function (39) is explicitly given by
4
X
1
ω(r) =
2π

s

k=1


1+

dζ
dr

2

ωreg (r)
= p
,
(a − r)(r − b)

(44)

where the regular part of (44) is given by
2
ωreg (r) =
π

s

1
+ (a − r)(r − b) .
r2 (a + r)(r + b)

4.6 Simulation: Randomly Oriented Elliptic Cylinder
In this subsection we explore numerically the randomly oriented elliptic cylinder
example, that was introduced in the previous subsection. First, let us setup an
inner-product preserving discretization of (43). Applying the linear change of varia+b
ables on the radial variable, r(σ) = a−b
2 σ + 2 , which maps (b, a) onto (−1, 1),
we obtain
1
2
1
p
√
=
, σ ∈ (−1, 1) .
a
−
b
1
−
σ2
(a − r)(r − b)
Hence, we can employ the Chebyshev-Gauss quadrature in terms of σ for the
integration in the radial direction,


2i − 1
a − b (i) a + b
(i)
(i)
σ +
, σ = cos
π
i = 1, 2, . . . , M .
r =
2
2
2M
For the angular variable we employ the composite trapezoidal rule,
θ(i,j) = 2

j−1
π j = 1, 2, . . . , N (i) = 10br(i) c ,
N (i)

where N (i) is proportional to r(i) to accommodate for the integration over the
circumference 2πr(i) . The corresponding spatial cubature formula is
Z

2π

θ=0

Z

b

r=a

M N (i)

XX
ψ(r, θ) dr
dθ
2π
p
≈
ψ(r(i) , θ(i,j) ) ·
.
M
·
N (i)
(a − r)(r − b) 2π
i=1 j=1

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

19

Finally, using (41) and (42) we obtain the following expression for the corresponding cubature points in terms of z,
z (i,j,1) = θ(j) − ζ (i)

, z (i,j,2) = θ(j) + ζ (i) ,

z (i,j,3) = θ(j) − ζ (i) + π , z (i,j,4) = θ(j) + ζ (i) − π ,
where


ζ

(i)

a
= arccos  (i)
r

s


r(i)2 − b2 
.
a2 − b2

The resulting spatial grid {(r(i) cos θ(i,j) , r(i) sin θ(i,j) )} in Rtra and the corresponding parametric grid {(θ(i,j) , z (i,j,k) )} in the (θ, z)-plane are displayed in
Figure (3). For comparison we show in Figure (4) a naive discretization of the random surface, whose parametric grid is uniformly distributed in the (θ, z)-plane,
j−1
i−1
(i)
π , j = 1, 2, . . . , N , znaive = 2
π , i = 1, 2, . . . , M .
N
M
Evidently, the corresponding naive spatial grid does a poor job in properly covering
the transition region, Rtra .
(j)

θnaive = 2

(a) the annular (Rtra ) cubature grid.

(b) the (θ, z)-plane cubature grid.

Fig. 3: Elliptic Cylinder: Coarea discretization. Displaying the distribution
of grid points in the spatial and the parametric domains. Gridpoints intensity .
shifts from light at ρmin = b to dark at ρmax = a. (a) spatial gridpoints in the .
annulus Rtra ⊂ R2 . (b) parametric gridpoints in the (θ, z)-plane.
.
For the simulation we consider an elliptic cylinder whose semi-major axis is
a = 5 and whose semi-minor axis is b = 1. We assume an incident plane-wave,
uinc = eiκx = eiκr cos θ =

∞
X
m=−∞

im Jm (κr)eimθ ,

(45)

20

Yuval Harness

(a) the annular (Rtra ) cubature grid.

(b) the (θ, z)-plane cubature grid.

Fig. 4: Elliptic Cylinder: Naive discretization. Displaying the distribution
of grid points in the spatial and the parametric domains. Gridpoints intesity
shifts from light at ρmin = b to dark at ρmax = a. (a) spatial gridpoints in the
annulus Rtra ⊂ R2 . (b) parametric gridpoints in the (θ, z)-plane.

.
.
.
.

which is approximated by truncating the infinite sum in (45) to a finite sum over
the modes |m| ≤ µ where µ satisfies (9). For the discretization we have used
M = 15 and N = 10, and the following thresholds for the reconstruction
evl = 10−4 ,

dis = 10−8 .

b 2 , where G is a matrix whose
Figure (5) displays the construction error kG − Gk
rows are the discretized target functionals,
gm (r) = Jm (κ|r|)e−imθ P0 (z) = Jm (κ|r|)e−imθ ,

|m| ≤ µ ,

b is a matrix whose rows are the corresponding reconstructed target funcand G
tionals from the following set of information functionals
f`,m = Hm (κ |r − r` |)eim^(r−r` ) ,

|m| ≤ M .

The singular points r` ∈ B(Z) are uniformly distributed along the family of random
surface, S(Z),
r` (Z) = 0.95ρ(θ` ; z) · (cos θ` , sin θ` ) ,

θ` = 2

`−1
π,
L

` = 1, 2, . . . , L .

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

21

Fig. 5: Elliptic Cylinder: Reconstruction Error. Displaying the
b vs. log (κ)ρmax for various values of the
reconstruction error kG − Gk
10
parameter defining the number of target functionals, L = 25, 30, . . . , 50.

5 Randomly Shaped Polygonal Cylinders
In this section we consider the application of the theory to the general class of
randomly shaped polygonal cylinders. This class is characterized by non-smooth
randomly varying geometry, and thus serves a proof of concept that the proposed
method can, indeed, be applied to complex shapes. The generalizations to more
complex geometries is discussed in Section (6).

5.1 Piecewise Smooth Polar Form Representation
We consider a star-shape polygonal cylinder in 2D which is given by an ordered
set of points,
ρ0 = (x0 , y0 ) , . . . , ρQ−1 = (xQ−1 , yQ−1 ) , ρQ = ρ0 ,
counter-clockwise distributed in

R2 ,

0 ≤ θq−1 = arctan(yq−1 /xq−1 ) < arctan(yq /xq ) = θq ≤ 2π ,

(46)

that describe the boundary of the polygon where Sq is the line connecting ρq−1
and ρq ,
Q
[

S=
Sq , Sq = ρq−1 + t(ρq − ρq−1 ) t ∈ [0, 1] ,
`=1

for all q = 1, 2, . . . , Q. We also assume that the points ρq are smooth functions of
a random vector Z ∈ IZ , defining simple open differentiable curves that do not
intersect in 2D space. In particular each curve does not cross itself.

22

Yuval Harness

For a polar representation we consider an arbitrary Sq . The equality
(ξ, ψ) = ρ = ρq−1 + t(ρq − ρq−1 ) ,
is equivalent to the system
ξ(z, t) = ρ cos θ = xq−1 + t · (xq − xq−1 ) ,
ψ(z, t) = ρ sin θ = yq−1 + t · (yq − yq−1 ) .
where for every z ∈ IZ , θ(z, t) ∈ [θq−1 , θq ]. Assuming t = t(θ; z) ∈ [0, 1], we obtain
t(θ; z) =

xq−1 sin θ − yq−1 cos θ
ρ(θ; z) sin(θ − θq−1 )
=
·
,
(yq − yq−1 ) cos θ − (xq − xq−1 ) sin θ
ρq (z) sin(θq − θq−1 )

which leads to the following expressions
ρq ρq−1 sin(θq − θq−1 )
,
ρq sin(θq − θ) + ρq−1 sin(θ − θq−1 )
∂ρ
ρq cos(θq − θ) − ρq−1 cos(θ − θq−1 )
(θ; z) = ρ(θ; z) ·
,
∂θ
ρq sin(θq − θ) + ρq−1 sin(θ − θq−1 )
p
where ρq (z) = x2q + yq2 .
ρ(θ; z) =

(47)
(48)

5.2 Considerations for Setting the Spatial Cubature for a Single random Variable
Let Z = Z ∈ [−1, 1] with a given probability density function dFdz(z) , and consider
the functional (38) represented as sum of integrals on the sides of the polygon,


Q ZZ
X
X
dF
Z
gm (r) 
h(θ; zr )Pn (zr ) |∇zr |
bm,n =
(zr ) r dθ dr ,
dz
Rq
q=1
zr ∈IZr ,r∈Rq

=

Q ZZ
X
q=1

(q)
gm (r)φ(q)
(r)r dθ dr ,
n (h)ω

(49)

Rq

where the subdomains of integration are

Rq = r ∈ Rtra ∃z ∈ IZ : r ∈ Sq (z) ⊂ R2 .

(50)

The corresponding weight functions are
ω (q) (r) =

X
zr ∈IZr

|∇zr |

dFZ
(zr ) ,
dz

r ∈ Rq ,

and the bounded linear functional (40) is given by
−1

X
X
dF
Z
(q)
(yr ) ,
φn (h) =
h(θ; zr )Pn (zr ) 
|∇yr |
dz
zr ∈IZr

r ∈ Rq .

yr ∈IZr \{zr }

Note, that by our definitions for any q 6= q 0 , Rq ∩ Rq0 is a zero measure Jordan
set and ∪q Rq is a subset of Rtra but not equal to Rtra (22).

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

23

To define a proper cubature rule we must identify the zeros of ∂z ρ(θ; z) and
their local behavior, i.e., Taylor expansion. Note, that (37) implies that when |∂θ zr |
exists (i.e., a finite and real) then any zero of ∂z ρ(θ; z) is also a zero of ∂θ ρ(θ; z),
which is easier to compute. Indeed, by (48) it is sufficient to solve the equalities
ρq cos(θq − θ) = ρq−1 cos(θ − θq−1 ) ,

q = 1, 2, . . . , Q .

Finally, we note that, in general, ρ(θ; z) as well as ∂z ρ(θ; z) are not smooth as
functions of θ across the angels θq (z) (46). Thus, we must construct a separate
cubature rule for each integral associated with each sub-domain Rq in (49).
5.3 Numerical Example
Consider a star shaped octagonal cylinder whose vertices are given in Table (1),
where a > b > a/2 > 0 and the random variable Z is uniformly distributed in
[−1, 1]. When z = −1 the shape is non-convex star polygon. As z increases the
shape eventually becomes convex and gradually transforms to a full square at
z = 1. See Figure (6) for an illustration. Note, that the radial bounds are
√
√
ρmin = 2(a − b) < 2(a + b) = ρmax .

q

ρq−1 = (xq−1 , yq−1 )

ρq = (xq , yq )

[θq−1 , θq ]

1
2
3
4
5
6
7
8

(a + b(z + 1)/2, 0)
(a + bz, a + bz)
(0, a + b(z + 1)/2)
(−a − bz, a + bz)
(−a − b(z + 1)/2, 0)
(−a − bz, −a − bz)
(0, −a − b(z + 1)/2)
(a + bz, −a − bz)

(a + bz, a + bz)
(0, a + b(z + 1)/2)
(−a − bz, a + bz)
(−a − b(z + 1)/2, 0)
(−a − bz, −a − bz)
(0, −a − b(z + 1)/2)
(a + bz, −a − bz)
(a + b(z + 1)/2, 0)

[0, π/4]
[π/4, π/2]
[π/2, 3π/4]
[3π/4, π]
[π, 5π/4]
[5π/4, 3π/2]
[3π/2, 7π/4]
[7π/4, 2π]

Table 1: Star Shaped Octagonal Cylinder Polar Form Parameters.

The data in Table (6), in fact, satisfies the following formula


π
1 + (−1)q
b
1 − (−1)q √
θq = q , ρq =
2(a + bz) ,
a + (z + 1) +
4
2
2
2
where q ∈ {1, 2, . . . , 8} denotes the index of the segment. Using (47) we have

(a + bz) · a + 2b + 2b z
√

ρ(θ; z) =
,
2(a + bz) · αq (θ) + a + 2b + 2b z · βq (θ)
where

 π

1 + (−1)q
π
π  1 − (−1)q
sin θ − q +
+
sin q − θ ,
2
4
4
2
4
 π
 1 − (−1)q

1 + (−1)q
π
π
βq =
sin q − θ +
sin θ − q +
.
2
4
2
4
4

αq =

24

Yuval Harness

Fig. 6: Polygonal Cylinder: random shape at z = ±1.

By direct calculations we obtain
∂ρ
(θ; z) = bρ2 (θ; z) ·
∂z



βq (θ)
αq (θ)
+
a + bz
2a + b + bz


,

which is strictly positive for all θ ∈ (θq−1 , θq ) and any q = 1, 2, . . . , 8. Hence,
∂z ρ 6= 0 for all (z, θ) ∈ [−1, 1] × [0, 2π].
We observe that in the current example each sub-domain of integration Rq (50)
is a quadrilateral in R2 whose corners are given by ρ(θq−1 , ±1)·(cos θq−1 , sin θq−1 ),
and ρ(θq , ±1) · (cos θq , sin θq ), as illustrated in Figure (7). Hence, we setup the

Fig. 7: Polygonal Cylinder: sub-domains of integration, Rq .

cubature rule by mapping each quadrilateral Rq onto [−1, 1] × [−1, 1] followed by
applying Gauss-Legendre quadratures. The explicit implementation is as follows:

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

25

1. Linearly map the angular segment [θq−1 , θq ] onto [−1, 1],
π
θ(τ ) = (τ + 2q − 1) , τ ∈ [−1, 1] .
8
2. Evaluate the Nq ∈ N Gauss-Legendre quadrature nodes in τ ,
−1 < τ (q,1) < τ (q,2) < · · · < τ (q,Nq ) < 1 .
3. For each θ(q,j) = θ(τ (q,j) ) apply linear map in the radial direction
ρ(θ(q,j) ; 1) − ρ(θ(q,j) ; −1)
ρ(θ(q,j) ; 1) + ρ(θ(q,j) ; −1)
·σ+
.
2
2
4. Evaluate the Mq ∈ N Gauss-Legendre quadrature nodes in σ,
r(q,j) (σ) =

−1 < σ (q,1) < σ (q,2) < · · · < σ (q,Mq ) < 1 .
The overall discretized approximation of (49) becomes
Mq Nq
8 X
X
X

[h(r, zr ) |∂r zr |]r=r(q,i,j) Pn (zr(q,i,j) )r(q,i,j) ω (q,i,j) ,

q=1 i=1 j=1

where
ρ(θ(q,j) ; 1) − ρ(θ(q,j) ; −1) (q,i) ρ(θ(q,j) ; 1) + ρ(θ(q,j) ; −1)
·σ
+
,
2
2
= (r(q,i,j) cos θ(q,j) , r(q,i,j) sin θ(q,j) ) ,

r(q,i,j) =
r(q,i,j)

ω (q,i,j) =

π
ρ(θ(q,j) ; 1) − ρ(θ(q,j) ; −1)
i
i
h
·h
4
0 (σ)
2 )P 0 (τ )
·
(1
−
τ
(1 − s2 )PN
Mq
q

.

τ =τ (q,j)

σ=σ (q,i)

An illustration of the spatial integration nodes admitted by the procedure suggested above is given in Figure (8).
For the simulation we consider the star-shaped polygonal as given in Table (6),
whose parameters are a = 5 and b = 4. We assume an incident plane-wave,
uinc = eiκx = eiκr cos θ =

∞
X

im Jm (κr)eimθ ,

(51)

m=−∞

which is approximated by truncating the infinite sum in (51) to a finite sum over
the modes |m| ≤ µ where µ satisfies (9). For the discretization we have used
Mq ≡ 15 and Nq ≡ 12, and the thresholds evl = 10−4 and dis = 10−8 for the
b 2 , where G is a
reconstruction. Figure (9) displays the construction error kG − Gk
matrix whose rows are the discretized target functionals,
gm (r) = Jm (κ|r|)e−imθ P0 (z) = Jm (κ|r|)e−imθ ,

|m| ≤ µ ,

b is a matrix whose rows are the corresponding reconstructed target funcand G
tionals. For the reconstruction We employed the following information functionals
f`,m = Hm (κ |r − r` |)eim^(r−r` ) ,

|m| ≤ Mq

whose singularities r` = 0.95ρ(θ` ; z) · (cos θ` , sin θ` ) ∈ B(Z) are uniformly distributed along the family of random surface, S(Z); θ` = 2 `−1
L π, ` = 1, 2, . . . , L.

26

Yuval Harness

Fig. 8: Polygonal Cylinder: Coarea discretization. Gridpoints intensity
interchages between light or dark in adjacent subdomains Rq (q = 1, . . . , 8).

Fig. 9: Polygonal Cylinder: Reconstruction Error. Displaying the
b vs. log (κρmax ) for varous numbers of the
reconstruction error kG − Gk
10
parameter defining the number of target functionals, L = 75, 80, . . . , 100.

6 Summary and Future Study
The present paper introduced an alternative approach for quantifying the effects
of random shape in acoustic scattering problems. The core idea of the proposed
method is to construct a spatial embedding within the family of random surfaces
which facilitates a natural choice of a spatial low-dimensional integration rule.
The chosen integration weight function is positive, encompasses random surface
irregularities and of minimal variance. This essentially avoids the fundamental

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

27

problems associated with random surface discretizations, namely lack of smoothness in the proximity of the surface when using a level-set method and strong
non-linear dependence on the variation of the boundary when utilizing random
domain mapping. The efficacy of the proposed approach was demonstrated on
model problems in R2 .
Generalizing the new approach to non star-shaped obstacles as well as to 3D
spatial setting is straightforward. Indeed, given a non star-shaped obstacle we can
represent its random boundary as a union of star-shaped sub-surfaces each with a
local origin. For a star-shaped obstacle in R3 the surface of the obstacle possesses
a spherical representation,
S 3D (Z) = { (ρ cos θ cos ϕ, ρ cos θ sin ϕ, ρ sin θ)| θ ∈ [0, π] , ϕ ∈ [0, 2π]} ,
where ρ = ρ(θ, ϕ; Z). Thus, the gPC expansion of the scattering coefficients (8)
can be represented by
 π 2π

Z
Z
Z Z
P
Y
m
bm
··· 
h(ρ) · gn
(ρ)Pn (z)S(θ, ϕ; z) dθ dϕ
dFZp ,
n,n =
z∈IZ

p=1

θ=0 ϕ=0
|m|

m
where gn
(r) = −jn (κr)Pn (cos ϕ)eimθ , ρ(θ, ϕ; z) = (ρ cos θ cos ϕ, ρ cos θ sin ϕ, ρ sin θ)
and the metric tensor for the spherical case is given by
v
u

2 !
 2
u
∂ρ
∂ρ
t
2
ρ +
sin2 ϕ +
(ρ = ρ(θ, ϕ; z)) .
S(θ, ϕ; Z) = ρ
∂ϕ
∂θ

Employing the Coarea formula with respect to the level sets of Φ(z, θ, ϕ) =
ρ(θ, ϕ; z) yields an expression similar to (28)
ZZZ
m
=
(r)φn (h)ω(r)r2 sin θ dθ dr ,
bm
gn
n,n
Rtra

where the 3D transition region is
(
tra

R

=

r∈R

3

)

inf ρ(θ, ϕ; z) < |r| < sup ρ(θ, ϕ; z)

θ,ϕ,z

θ,ϕ,z

⊂ R3 ,

and the weight function is proportional to the conditional expectetation of the
normalized metric tensor,


S(θ, ϕ; z)
ρ(θ,
ϕ;
z)
=
r
· fρ (z) .
ω(r) = E
ρ2 sin θ
The expressions are, however, technically more complicated to work with. The
usage of automatic integration as well as optimization methods for obtaining the
spatial grid may prove to be a necessity.
Another challenge is to efficiently deal with several random variables. This
can be achieved by employing the single random variable formula (35) as a basic
building block. Indeed, in the 2D case, choosing a single random variable, Zp , and

28

Yuval Harness

applying the coarea formula (28) with respect to the corresponding integration
variable zp yields the following equivalent representation
ZZ
bm,n =
gm (r)φn (h)ω(r)r dθ dr ,
Rtra

where the functional operating on h explicitly satisfies

Z
Z
X
s(θ; z)
φn (h)ω(r) = · · · 
h(z, θ)Pn (z)
∂zp ρ z =z
z ∈I
IZ \IZp

p,r

Zp,r

p

p,r


Y
dFZp
dFZq ,
(zp,r )
dzp
q6=p

and the zero dimensional fiber set corresponding to Zp is

IZp,r = z ∈ IZp ρ(θ; z) = r ⊂ IZp .
Note that zp,r is a function of r as well as z1 , . . . , zp−1 , zp+1 , . . . , zP . In general,
globally using the last formula is not expected to produce an optimal result. A
more sophisticated approach is to partition Rtra into subregions where each subregion is associated with one significant random variable. This, requires performing
sensitivity analysis, similar to the analysis of variance (ANOVA) method [15], that
is also dependent on the spatial coordinates, r. The exploration of this approach
needs a separate extensive study, and deferred to future work.

References
1. A. Ditkowski, Y.H.: Wave scattering by randomly shaped objects. J. Applied Numerical
Mathematics 62, 1819–1836 (2012)
2. Bejan, A.: Shape and Structure, from Engineering to Nature. Shape and Structure, from
Engineering to Nature. Cambridge University Press (2000)
3. Canuto, C., Kozubek, T.: A fictitious domain approach to the numerical solutions of PDEs
in stochastic domains. Numer. Math. 107, 257–293 (2007)
4. Castrillón-Candás, J.E., Nobile, F., Tempone, R.F.: Analytic regularity and collocation
approximation for elliptic pdes with random domain deformations. Computers & Mathematics with Applications 71(6), 1173 – 1197 (2016)
5. Colton, D., Kress, R.: Inverse Acoustic and Electromagnetic Scattering Theory, Second
Edition. Springer (1998)
6. Davis, P.J., Rabinowitz, P.: Methods of numerical integration. Courier Dover Publications
(2007)
7. Dini, U.: Lezioni di analisi infinitesimale. Universitá di Pisa, Pisa, Italy (1887)
8. Federer, H.: Curvature measures. Transactions of the American Mathematical Society
93(3), 418–491 (1959)
9. Ghanem, R.G., Spanos, P.: Stochastic Finite Elements: A Spectral Approach. Springer
Verlag (2002)
10. Golomb, M., Weinberger, H.F.: Optimal approximation and error bounds in Numerical
Approximation. University of Wisconsin Press, Wisconsin Madison (1959)
11. Harbrecht, H., Peters, M., Siebenmorgen, M.: Analysis of the domain mapping method for
elliptic diffusion problems on random domains. Numerische Mathematik 134(4), 823–856
(2016)
12. Harness, Y.: Wave Scattering by Randomly Shaped Obstacles. Thesis, Tel Aviv University
(2013)
13. Harness, Y., Ditkowski, A.: The null-field method: a reconstruction kernel approach. J.
Comput. Phys. 248, 127–146 (2013)
14. Hiptmair, R., Scarabosio, L., Schillings, C., Schwab, C.: Large deformation shape uncertainty quantification in acoustic scattering. Tech. rep., Zürich (2015)

Low-Dimensional Spatial Embedding Method for Shape UQ in Acoustic Scattering

29

15. Holtz, M.: Sparse Grid Quadrature in High Dimensions with Applications in Finance and
Insurance. Lecture Notes in Computational Science and Engineering. Springer Berlin
Heidelberg (2010)
16. Kress, R.: Linear Integral Equations. Applied Mathematical Sciences. Springer New York
(2013)
17. Louis, A.K.: Feature reconstruction in inverse problems. Inverse Problems 27(6)
18. Louis, A.K.: A unified approach to regularization methods for linear ill-posed problems.
Inverse Problems 15, 489–498 (1999)
19. Louis, A.K., Maass, P.: A mollifier method for linear operator equations of the first kind.
Inverse Problems 6, 427–440 (1990)
20. Martin, P.: Multiple Scattering: Interaction of Time-Harmonic Waves with N Obstacles.
Cambridge University Press (2006)
21. Micchelli, C.A., Rivlin, T.J.: A survey of optimal recovery in Optimal Estimation in Approximation Theory. Plenum Press, New York (1985)
22. Nouy, A., Clment, A., Schoefs, F., Moës, N.: An extended stochastic finite element method
for solving stochastic partial differential equations on random domains. Computer Methods
in Applied Mechanics and Engineering 197(51-52), 4663 – 4682 (2008)
23. Nouy, A., Schoefs, F., Moës, N.: X-sfem, a computational technique based on x-fem to
deal with random shapes. European Journal of Computational Mechanics 16(2), 277–293
(2007)
24. Sauter, S., Schwab, C.: Boundary Element Methods. Springer Series in Computational
Mathematics. Springer Berlin Heidelberg (2010)
25. Sethian, J.: Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science. Cambridge
Monographs on Applied and Computational Mathematics. Cambridge University Press
(1999)
26. Tartakovsky, D., Xiu, D.: Numerical methods for differential equations in random domain.
SIAM J. Sci. Comput. 28(3), 1167–1185 (2006)
27. Tartakovsky, D., Xiu, D.: Stochastic analysis of transport in tubes with rough walls. J.
Comput. Phys. 217(1), 248–259 (2006)
28. Wan, X., Xiu, D., Karnidakis, G.E.: Modeling uncertainty in flow simulations via generalized polynomial chaos. J. Comput. Phys. 187, 137–167 (2003)
29. Warnick, K., Chew, W.: Numerical simulation methods for rough surface scattering. Waves
in Random & Complex Media 11(1), 1–30 (2001)
30. Waterman, P.C.: Matrix formulation of electromagnetic scattering. Proceedings of the
IEEE 53(8), 805–812 (1965)
31. Waterman, P.C.: New formulation of acoustic scattering. J. Acoust. Soc. Am. 45(6),
1417–1429 (1969)
32. Wiener, N.: The homogeneous chaos. Amer. J. Math. 60, 897–936 (1938)
33. Wriedt, T.: Light scattering theory and programs: discussion of latest advances and open
problems. Journal of Quantitative Spectroscopy and Radiative Transfer 113, 2465–2469
(2012)
34. Xiu, D.: Numerical Methods for Stochastic Computations: A Spectral Method Approach.
Princeton University Press (2010)
35. Xiu, D., Karnidakis, G.E.: Modeling uncertainty in steady state diffusion problems via
generalized polynomial chaos. Comput. Methods Appl. Math. Eng. 11, 4927–4948 (2002)
36. Xiu, D., Karnidakis, G.E.: The wiener-askey polynomial chaos for stochastic differential
equations. SIAM J. Sci. Comput. 24(2), 26–40 (2002)

