# Effective Parameter Dimension via Bayesian Model Selection in the Inverse Acoustic Scattering Problem

Authors: Abel Palafox; Marcos A. Capistrán; J. Andrés Christen
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We address a prototype inverse scattering problem in the interface of applied mathematics, statistics, and scientific computing. We pose the acoustic inverse scattering problem in a Bayesian inference perspective and simulate from the posterior distribution using MCMC. The PDE forward map is implemented using high performance computing methods. We implement a standard Bayesian model selection method to estimate an effective number of Fourier coefficients that may be retrieved from noisy data within a standard formulation.


---

## Body

## 1. Introduction

Inverse problems are typically ill-posed and analytical solutions are seldom available. Approaches to inverse problems in the interface of applied mathematics, statistics, and scientific computing represent a setting with a myriad of tools for robust solutions, for a number of reasons including its sound theoretical setting for uncertainty quantification [1, 2], prediction, and decision making. See [3] for a recent review. Although these interdisciplinary approaches are becoming prevalent [4–8], topics such as knowledge representation in complex systems [9] and effective dimension are not as mature yet. Sparsity promoting regularization and Bayesian methods have received attention recently [8, 10].In this paper we consider a nonlinear acoustic scattering inverse problem as a case study. The rationale is that the well posedness of the direct problem and robust numerical methods has been comprehensively studied [11, 12]. This fact allows us to address the inverse problem with a reliable direct problem solver.The propagation of acoustic waves in a homogeneous isotropic medium with constant speed of sound is governed by the linear wave equation:(1)
U
t
t
=
c
2
Δ
U
in
R
n
,
n
=
2,3
for a velocity potential U and c the speed of sound in the medium. For monochromatic time-harmonic waves with frequency ω we have
(2)
U
(
x
,
t
)
=
Re
(
e
i
ω
t
u
(
x
)
)
,
where the space dependent term u
(
x
) satisfies the Helmholtz equation. Consider the following:
(3)
Δ
u
+
k
2
u
=
0
,
where k
=
ω
/
c is called thewave number. Given an obstacle of compact support D
⊂
R
n (n
=
2,3), its forward scattering problem is governed by the Helmholtz equation in R
n
-
D
-. The total wave u
(
x
)
=
exp
⁡
(
i
k
x
·
d
)
+
u
s
(
x
) is a superposition of the incident wave exp
⁡
(
i
k
x
·
d
) and the scattered wave u
s
(
x
), and it is subject to a boundary condition on Γ
=
∂
D. The boundary condition may be of type u
=
0 (Dirichlet), ∂
u
/
∂
ν
=
0 (Neumann), or ∂
u
/
∂
ν
+
i
k
λ
u
=
0 (impedance). Also, the scattered waves u
s satisfy the Sommerfeld radiation condition
(4)
lim
⁡
r
→
∞
r
(
n
-
1
)
/
2
(
∂
u
s
∂
r
-
i
k
u
s
)
=
0
,
r
=
∥
x
∥
2
and are referred to asradiating solutions of the Helmholtz equation.Given an incident waveexp
⁡
(
i
k
x
·
d
), the obstacle boundary Γ uniquely determines the scattered wave u
s
(
x
), and every scattered wave u
s
(
x
) which is a radiating solution of the Helmholtz equation has the asymptotic behavior of an outgoing spherical wave. Consider the following:
(5)
u
(
x
)
=
e
i
k
|
x
|
|
x
|
{
u
∞
(
x
^
)
+
O
(
1
|
x
|
)
}
,
where x
^
=
x
/
|
x
| and u
∞
(
x
^
) are the scattering amplitude orfar field pattern. We denote by F
Γ the boundary to far field mapping:
(6)
F
Γ
(
x
^
)
=
u
∞
(
x
^
)
.
The mapping from the obstacle boundary Γ to the far field scattering amplitude u
∞ is injective [11, 12].In this work we address the problem of estimating the boundaryΓ given noisy measurements of the far field data u
∞. This paper is organized as follows. For the sake of making this paper self-contained, in Section 2 we briefly review useful results used throughout the paper. The acoustic layer potential and integral equation method are presented in Section 2.1. We discuss the numerical solution of the forward mapping and its implementation in Section 2.2. In Section 2.3 we present in detail the Bayesian formulation of the inverse problem along with the probability distributions involved. The inverse problem results and the effective dimension analysis are presented in Section 3. We consider the effective dimension exploration as the main contribution of this work.


## 2. Materials and Methods

In this section we describe the forward mapping evaluation using the integral equations method and the layer acoustic potentials approach as in [11, 12].


### 2.1. Single Layer and Double Layer Potentials

The fundamental solution of the Helmholtz equation (3) is
(7)
Φ
(
x
,
y
)
=
{
1
4
π
e
i
k
|
x
-
y
|
|
x
-
y
|
n
=
3
i
4
H
0
(
1
)
(
k
|
x
-
y
|
)
n
=
2
,
for x
≠
y, where |
·
| denotes the L
2 norm and H
0
(
1
) is the Hankel function of first kind and order zero. Given an integrable function φ, the integral operators
(8)
(
S
φ
)
(
x
)
=
∫
Γ
‍
φ
(
y
)
Φ
(
x
,
y
)
d
s
(
y
)
,
x
∈
R
n
∖
Γ
,
(9)
(
K
φ
)
(
x
)
=
∫
Γ
‍
φ
(
y
)
∂
Φ
(
x
,
y
)
∂
ν
(
y
)
d
s
(
y
)
,
x
∈
R
n
∖
Γ
,
with ν being the unit normal vector to Γ directed to the exterior of D, are called acousticsingle layer and acousticdouble layer potentials for density φ, respectively. Both, the single and double layer potentials are radiant solutions of the Helmholtz equation in R
n
∖
D. The single layer potential is continuous in R
n
∖
D (including the boundary Γ) and consequently it is defined as in (8) also for x
∈
Γ. The double layer potential is only continuous in R
n
∖
D
- (not including Γ); nonetheless (
K
φ
)
(
x
) is defined in the boundary as its corresponding limit when x approaches Γ, namely,
(10)
(
K
φ
)
(
x
)
=
1
2
φ
(
x
)
-
∫
Γ
‍
∂
Φ
(
x
,
y
)
∂
ν
(
y
)
φ
(
y
)
d
s
(
y
)
,
x
∈
Γ
;
see [12] for details.The direct problem is formulated in the form of acombined potential using the Dirichlet boundary condition:
(11)
1
2
φ
+
K
φ
-
i
η
S
φ
=
-
u
inc
,
where u
inc
=
exp
⁡
(
i
k
x
·
d
) is the incident wave, η is a coupling factor (in our case we set η
=
k), and S and K are the acoustic potentials defined above.The three potentials (single layer, double layer, and combined potential) reproduce the same far field pattern given a boundaryΓ. However, the combined potential operator in (11) is the sum of the identity and a compact operator, and therefore its singular values are bounded away from zero. Consequently, the combined potential is more stable regarding numerical implementation and then we use the combined potential in the remainder.The far field pattern for the combined potential is given by(12)
u
∞
(
x
^
)
=
e
-
i
(
π
/
4
)
8
π
k
∫
Γ
‍
{
k
ν
(
y
)
·
x
^
+
η
}
e
-
i
k
x
^
·
y
φ
(
y
)
d
s
(
y
)
.
The corresponding far field pattern for the single-layer potential is given by
(13)
u
∞
(
x
^
)
=
e
i
(
π
/
4
)
8
π
k
∫
Γ
‍
e
-
i
k
x
^
·
y
φ
(
y
)
d
s
(
y
)
.The integral (12) can be evaluated numerically using the trapezoidal rule after solving the integral equation (11) for φ. It is known [13] that the trapezoidal rule has high accuracy when integrating over periodic functions and therefore is a sensible method to use in this case.


### 2.2. Numerical Solution of the Forward Map

As a test case we consider the two-dimensional kite-shaped domain shown in Figure1. The parametric representation for the boundary is given by Γ
(
t
)
=
(
cos
⁡
t
+
0.65
cos
⁡
2
t
-
0.65,1.5
sin
t
)
,
0
≤
t
≤
2
π.Figure 1
Synthetic example: a kite-shaped object. Numerical results for forward mapping evaluations are obtained by generating synthetic far field pattern measurements for this smooth, nonconvex boundary.The combined potential equation (11) gives rise to a linear system upon discretization:
(14)
A
φ
=
g
,
where φ is the density, g
=
-
2
u
inc, and the matrix A has the form
(15)
A
=
I
-
L
-
i
k
M
,
where I is the N
×
N identity matrix and N is the number of knots used to discretize the boundary Γ. L and M are the discrete kernels corresponding to single layer potential and double layer potential, respectively (for details see [12]).MatrixA in (15) is square, nonsymmetric, dense, complex, and well conditioned. Mathematical software such as Matlab or Python can be used to evaluate the direct problem. However, it must be taken into account that a big number of evaluations of the forward mapping must be done in a Bayesian statistics approach.Consequently, in order to have an efficient numerical machinery to evaluate the direct problem we have implemented a parallel version of the general conjugate residuals method (GCR) to solve the linear system (14). The GCR method was programmed in C on a graphics processing unit (GPU). Of note, a serial version of the GCR method programmed in C is slow, even compared with matlab and python solvers of the linear systems; see Figure 2. The implementation details for parallel computing scheme are presented on Appendix B. Although the evaluation of the far field pattern (12) is not computationally expensive, it can be done in parallel. Details are shown also in Appendix B.Figure 2
Elapsed times of serial and parallel numerical implementations for a single forward mapping evaluation using the combined potential.The far field pattern valueu
∞
(
d
) obtained for the kite in the directions d
=
(
1,0
) closely resembles results previously reported in [12]. The convergence of Nÿstrom method, used to numerically evaluate the integral operators, is illustrated on Figure 3. We use as a reference solution the far field pattern value u
∞
(
1,0
) with wave number k
=
1, obtained with the combined potential for N
=
1024 points, and we plot in logarithmic scale for y-axis the relative errors for far field patterns obtained with combined potential (solid line in Figure 3) and single layer potential (dashed line in Figure 3) for number of points N
=
8
,
16
,
32
,
64
,
128
,
256.Figure 3
Relative errors of the forward mapping evaluation on directiond
=
(
1,0
), for the smooth and periodic kite-shaped boundary, varying the number of points N and wave number k
=
1 fixed. This plot with logarithmic scale for y-axis illustrates the exponential convergence of Nÿstrom method for combined and single layer potentials.


### 2.3. Bayesian Formulation

First let us remember the Bayes rule:(16)
π
(
θ
∣
u
∞
)
=
π
(
u
∞
∣
θ
)
π
(
θ
)
π
(
u
∞
)
.
Theposterior density
π
(
θ
∣
u
∞
) quantifies the uncertainty of the parameter θ to be recovered from data u
∞. Theprior density
π
(
θ
) expresses the information regarding the unknown parameter θ. π
(
u
∞
∣
θ
) is thelikelihood function, that is, the measurement error distribution assumed, conditional on θ. The density π
(
u
∞
) is a normalizing constant. Below we describe these density functions for our inverse problem.We consider the inverse problem defined by the nonlinear forward mapping (6) from the point of view of Bayesian statistics. In Bayesian statistical inversion theory, the solution of the inverse problem is the conditional probability of parameters, given the data (posterior probability). We assume that the data are noisy measurements of far field pattern u
∞ and that the boundary Γ may be described by a vector of parameters θ as follows.We construct a prior model based on a parameterization forΓ considering that it is a simple 2
π-periodic curve and Γ
∈
C
2
(
R
2
). We express Γ as follows:
(17)
Γ
(
t
)
=
r
(
t
)
(
cos
⁡
t
,
sin
t
)
,
where t
∈
[
0,2
π
) and
(18)
r
(
t
)
=
a
0
2
+
∑
i
=
1
∞
‍
a
i
cos
⁡
(
i
t
)
+
b
i
sin
(
i
t
)
.
Thus the radius of Γ is expressed in terms of a Fourier series. The vector of parameters
(19)
θ
=
(
a
0
,
a
1
,
…
,
a
n
,
…
,
b
1
,
b
2
,
…
,
b
n
,
…
)
contains the Fourier series coefficients in (18). For a boundary Γ that is defined by a vector of coefficients θ we denote the forward mapping by F
θ. We refer to
(20)
θ
n
=
(
a
0
,
a
1
,
…
,
a
n
,
b
1
,
b
2
,
…
,
b
n
)
as the finite vector of c
=
2
n
+
1 Fourier coefficients to represent a boundary Γ
n.


#### 2.3.1. Posterior Distribution

The posterior distribution does not have to our knowledge a closed form. Consequently, we resort toMarkov chain Monte Carlo (MCMC) simulation methods to analyze the posterior distribution. We use the t-walk [14], a general purpose sampling algorithm for continuous distributions, in order to sample satisfactorily from the posterior distribution. We describe briefly the t-walk algorithm and implementation details in Appendix A.The t-walk algorithm evaluates the energy function (minus log of the nonnormalized posterior). Consider the following:(21)
ε
(
θ
)
=
E
-
log
⁡
L
(
θ
,
u
∞
)
-
log
⁡
π
Θ
(
θ
)
,
for an arbitrary constant E (i.e., minus log of the nonnormalized posterior) on each iteration. In practice the constant E is not used and the energy is evaluated as minus log of the product likelihood times the prior one.The output simulations of the t-walk are vectors of Fourier coefficientsθ
n for n fixed, each one of which defines a boundary Γ
n through (18). We refer subsequently to the simulations θ
n, obtained with the t-walk, as MCMC simulations or posterior samples and we refer to asimulated boundary as the curve Γ
n which is defined by a simulated vector θ
n.


#### 2.3.2. Prior Distribution

We establish the prior model as follows. Suppose thatΓ
∈
C
2
(
R
2
) corresponds to a star shaped domain. Then Γ admits a representation (17)-(18). Furthermore, the mean square error for approximating Γ by Γ
~
n is given by
(22)
∥
Γ
-
Γ
~
n
∥
2
=
∑
m
>
n
‍
|
(
a
m
,
b
m
)
|
2
.
Further details of this consequence of Parseval inequality are described in [15].Due to smoothness ofΓ, the amplitude of the Fourier coefficients decays quickly along the series. In fact
(23)
|
(
a
m
,
b
m
)
|
=
O
(
m
-
2
)
.
Note that m is the position of the coefficient in the Fourier series. Therefore the norm |
(
a
m
,
b
m
)
| of the coefficients a
m
,
b
m is by definition of O
(
m
-
2
) bounded by
(24)
|
(
a
m
,
b
m
)
|
<
C
m
2
,
for some positive constant C
∈
R.We pose a normal distribution for each paira
m
,
b
m of coefficients as follows:
(25)
a
m
,
b
m
~
N
(
0
,
C
4
m
2
)
,
with a
m
,
b
m mutually independent. We note that the variance was set to C
/
4
m
2. The rationale is that with the variance in this way, we guarantee that the pair a
m
,
b
m satisfies the inequality (24) with probability higher than 0.99. Moreover the variance decreases when increasing position m, that is, for high order coefficients along the Fourier series. Consequently, this modeling incorporates the smoothness of Γ on the prior distribution.The coefficienta
0 in (18) is related to the size of the scattering object. We assume that
(26)
a
0
~
N
(
2.5
,
0.14
)
.It is not clear how to set the constantC in distribution (25). This constant can be seen as a scaling factor and the parameter m controls the spread of the distribution by 1
∖
4
m
2. A single constant C can be chosen for all prior distributions (25) independently of the value of the index m
=
1
,
…
,
n. For this aim we sample coefficients θ
~
n from the distributions (25) varying the value of C. We add the condition
(27)
r
(
t
)
>
0.75
∀
t
∈
[
0,2
π
)
on the radius (18) to guarantee that the sample curves do not intersect themselves. The sampled curves are presented on Figure 4. As we observe, using the values of C
=
0.5 and C
=
1 (see Figures 4(a) and 4(b)) we obtain very smooth curves. On the other hand, Figures 4(c) and 4(d) show curves with more oscillations. For our purposes, a value of C
=
2 is chosen.Prior distribution samples for coefficients on our boundary representation. The sample is drawn from a Gaussian distribution with meanμ
=
0 and variance σ
2
=
C
/
4
m
2, where m is the position of the coefficient on a Fourier series and C is constant. We vary the variance and set the value of C depending on the expected smoothness of the boundary to recover.
(a)
C
=
0.5
(b)
C
=
1.5
(c)
C
=
2
(d)
C
=
2.5


#### 2.3.3. Likelihood

Assume that the measurements of the far field pattern have additive Gaussian noise(28)
u
∞
(
x
^
)
=
F
Γ
(
θ
)
+
η
x
^
,
with η
x
^
~
N
(
0
,
σ
2
) for a boundary Γ. The likelihood is given by the noise distribution
(29)
L
(
θ
,
u
∞
)
=
(
2
π
σ
)
-
N
/
2
exp
⁡
{
-
1
2
∑
x
^
‍
∥
u
∞
(
x
^
)
-
F
Γ
(
x
^
)
σ
∥
2
}
,
assuming independent measurements. We assume that the measurements u
∞ are taken on evenly spaced directions x
^.


## 2.1. Single Layer and Double Layer Potentials

The fundamental solution of the Helmholtz equation (3) is
(7)
Φ
(
x
,
y
)
=
{
1
4
π
e
i
k
|
x
-
y
|
|
x
-
y
|
n
=
3
i
4
H
0
(
1
)
(
k
|
x
-
y
|
)
n
=
2
,
for x
≠
y, where |
·
| denotes the L
2 norm and H
0
(
1
) is the Hankel function of first kind and order zero. Given an integrable function φ, the integral operators
(8)
(
S
φ
)
(
x
)
=
∫
Γ
‍
φ
(
y
)
Φ
(
x
,
y
)
d
s
(
y
)
,
x
∈
R
n
∖
Γ
,
(9)
(
K
φ
)
(
x
)
=
∫
Γ
‍
φ
(
y
)
∂
Φ
(
x
,
y
)
∂
ν
(
y
)
d
s
(
y
)
,
x
∈
R
n
∖
Γ
,
with ν being the unit normal vector to Γ directed to the exterior of D, are called acousticsingle layer and acousticdouble layer potentials for density φ, respectively. Both, the single and double layer potentials are radiant solutions of the Helmholtz equation in R
n
∖
D. The single layer potential is continuous in R
n
∖
D (including the boundary Γ) and consequently it is defined as in (8) also for x
∈
Γ. The double layer potential is only continuous in R
n
∖
D
- (not including Γ); nonetheless (
K
φ
)
(
x
) is defined in the boundary as its corresponding limit when x approaches Γ, namely,
(10)
(
K
φ
)
(
x
)
=
1
2
φ
(
x
)
-
∫
Γ
‍
∂
Φ
(
x
,
y
)
∂
ν
(
y
)
φ
(
y
)
d
s
(
y
)
,
x
∈
Γ
;
see [12] for details.The direct problem is formulated in the form of acombined potential using the Dirichlet boundary condition:
(11)
1
2
φ
+
K
φ
-
i
η
S
φ
=
-
u
inc
,
where u
inc
=
exp
⁡
(
i
k
x
·
d
) is the incident wave, η is a coupling factor (in our case we set η
=
k), and S and K are the acoustic potentials defined above.The three potentials (single layer, double layer, and combined potential) reproduce the same far field pattern given a boundaryΓ. However, the combined potential operator in (11) is the sum of the identity and a compact operator, and therefore its singular values are bounded away from zero. Consequently, the combined potential is more stable regarding numerical implementation and then we use the combined potential in the remainder.The far field pattern for the combined potential is given by(12)
u
∞
(
x
^
)
=
e
-
i
(
π
/
4
)
8
π
k
∫
Γ
‍
{
k
ν
(
y
)
·
x
^
+
η
}
e
-
i
k
x
^
·
y
φ
(
y
)
d
s
(
y
)
.
The corresponding far field pattern for the single-layer potential is given by
(13)
u
∞
(
x
^
)
=
e
i
(
π
/
4
)
8
π
k
∫
Γ
‍
e
-
i
k
x
^
·
y
φ
(
y
)
d
s
(
y
)
.The integral (12) can be evaluated numerically using the trapezoidal rule after solving the integral equation (11) for φ. It is known [13] that the trapezoidal rule has high accuracy when integrating over periodic functions and therefore is a sensible method to use in this case.


## 2.2. Numerical Solution of the Forward Map

As a test case we consider the two-dimensional kite-shaped domain shown in Figure1. The parametric representation for the boundary is given by Γ
(
t
)
=
(
cos
⁡
t
+
0.65
cos
⁡
2
t
-
0.65,1.5
sin
t
)
,
0
≤
t
≤
2
π.Figure 1
Synthetic example: a kite-shaped object. Numerical results for forward mapping evaluations are obtained by generating synthetic far field pattern measurements for this smooth, nonconvex boundary.The combined potential equation (11) gives rise to a linear system upon discretization:
(14)
A
φ
=
g
,
where φ is the density, g
=
-
2
u
inc, and the matrix A has the form
(15)
A
=
I
-
L
-
i
k
M
,
where I is the N
×
N identity matrix and N is the number of knots used to discretize the boundary Γ. L and M are the discrete kernels corresponding to single layer potential and double layer potential, respectively (for details see [12]).MatrixA in (15) is square, nonsymmetric, dense, complex, and well conditioned. Mathematical software such as Matlab or Python can be used to evaluate the direct problem. However, it must be taken into account that a big number of evaluations of the forward mapping must be done in a Bayesian statistics approach.Consequently, in order to have an efficient numerical machinery to evaluate the direct problem we have implemented a parallel version of the general conjugate residuals method (GCR) to solve the linear system (14). The GCR method was programmed in C on a graphics processing unit (GPU). Of note, a serial version of the GCR method programmed in C is slow, even compared with matlab and python solvers of the linear systems; see Figure 2. The implementation details for parallel computing scheme are presented on Appendix B. Although the evaluation of the far field pattern (12) is not computationally expensive, it can be done in parallel. Details are shown also in Appendix B.Figure 2
Elapsed times of serial and parallel numerical implementations for a single forward mapping evaluation using the combined potential.The far field pattern valueu
∞
(
d
) obtained for the kite in the directions d
=
(
1,0
) closely resembles results previously reported in [12]. The convergence of Nÿstrom method, used to numerically evaluate the integral operators, is illustrated on Figure 3. We use as a reference solution the far field pattern value u
∞
(
1,0
) with wave number k
=
1, obtained with the combined potential for N
=
1024 points, and we plot in logarithmic scale for y-axis the relative errors for far field patterns obtained with combined potential (solid line in Figure 3) and single layer potential (dashed line in Figure 3) for number of points N
=
8
,
16
,
32
,
64
,
128
,
256.Figure 3
Relative errors of the forward mapping evaluation on directiond
=
(
1,0
), for the smooth and periodic kite-shaped boundary, varying the number of points N and wave number k
=
1 fixed. This plot with logarithmic scale for y-axis illustrates the exponential convergence of Nÿstrom method for combined and single layer potentials.


## 2.3. Bayesian Formulation

First let us remember the Bayes rule:(16)
π
(
θ
∣
u
∞
)
=
π
(
u
∞
∣
θ
)
π
(
θ
)
π
(
u
∞
)
.
Theposterior density
π
(
θ
∣
u
∞
) quantifies the uncertainty of the parameter θ to be recovered from data u
∞. Theprior density
π
(
θ
) expresses the information regarding the unknown parameter θ. π
(
u
∞
∣
θ
) is thelikelihood function, that is, the measurement error distribution assumed, conditional on θ. The density π
(
u
∞
) is a normalizing constant. Below we describe these density functions for our inverse problem.We consider the inverse problem defined by the nonlinear forward mapping (6) from the point of view of Bayesian statistics. In Bayesian statistical inversion theory, the solution of the inverse problem is the conditional probability of parameters, given the data (posterior probability). We assume that the data are noisy measurements of far field pattern u
∞ and that the boundary Γ may be described by a vector of parameters θ as follows.We construct a prior model based on a parameterization forΓ considering that it is a simple 2
π-periodic curve and Γ
∈
C
2
(
R
2
). We express Γ as follows:
(17)
Γ
(
t
)
=
r
(
t
)
(
cos
⁡
t
,
sin
t
)
,
where t
∈
[
0,2
π
) and
(18)
r
(
t
)
=
a
0
2
+
∑
i
=
1
∞
‍
a
i
cos
⁡
(
i
t
)
+
b
i
sin
(
i
t
)
.
Thus the radius of Γ is expressed in terms of a Fourier series. The vector of parameters
(19)
θ
=
(
a
0
,
a
1
,
…
,
a
n
,
…
,
b
1
,
b
2
,
…
,
b
n
,
…
)
contains the Fourier series coefficients in (18). For a boundary Γ that is defined by a vector of coefficients θ we denote the forward mapping by F
θ. We refer to
(20)
θ
n
=
(
a
0
,
a
1
,
…
,
a
n
,
b
1
,
b
2
,
…
,
b
n
)
as the finite vector of c
=
2
n
+
1 Fourier coefficients to represent a boundary Γ
n.


### 2.3.1. Posterior Distribution

The posterior distribution does not have to our knowledge a closed form. Consequently, we resort toMarkov chain Monte Carlo (MCMC) simulation methods to analyze the posterior distribution. We use the t-walk [14], a general purpose sampling algorithm for continuous distributions, in order to sample satisfactorily from the posterior distribution. We describe briefly the t-walk algorithm and implementation details in Appendix A.The t-walk algorithm evaluates the energy function (minus log of the nonnormalized posterior). Consider the following:(21)
ε
(
θ
)
=
E
-
log
⁡
L
(
θ
,
u
∞
)
-
log
⁡
π
Θ
(
θ
)
,
for an arbitrary constant E (i.e., minus log of the nonnormalized posterior) on each iteration. In practice the constant E is not used and the energy is evaluated as minus log of the product likelihood times the prior one.The output simulations of the t-walk are vectors of Fourier coefficientsθ
n for n fixed, each one of which defines a boundary Γ
n through (18). We refer subsequently to the simulations θ
n, obtained with the t-walk, as MCMC simulations or posterior samples and we refer to asimulated boundary as the curve Γ
n which is defined by a simulated vector θ
n.


### 2.3.2. Prior Distribution

We establish the prior model as follows. Suppose thatΓ
∈
C
2
(
R
2
) corresponds to a star shaped domain. Then Γ admits a representation (17)-(18). Furthermore, the mean square error for approximating Γ by Γ
~
n is given by
(22)
∥
Γ
-
Γ
~
n
∥
2
=
∑
m
>
n
‍
|
(
a
m
,
b
m
)
|
2
.
Further details of this consequence of Parseval inequality are described in [15].Due to smoothness ofΓ, the amplitude of the Fourier coefficients decays quickly along the series. In fact
(23)
|
(
a
m
,
b
m
)
|
=
O
(
m
-
2
)
.
Note that m is the position of the coefficient in the Fourier series. Therefore the norm |
(
a
m
,
b
m
)
| of the coefficients a
m
,
b
m is by definition of O
(
m
-
2
) bounded by
(24)
|
(
a
m
,
b
m
)
|
<
C
m
2
,
for some positive constant C
∈
R.We pose a normal distribution for each paira
m
,
b
m of coefficients as follows:
(25)
a
m
,
b
m
~
N
(
0
,
C
4
m
2
)
,
with a
m
,
b
m mutually independent. We note that the variance was set to C
/
4
m
2. The rationale is that with the variance in this way, we guarantee that the pair a
m
,
b
m satisfies the inequality (24) with probability higher than 0.99. Moreover the variance decreases when increasing position m, that is, for high order coefficients along the Fourier series. Consequently, this modeling incorporates the smoothness of Γ on the prior distribution.The coefficienta
0 in (18) is related to the size of the scattering object. We assume that
(26)
a
0
~
N
(
2.5
,
0.14
)
.It is not clear how to set the constantC in distribution (25). This constant can be seen as a scaling factor and the parameter m controls the spread of the distribution by 1
∖
4
m
2. A single constant C can be chosen for all prior distributions (25) independently of the value of the index m
=
1
,
…
,
n. For this aim we sample coefficients θ
~
n from the distributions (25) varying the value of C. We add the condition
(27)
r
(
t
)
>
0.75
∀
t
∈
[
0,2
π
)
on the radius (18) to guarantee that the sample curves do not intersect themselves. The sampled curves are presented on Figure 4. As we observe, using the values of C
=
0.5 and C
=
1 (see Figures 4(a) and 4(b)) we obtain very smooth curves. On the other hand, Figures 4(c) and 4(d) show curves with more oscillations. For our purposes, a value of C
=
2 is chosen.Prior distribution samples for coefficients on our boundary representation. The sample is drawn from a Gaussian distribution with meanμ
=
0 and variance σ
2
=
C
/
4
m
2, where m is the position of the coefficient on a Fourier series and C is constant. We vary the variance and set the value of C depending on the expected smoothness of the boundary to recover.
(a)
C
=
0.5
(b)
C
=
1.5
(c)
C
=
2
(d)
C
=
2.5


### 2.3.3. Likelihood

Assume that the measurements of the far field pattern have additive Gaussian noise(28)
u
∞
(
x
^
)
=
F
Γ
(
θ
)
+
η
x
^
,
with η
x
^
~
N
(
0
,
σ
2
) for a boundary Γ. The likelihood is given by the noise distribution
(29)
L
(
θ
,
u
∞
)
=
(
2
π
σ
)
-
N
/
2
exp
⁡
{
-
1
2
∑
x
^
‍
∥
u
∞
(
x
^
)
-
F
Γ
(
x
^
)
σ
∥
2
}
,
assuming independent measurements. We assume that the measurements u
∞ are taken on evenly spaced directions x
^.


## 2.3.1. Posterior Distribution

The posterior distribution does not have to our knowledge a closed form. Consequently, we resort toMarkov chain Monte Carlo (MCMC) simulation methods to analyze the posterior distribution. We use the t-walk [14], a general purpose sampling algorithm for continuous distributions, in order to sample satisfactorily from the posterior distribution. We describe briefly the t-walk algorithm and implementation details in Appendix A.The t-walk algorithm evaluates the energy function (minus log of the nonnormalized posterior). Consider the following:(21)
ε
(
θ
)
=
E
-
log
⁡
L
(
θ
,
u
∞
)
-
log
⁡
π
Θ
(
θ
)
,
for an arbitrary constant E (i.e., minus log of the nonnormalized posterior) on each iteration. In practice the constant E is not used and the energy is evaluated as minus log of the product likelihood times the prior one.The output simulations of the t-walk are vectors of Fourier coefficientsθ
n for n fixed, each one of which defines a boundary Γ
n through (18). We refer subsequently to the simulations θ
n, obtained with the t-walk, as MCMC simulations or posterior samples and we refer to asimulated boundary as the curve Γ
n which is defined by a simulated vector θ
n.


## 2.3.2. Prior Distribution

We establish the prior model as follows. Suppose thatΓ
∈
C
2
(
R
2
) corresponds to a star shaped domain. Then Γ admits a representation (17)-(18). Furthermore, the mean square error for approximating Γ by Γ
~
n is given by
(22)
∥
Γ
-
Γ
~
n
∥
2
=
∑
m
>
n
‍
|
(
a
m
,
b
m
)
|
2
.
Further details of this consequence of Parseval inequality are described in [15].Due to smoothness ofΓ, the amplitude of the Fourier coefficients decays quickly along the series. In fact
(23)
|
(
a
m
,
b
m
)
|
=
O
(
m
-
2
)
.
Note that m is the position of the coefficient in the Fourier series. Therefore the norm |
(
a
m
,
b
m
)
| of the coefficients a
m
,
b
m is by definition of O
(
m
-
2
) bounded by
(24)
|
(
a
m
,
b
m
)
|
<
C
m
2
,
for some positive constant C
∈
R.We pose a normal distribution for each paira
m
,
b
m of coefficients as follows:
(25)
a
m
,
b
m
~
N
(
0
,
C
4
m
2
)
,
with a
m
,
b
m mutually independent. We note that the variance was set to C
/
4
m
2. The rationale is that with the variance in this way, we guarantee that the pair a
m
,
b
m satisfies the inequality (24) with probability higher than 0.99. Moreover the variance decreases when increasing position m, that is, for high order coefficients along the Fourier series. Consequently, this modeling incorporates the smoothness of Γ on the prior distribution.The coefficienta
0 in (18) is related to the size of the scattering object. We assume that
(26)
a
0
~
N
(
2.5
,
0.14
)
.It is not clear how to set the constantC in distribution (25). This constant can be seen as a scaling factor and the parameter m controls the spread of the distribution by 1
∖
4
m
2. A single constant C can be chosen for all prior distributions (25) independently of the value of the index m
=
1
,
…
,
n. For this aim we sample coefficients θ
~
n from the distributions (25) varying the value of C. We add the condition
(27)
r
(
t
)
>
0.75
∀
t
∈
[
0,2
π
)
on the radius (18) to guarantee that the sample curves do not intersect themselves. The sampled curves are presented on Figure 4. As we observe, using the values of C
=
0.5 and C
=
1 (see Figures 4(a) and 4(b)) we obtain very smooth curves. On the other hand, Figures 4(c) and 4(d) show curves with more oscillations. For our purposes, a value of C
=
2 is chosen.Prior distribution samples for coefficients on our boundary representation. The sample is drawn from a Gaussian distribution with meanμ
=
0 and variance σ
2
=
C
/
4
m
2, where m is the position of the coefficient on a Fourier series and C is constant. We vary the variance and set the value of C depending on the expected smoothness of the boundary to recover.
(a)
C
=
0.5
(b)
C
=
1.5
(c)
C
=
2
(d)
C
=
2.5


## 2.3.3. Likelihood

Assume that the measurements of the far field pattern have additive Gaussian noise(28)
u
∞
(
x
^
)
=
F
Γ
(
θ
)
+
η
x
^
,
with η
x
^
~
N
(
0
,
σ
2
) for a boundary Γ. The likelihood is given by the noise distribution
(29)
L
(
θ
,
u
∞
)
=
(
2
π
σ
)
-
N
/
2
exp
⁡
{
-
1
2
∑
x
^
‍
∥
u
∞
(
x
^
)
-
F
Γ
(
x
^
)
σ
∥
2
}
,
assuming independent measurements. We assume that the measurements u
∞ are taken on evenly spaced directions x
^.


## 3. Results and Discussion

In order to avoid committing an “inverse crime” [5] in our numerical experiments we proceed as follows. First we generate synthetic noisy far field measurements from the kite-shaped object (see Figure 1). Second, this kite is not obtained from any finite Fourier series expansion and it is therefore not included in any of our models. This forces our methodology to cope with diverse shapes, being the finite Fourier series only an approximation. Third, synthetic data generation is done using the single-layer potential approach, whereas for MCMC steps we evaluate the forward mapping using the combined potential. Of note, both the single and the combined potential approaches approximate the same far field pattern when N is large enough, however they are numerically different.We present results of the MCMC simulations, varying the number of Fourier coefficients in Section3.1. Results representing the main contribution of this paper are presented on Section 3.2.


### 3.1. Inverse Problem Results

For our experiments we discretize the boundary of the kite-shaped object with 256 points. The far field pattern was computed on 16 evenly spaced points. We use 16 incident waves and we set the coupling parameterη to be equal to the wave number k
=
1. The synthetic data are generated with additive Gaussian noise N
(
0
,
σ
2
) whit σ
=
0.1 (equivalently SNR = 0.97). We generate 100,000 MCMC simulations of the posterior distribution for each number of coefficients (c
=
3,5
,
…
,
31,33). The transient stage of the MCMC is taken into account and we remove the first 5,000 iterations in each case. We highlight that it takes about 7 hours to produce 100,000 MCMC simulations on a Python implementation. With our parallel implementation the execution time is reduced to about 1 hour.Due to space constrains, we present only marginal posterior probability distributions of cosine terms forc
=
11 coefficients in Figure 5. As we observe, for each coefficient, the corresponding distribution is unimodal with low variation. The same feature is presented for coefficients for cases c
∈
{
3,5
,
…
,
31,33
}
,
c
≠
11 (not shown).Marginal posterior distributions fora
i’s coefficients for c
=
11. These unimodal distributions are used to estimate coefficients of the inverse problem.
(a)
a
0
(b)
a
1
(c)
a
2
(d)
a
3
(e)
a
4
(f)
a
5We present in Figure6 the probability region (in gray) for the number of Fourier coefficients for cases c
∈
{
3,5
,
…
,
31,33
}. For this aim we draw 1000 simulated boundaries subsampled with their corresponding IAT factor in each case. The true kite-shaped boundary is shown as a red dashed line. From Figure 6 it is apparent that as we increase the number of coefficients the samples have more oscillations and the probability region becomes imprecise. This fact exhibits the numerical instability of Fourier-based representation for the solutions of the inverse problem of interest.Probability Regions from the MCMC simulations (gray). The boundaries defined from mean of marginal posterior distribution are shown in blue. The corresponding MAP boundaries are shown in green. The results are compared to the kite object in dashed red line.
(a)
(b)Also, in Figure6 we show what we refer to as MAPboundary (dashed line) andmean boundary (continuous line) for each value of c. The MAP boundary is the curve defined by the simulated vector θ
n that has the lowest energy value (maximum a posteriori). The mean boundary is the curve defined by the mean of each coefficient of simulated vectors θ
n. MAP boundaries and mean boundaries seem to approximate the true curve with approximately the same quality. In fact, we cannot qualitatively distinguish the best approximation from c
=
9 coefficients to c
=
33. This issue gives rise to the analysis presented on the next section.


### 3.2. Effective Dimension

We refer toeffective dimension as the number of parameters that can be properly retrieved from a noisy data set. In the context of this paper, effective dimension is related to the uncertainty quantification of parameters, that is, identifying the number of Fourier coefficients on the representation given by (17)–(20), required to approximate properly the boundary of the kite (see Figure 1) from synthetic far field measurements. As an aside comment we want to mention that effective dimension can be seen as an indicator for a parsimonious approximation regarding computational cost. However this feature is not used here due to the fact that effective dimension analysis is performed as a posterior procedure. In this section we have changed the notation for the distributions to be consistent with [14].We define a collection ofk models given by
(30)
π
u
∞
∣
Θ
n
n
(
u
∞
∣
θ
n
)
,
with n
∈
{
1,2
,
…
,
k
}. Θ
n is a vector of random variables. A realization of Θ
n is a vector θ
n defined in (20), and π
u
∞
∣
Θ
n
n is the corresponding posterior distribution for model n. For simplicity we use as index n to indicate the model corresponding to c
=
2
n
+
1 Fourier coefficients.We use the super-model approach; that is, a new model is defined(31)
π
u
∞
∣
Θ
,
M
(
u
∞
∣
θ
,
n
)
=
π
u
∞
∣
Θ
n
n
(
u
∞
∣
θ
n
)
,
where M
=
n is the indicator of the model and Θ
=
(
Θ
1
,
Θ
2
,
…
,
Θ
k
).The posterior probability for modelM
=
n of explaining the observed far field pattern is given by
(32)
π
(
M
=
n
∣
u
∞
)
∝
π
M
(
n
)
π
u
∞
n
(
u
∞
)
,
where π
M
(
n
) is the prior probability of model M
=
n and π
u
∞
n
(
u
∞
) is the corresponding normalizing constant (also known as the marginal likelihood) for model n. A straightforward way to perform model selection is by computing (32) for all the models and choosing the model with highest probability. Indeed, the difficulty regarding the computation of (32) is evaluating the normalizing constant π
u
∞
n
(
u
∞
). To perform our MCMC only the nonnormalized posterior is needed and therefore π
u
∞
n
(
u
∞
) was not previously calculated (this is the usual case when doing MCMC in a fixed number of parameters in a Bayesian application). Nevertheless, the MCMC simulations may be used to estimate π
u
∞
n
(
u
∞
) and provide information to explore the effective dimension. This is explained below.The ratio of two normalizing constants is defined as follows:(33)
r
=
π
u
∞
l
(
u
∞
)
π
u
∞
s
(
u
∞
)
,
where s
,
l
∈
{
1,2
,
…
,
k
}
,
s
≠
l. Considering a uniform prior distribution for the models, that is, π
M
(
n
)
=
1
∖
n, model s is better than model l when r
>
1, and they are indistinguishable if r
=
1 and l is better than s when r
<
1. Then, the best model has the highest normalizing constant value.A consistent estimator of the ratio is(34)
r
^
=
(
1
/
m
l
)
∑
i
=
1
m
l
‍
π
u
∞
∣
Θ
l
l
(
u
∞
∣
θ
l
,
i
)
π
Θ
l
l
(
θ
l
,
i
)
/
π
l
I
(
θ
l
,
i
)
(
1
/
m
s
)
∑
i
=
1
m
s
‍
π
u
∞
∣
Θ
s
s
(
u
∞
∣
θ
s
,
i
)
π
Θ
s
s
(
θ
s
,
i
)
/
π
s
I
(
θ
s
,
i
)
,
where π
l
I
,
π
s
I are completely known importance sampling densities for π
Θ
l
|
u
∞
l and π
Θ
s
|
u
∞
s, respectively, and the sets {
θ
l
,
1
,
θ
l
,
2
,
…
,
θ
l
,
m
l
} and {
θ
s
,
1
,
θ
s
,
2
,
…
,
θ
s
,
m
s
} are samples of the posterior distribution for models l and s, respectively, (see [14] for details).We observe from (34) that for a model s, the product
(35)
π
u
∞
∣
Θ
s
s
(
u
∞
∣
θ
s
,
i
)
π
Θ
s
s
(
θ
s
,
i
)
for
each
θ
s
,
i
∈
{
θ
s
,
1
,
θ
s
,
2
,
…
,
θ
s
,
m
s
}
,
is the corresponding likelihood times prior distribution. This product is in fact the energy (21) using a suitable constant E. Then the resulting sample from MCMC posterior simulations is used to compute (34), as a byproduct of the former with marginal computational cost added. Some details about approximating normalizing constant from MCMC simulations are discussed in [16].The importance sampling densityπ
s
I must be comparable to the product (35). We propose to numerically estimate the density π
s
I for every model using a Kernel density estimation (KDE) [17] of the posterior distribution using the MCMC sample. A rough estimate may be used, for the sole purpose of estimating r.The kernel density estimator has the form(36)
π
^
s
(
θ
s
*
∣
u
∞
)
=
1
m
s
h
m
s
∑
i
=
1
m
s
‍
K
s
(
θ
s
*
-
θ
s
,
i
h
m
s
i
)
,
where the kernel K is a bounded density on R
s and h
m
s is the bandwidth. For this estimator we define a Gaussian kernel
(37)
K
s
(
x
)
=
(
2
π
)
-
s
/
2
e
-
x
2
/
2
and a bandwidth
(38)
h
m
s
=
∏
i
=
1
m
s
‍
h
m
s
i
,
where h
m
s
i
=
1.06
σ
i
m
s
-
1
/
5 and σ
i is the standard deviation for each coefficient on the model s. We compute
(39)
c
s
=
(
1
m
s
)
∑
i
=
1
m
s
‍
π
u
∞
∣
Θ
s
s
(
u
∞
∣
θ
s
,
i
)
π
Θ
s
s
(
θ
s
,
i
)
π
s
I
(
θ
s
,
i
)
,
for each model s
∈
{
1,2
,
…
,
k
} then all the models can be compared together. Then effective dimension corresponds to the model with the highest normalizing constant, equivalently the lowest value of minus logarithm of the normalizing constant.We present in Figure7 the minus logarithm of the estimated normalizing constants (39) for all models (red line). As we see in this figure, the lowest value is obtained with n
=
11 coefficients. Also we present in the same Figure the distance (norm L
2 of the pointwise difference) between the radius of the kite and the radius of the simulated boundaries. The continuous line shows the difference with respect to the mean boundaries. The dashed line corresponds to the distances with respect to the map boundaries. The boxplot in the same figure corresponds to distances of the last 1000 simulations of the posterior. These last simulations are taken and subsampled with their corresponding IAT factor. These boundaries are in fact the probability region shown in Figure 6.Figure 7
Minus logarithm of estimated normalizing constants for Bayesian model comparisons. Best model refers to the highest normalizing constant (c
=
11 for this case). Continuous line shows the distance between the kite and the mean boundaries. Dashed line corresponds to the distance between the kite and the MAP. The boxplot corresponds to distances of the last 1000 simulations of the posterior to the kite. Red scale corresponds to minus logarithm of normalizing constants and blue scale is for the distances.Based on the Bayes factors, the best model corresponds ton
=
11 (we recall that c
=
2
n
+
1). Of note, there are MAP boundaries (e.g. c
∈
{
13
,
19
,
25
,
29
}) with lower mean squared error compared with the MAP for c
=
11. Also the mean boundaries from c
=
15 to c
=
33 have lower distances to the MAP. However as we observe on the boxplot, the mean of the error decreases from 3 to 11 Fourier coefficients and then oscillations appear, increasing spread and increasing error from 13 to 33 coefficients. Despite the low error of MAP and mean boundaries for models with 13 to 33 coefficients, due to the spread shown in the boxplot, the best model is certainly given by Bayes’ factor (i.e., for c
=
11). A qualitative analysis that agrees with these results is also discussed in Section 3.Fourier series is a classical representation for smooth periodic closed curves. Furthermore, we have used a well understood and high order numerical method for forward mapping evaluation. However, these two elements by themselves are not enough to have confidence regarding the solution of the inverse problem. The quality of the solution depends also on the number of coefficients of the Fourier series used for the representation.Our effective dimension discussion relies on MCMC methods to provide a quantitative strategy to determine the number of parameters that can be retrieved given a data set.


## 3.1. Inverse Problem Results

For our experiments we discretize the boundary of the kite-shaped object with 256 points. The far field pattern was computed on 16 evenly spaced points. We use 16 incident waves and we set the coupling parameterη to be equal to the wave number k
=
1. The synthetic data are generated with additive Gaussian noise N
(
0
,
σ
2
) whit σ
=
0.1 (equivalently SNR = 0.97). We generate 100,000 MCMC simulations of the posterior distribution for each number of coefficients (c
=
3,5
,
…
,
31,33). The transient stage of the MCMC is taken into account and we remove the first 5,000 iterations in each case. We highlight that it takes about 7 hours to produce 100,000 MCMC simulations on a Python implementation. With our parallel implementation the execution time is reduced to about 1 hour.Due to space constrains, we present only marginal posterior probability distributions of cosine terms forc
=
11 coefficients in Figure 5. As we observe, for each coefficient, the corresponding distribution is unimodal with low variation. The same feature is presented for coefficients for cases c
∈
{
3,5
,
…
,
31,33
}
,
c
≠
11 (not shown).Marginal posterior distributions fora
i’s coefficients for c
=
11. These unimodal distributions are used to estimate coefficients of the inverse problem.
(a)
a
0
(b)
a
1
(c)
a
2
(d)
a
3
(e)
a
4
(f)
a
5We present in Figure6 the probability region (in gray) for the number of Fourier coefficients for cases c
∈
{
3,5
,
…
,
31,33
}. For this aim we draw 1000 simulated boundaries subsampled with their corresponding IAT factor in each case. The true kite-shaped boundary is shown as a red dashed line. From Figure 6 it is apparent that as we increase the number of coefficients the samples have more oscillations and the probability region becomes imprecise. This fact exhibits the numerical instability of Fourier-based representation for the solutions of the inverse problem of interest.Probability Regions from the MCMC simulations (gray). The boundaries defined from mean of marginal posterior distribution are shown in blue. The corresponding MAP boundaries are shown in green. The results are compared to the kite object in dashed red line.
(a)
(b)Also, in Figure6 we show what we refer to as MAPboundary (dashed line) andmean boundary (continuous line) for each value of c. The MAP boundary is the curve defined by the simulated vector θ
n that has the lowest energy value (maximum a posteriori). The mean boundary is the curve defined by the mean of each coefficient of simulated vectors θ
n. MAP boundaries and mean boundaries seem to approximate the true curve with approximately the same quality. In fact, we cannot qualitatively distinguish the best approximation from c
=
9 coefficients to c
=
33. This issue gives rise to the analysis presented on the next section.


## 3.2. Effective Dimension

We refer toeffective dimension as the number of parameters that can be properly retrieved from a noisy data set. In the context of this paper, effective dimension is related to the uncertainty quantification of parameters, that is, identifying the number of Fourier coefficients on the representation given by (17)–(20), required to approximate properly the boundary of the kite (see Figure 1) from synthetic far field measurements. As an aside comment we want to mention that effective dimension can be seen as an indicator for a parsimonious approximation regarding computational cost. However this feature is not used here due to the fact that effective dimension analysis is performed as a posterior procedure. In this section we have changed the notation for the distributions to be consistent with [14].We define a collection ofk models given by
(30)
π
u
∞
∣
Θ
n
n
(
u
∞
∣
θ
n
)
,
with n
∈
{
1,2
,
…
,
k
}. Θ
n is a vector of random variables. A realization of Θ
n is a vector θ
n defined in (20), and π
u
∞
∣
Θ
n
n is the corresponding posterior distribution for model n. For simplicity we use as index n to indicate the model corresponding to c
=
2
n
+
1 Fourier coefficients.We use the super-model approach; that is, a new model is defined(31)
π
u
∞
∣
Θ
,
M
(
u
∞
∣
θ
,
n
)
=
π
u
∞
∣
Θ
n
n
(
u
∞
∣
θ
n
)
,
where M
=
n is the indicator of the model and Θ
=
(
Θ
1
,
Θ
2
,
…
,
Θ
k
).The posterior probability for modelM
=
n of explaining the observed far field pattern is given by
(32)
π
(
M
=
n
∣
u
∞
)
∝
π
M
(
n
)
π
u
∞
n
(
u
∞
)
,
where π
M
(
n
) is the prior probability of model M
=
n and π
u
∞
n
(
u
∞
) is the corresponding normalizing constant (also known as the marginal likelihood) for model n. A straightforward way to perform model selection is by computing (32) for all the models and choosing the model with highest probability. Indeed, the difficulty regarding the computation of (32) is evaluating the normalizing constant π
u
∞
n
(
u
∞
). To perform our MCMC only the nonnormalized posterior is needed and therefore π
u
∞
n
(
u
∞
) was not previously calculated (this is the usual case when doing MCMC in a fixed number of parameters in a Bayesian application). Nevertheless, the MCMC simulations may be used to estimate π
u
∞
n
(
u
∞
) and provide information to explore the effective dimension. This is explained below.The ratio of two normalizing constants is defined as follows:(33)
r
=
π
u
∞
l
(
u
∞
)
π
u
∞
s
(
u
∞
)
,
where s
,
l
∈
{
1,2
,
…
,
k
}
,
s
≠
l. Considering a uniform prior distribution for the models, that is, π
M
(
n
)
=
1
∖
n, model s is better than model l when r
>
1, and they are indistinguishable if r
=
1 and l is better than s when r
<
1. Then, the best model has the highest normalizing constant value.A consistent estimator of the ratio is(34)
r
^
=
(
1
/
m
l
)
∑
i
=
1
m
l
‍
π
u
∞
∣
Θ
l
l
(
u
∞
∣
θ
l
,
i
)
π
Θ
l
l
(
θ
l
,
i
)
/
π
l
I
(
θ
l
,
i
)
(
1
/
m
s
)
∑
i
=
1
m
s
‍
π
u
∞
∣
Θ
s
s
(
u
∞
∣
θ
s
,
i
)
π
Θ
s
s
(
θ
s
,
i
)
/
π
s
I
(
θ
s
,
i
)
,
where π
l
I
,
π
s
I are completely known importance sampling densities for π
Θ
l
|
u
∞
l and π
Θ
s
|
u
∞
s, respectively, and the sets {
θ
l
,
1
,
θ
l
,
2
,
…
,
θ
l
,
m
l
} and {
θ
s
,
1
,
θ
s
,
2
,
…
,
θ
s
,
m
s
} are samples of the posterior distribution for models l and s, respectively, (see [14] for details).We observe from (34) that for a model s, the product
(35)
π
u
∞
∣
Θ
s
s
(
u
∞
∣
θ
s
,
i
)
π
Θ
s
s
(
θ
s
,
i
)
for
each
θ
s
,
i
∈
{
θ
s
,
1
,
θ
s
,
2
,
…
,
θ
s
,
m
s
}
,
is the corresponding likelihood times prior distribution. This product is in fact the energy (21) using a suitable constant E. Then the resulting sample from MCMC posterior simulations is used to compute (34), as a byproduct of the former with marginal computational cost added. Some details about approximating normalizing constant from MCMC simulations are discussed in [16].The importance sampling densityπ
s
I must be comparable to the product (35). We propose to numerically estimate the density π
s
I for every model using a Kernel density estimation (KDE) [17] of the posterior distribution using the MCMC sample. A rough estimate may be used, for the sole purpose of estimating r.The kernel density estimator has the form(36)
π
^
s
(
θ
s
*
∣
u
∞
)
=
1
m
s
h
m
s
∑
i
=
1
m
s
‍
K
s
(
θ
s
*
-
θ
s
,
i
h
m
s
i
)
,
where the kernel K is a bounded density on R
s and h
m
s is the bandwidth. For this estimator we define a Gaussian kernel
(37)
K
s
(
x
)
=
(
2
π
)
-
s
/
2
e
-
x
2
/
2
and a bandwidth
(38)
h
m
s
=
∏
i
=
1
m
s
‍
h
m
s
i
,
where h
m
s
i
=
1.06
σ
i
m
s
-
1
/
5 and σ
i is the standard deviation for each coefficient on the model s. We compute
(39)
c
s
=
(
1
m
s
)
∑
i
=
1
m
s
‍
π
u
∞
∣
Θ
s
s
(
u
∞
∣
θ
s
,
i
)
π
Θ
s
s
(
θ
s
,
i
)
π
s
I
(
θ
s
,
i
)
,
for each model s
∈
{
1,2
,
…
,
k
} then all the models can be compared together. Then effective dimension corresponds to the model with the highest normalizing constant, equivalently the lowest value of minus logarithm of the normalizing constant.We present in Figure7 the minus logarithm of the estimated normalizing constants (39) for all models (red line). As we see in this figure, the lowest value is obtained with n
=
11 coefficients. Also we present in the same Figure the distance (norm L
2 of the pointwise difference) between the radius of the kite and the radius of the simulated boundaries. The continuous line shows the difference with respect to the mean boundaries. The dashed line corresponds to the distances with respect to the map boundaries. The boxplot in the same figure corresponds to distances of the last 1000 simulations of the posterior. These last simulations are taken and subsampled with their corresponding IAT factor. These boundaries are in fact the probability region shown in Figure 6.Figure 7
Minus logarithm of estimated normalizing constants for Bayesian model comparisons. Best model refers to the highest normalizing constant (c
=
11 for this case). Continuous line shows the distance between the kite and the mean boundaries. Dashed line corresponds to the distance between the kite and the MAP. The boxplot corresponds to distances of the last 1000 simulations of the posterior to the kite. Red scale corresponds to minus logarithm of normalizing constants and blue scale is for the distances.Based on the Bayes factors, the best model corresponds ton
=
11 (we recall that c
=
2
n
+
1). Of note, there are MAP boundaries (e.g. c
∈
{
13
,
19
,
25
,
29
}) with lower mean squared error compared with the MAP for c
=
11. Also the mean boundaries from c
=
15 to c
=
33 have lower distances to the MAP. However as we observe on the boxplot, the mean of the error decreases from 3 to 11 Fourier coefficients and then oscillations appear, increasing spread and increasing error from 13 to 33 coefficients. Despite the low error of MAP and mean boundaries for models with 13 to 33 coefficients, due to the spread shown in the boxplot, the best model is certainly given by Bayes’ factor (i.e., for c
=
11). A qualitative analysis that agrees with these results is also discussed in Section 3.Fourier series is a classical representation for smooth periodic closed curves. Furthermore, we have used a well understood and high order numerical method for forward mapping evaluation. However, these two elements by themselves are not enough to have confidence regarding the solution of the inverse problem. The quality of the solution depends also on the number of coefficients of the Fourier series used for the representation.Our effective dimension discussion relies on MCMC methods to provide a quantitative strategy to determine the number of parameters that can be retrieved given a data set.


## 4. Conclusions

In this work we address the acoustic inverse scattering problem with a classical Fourier-based representation of the solution. We pose the inverse problem as a Bayesian inference problem and use the output of a MCMC method (namely, the t-walk) for our effective dimension results. For the corresponding direct problem we have used the classical layer potential approach, which was solved in a fast and reliable manner with a robust numerical method and parallel computing. Using Fourier series to represent solutions allows for a straightforward formulation that incorporates the smoothness of the solutions into the prior distribution. On the other hand, the finite Fourier representation is numerically unstable. Although other approaches to represent the scattering obstacle are applicable (e.g., wavelet basis which correspond to Besov priors), a fundamental question remains: How much information can be retrieved, within the representation, from a noisy data set?The main contribution of this paper is the effective dimension method, which is a quantitative method to estimate and quantify the uncertainty of the estimable parameters given a noisy data set. Given a parametric representation of the solution of the inverse problem, the effective dimension method is implemented via Bayesian model selection where the normalizing constant for each model is approximated using the MCMC output. Of note, the effective dimension method is applicable regardless of the parametric representation of the solution.



---

*Source: 427203-2014-09-07.xml*