SwitchNet: a neural network model for forward and inverse scattering problems

1. Introduction

In this paper, we study the forward and inverse scattering problems via the use of artificial neural networks (NNs). In order to simplify the discussion, we focus on the time-harmonic acoustic scattering in two dimensional space. The inhomogeneous media scattering problem with a fixed frequency $ω$ is modeled by the Helmholtz operator

(1)

L u := (- Δ - \frac{ω^{2}}{c^{2} (x)}) u,

where $c (x)$ is the velocity field. In many settings, there exists a known background velocity field $c_{0} (x)$ such that $c (x)$ is identical to $c_{0} (x)$ except in a compact domain $Ω$ . By introducing the scatterer $η (x)$ compactly supported in $Ω$

(2)

\frac{ω^{2}}{c (x)^{2}} = \frac{ω^{2}}{c_{0} (x)^{2}} + η (x),

we can equivalently work with $η (x)$ instead of $c (x)$ . Note that in this definition $η (x)$ scales quadratically with the frequency $ω$ . However, as $ω$ is assumed to be fixed throughout this paper, this scaling does not affect any discussion below.

In many real-world applications, $η (\cdot)$ is unknown. The task of the inverse problem is to recover $η (\cdot)$ based on some observation data $d (\cdot)$ . The observation data $d (\cdot)$ is often a quantity derived from the Green’s function $G = L^{- 1}$ of the Helmholtz operator $L$ and, therefore, it depends closely on $η (\cdot)$ . This paper is an exploratory attempt of constructing efficient approximations to the forward map $η \to d$ and the inverse map $d \to η$

using the modern tools from machine learning and artificial intelligence. Such approximations are highly useful for the numerical solutions of the scattering problems: an efficient map

η \to d

provides an alternative to expensive partial differential equation (PDE) solvers for the Helmholtz equation; an efficient map

d \to η

is more valuable as it allows us to solve the inverse problem of determining the scatterers from the scattering field, without going through the usual iterative process.

In the last several years, deep neural network has become the go-to method in computer vision, image processing, speech recognition and many other machine learning applications

[18, 28, 12, 9]. More recently, methods based on NN have also been applied to solving PDEs. Based on the way that the NN is used, these methods for solving the PDE can be roughly separated into two different categories. For the methods in the first category [16, 27, 3, 10, 14, 6], instead of specifying the solution space via the choice of basis (as in finite element method or Fourier spectral method), NN is used for representing the solution. Then an optimization problem, for example an variational formulation, is solved in order to obtain the parameters of the NN and hence the solution to the PDE. Similar to the use of an NN for regression and classification purposes, the methods in the second category such as [23, 11, 13, 8] use an NN to learn a map that goes from the coefficients in the PDE to the solution of the PDE. As in machine learning, the architecture design of an NN for solving PDE usually requires the incorporation of the knowledge from the PDE domain such that the NN architecture is able to capture the behavior of the solution process. Despite the abundance of the works in using the NN for solving PDE, none of the above mentioned methods have tried to obtain the solution to the wave equation.

This paper takes a deep learning approach to learn both the forward and inverse maps. For the Helmholtz operator (

1), we propose an NN architecture for determining the forward and inverse maps between the scatterer

η (\cdot)

and the observation data

d (\cdot)

generated from the scatterer. Although this task looks similar to the computer vision problems such as image segmentation, denoising, and super-resolution where the map between the two images has to be determined, the nature of the map in our problem is much more complicated. In many image processing tasks, the value of a pixel at the output generally only depends on a neighborhood of that pixel at the input layer. However, for the scattering problems, the input and output are often defined on different domains and, due to wave propagation, each location of the scatterer can influence every point of the scattered field. Therefore, the connectivity in the NN has to be wired in a non-local fashion, rendering typical NN with local connectivity insufficient. This leads to the development of the proposed

SwitchNet. The key idea is the inclusion of a novel low-complexity switch layer that sends information between all pairs of sites effectively, following the ideas from butterfly factorizations [21]. The same factorization was used earlier in the architecture proposed [19], but the network weights there are hardcoded and not trainable.

The paper is organized as followed. In Section 2, we discuss about some preliminary results concerning Helmholtz equation. In Section 3, we study the so called far field pattern of the scattering problem, where the sources and receivers can be regarded as placed at infinity. We propose SwitchNet to determine the maps between the far field scattering pattern and the scatterer. In Section 4, we turn to the setting of a seismic imaging problem. In this problem, the sources and receivers are at a finite distance, but yet well-separated from the scatterer.

2. Preliminary

The discussion of this paper shall focus on the two-dimensional case. Here, we summarize the mathematical tools and notations used in this paper. As mentioned above, the scatterer $η (x)$ is compactly supported in a domain $Ω$ , whose diameter is of $O (1)$ . For example, one can think of $Ω$ to be the unit square centered at the origin. In (1), the Helmholtz operator is defined on the whole space $R^{2}$ with the radiative (Sommerfeld) boundary condition [5] specified at infinity. Since the scatterer field $η (x)$ is localized in $Ω$ , it is convenient to truncate the computation domain to $Ω$ by imposing the perfectly matched layer [1] that approximates the radiative boundary condition.

In a typical numerical solution of the Helmholtz operator, $Ω$ is discretized by a Cartesian grid $X \subset Ω$ at the rate of a few points per wavelength. As a result, the number of grid points $N$ per dimension is proportional to the frequency $ω$ . We simply use ${x}_{x \in X}$ to denote the discretization points of this $N \times N$ grid $X$ . The Laplacian operator $- Δ$ in the Helmholtz operator is typically discretized with local numerical schemes, such as the finite difference method [17]. Via this discretization, we can consider the scatterer field $η$ , discretized at the points in $X$

, as a vector in

R^{N^{2}}

and the Helmholtz operator

L

as a matrix in

C^{N^{2} \times N^{2}}

Using the background velocity field $c_{0} (x)$ , we first introduce the background Helmholtz operator $L_{0} = - Δ - ω^{2} / c_{0}^{2}$ . With the help of $L_{0}$ , one can write $L$ in a perturbative way as

(3)

L = L_{0} - E, E = diag (η),

where $E$ is viewed as a perturbation. By introducing the background Green’s function

(4)

G_{0} := L_{0}^{- 1},

one can write down a formal expansion for the Green’s function $G = L^{- 1}$ of the $η$ -dependent Helmholtz operator $L$ :

(5)	$G$	$=$	$(L_{0} (I - G_{0} E))^{- 1}$
(6)		$\sim$	$(I + G_{0} E + G_{0} E G_{0} E + \dots) G_{0}$
(7)		$\sim$	$G_{0} + G_{0} E G_{0} + G_{0} E G_{0} E G_{0} + \dots$
(8)		$:=$	$G_{0} + G_{1} + G_{2} + \dots,$

which is valid when the scatterer field $η (x)$ is sufficiently small. The last line of the above equation serves as the definition of the successive terms of the expansion ( $G_{1}$ , $G_{2}$ , and so on). As $G_{0}$ can be computed from the knowledge of the background velocity field $c_{0} (x)$ , most data gathering processes (with appropriate post-processing) focus on the difference $G - G_{0} = G_{1} + G_{2} + \dots$ instead of $G$ itself.

A usual experimental setup consists of a set of sources $S$ and a set of receivers $R$ :

S = {s}_{s \in S}, R = {r}_{r \in R} .

The data gathering process usually involves three steps: (1) impose an external force or incoming wave field via some sources, (2) solve for the scattering field either computationally or physically, (3) gather the data with receivers at specific locations or directions. The second step is modeled by the difference of the Green’s function $G - G_{0}$ , as we mentioned above. As for the other steps, it is convenient at this point to model the first step with a source-dependent operator $Π_{S}$ and the third one with a receiver-dependent operator $Π_{R}$ . We shall see later how these operators are defined in more concrete settings. By putting these components together, one can set the observation data $d$ abstractly as

(9)

d = Π_{R} (G - G_{0}) Π_{S} = Π_{R} (G_{0} E G_{0} + G_{0} E G_{0} E G_{0} + \dots) Π_{S} = (Π_{R} G_{0}) (E + E G_{0} E + \dots) (G_{0} Π_{S}) .

In this paper, we focus on two scenarios: far field pattern and seismic imaging. We start with far field pattern first to motivate and introduce SwitchNet. We then move on to the seismic case by focusing on the main differences.

3. SwitchNet for far field pattern

3.1. Problem setup

In this section, we consider the problem of determining the map from the scatterer to the far field scattering pattern, along with its inverse map. Without loss of generality, we assume that the diameter of the domain $Ω$ is of $O (1)$ after appropriate rescaling. The background velocity $c_{0} (x)$ is assumed to be 1 since the far field pattern experiments are mostly performed in free space.

In this problem, both the sources and the receivers are indexed by a set of unit directions in $S^{1}$ . The source associated with a unit direction $s \in S \subset S^{1}$ is an incoming plane waves $e^{i ω s \cdot x}$ pointing at direction $s$ . It is well known that the scattered wave field, denoted by $u_{s} (x)$ , at a large distance takes the following form [5]

u_{s} (x) = \frac{e^{i ω | x |}}{\sqrt{| x |}} (u_{s}^{\infty} (\frac{x}{| x |}) + o (1)),

where the function $u_{s}^{\infty} (\cdot)$ is defined on the unit circle $S^{1}$ . The receiver at direction $r \in R \subset S^{1}$ simply records the quantity $u_{s}^{\infty} (r)$ for each $s$ . The set of observation data $d$ is then defined to be

d (r s) = u_{s}^{\infty} (r) .

Figure 1 provides an illustration of this experimental setup. Henceforth, we assume that both $R$ and $S$

are chosen to be a set of uniformly distributed directions on

S^{1}

. Their size, denoted by

M

, typically scales linearly with frequency

ω

Figure 1. Illustration of the incoming and outgoing waves for a far field pattern problem. The scatterer $η (x)$ is compactly supported in the domain $Ω$ . The incoming plane wave points at direction $s$ . The far field pattern is sampled at each receiver direction $r$ .

This data gathering process can be put into the framework of (9). First, one can think of the source prescription as a limiting process that produces in the limit the incoming wave $e^{i ω s \cdot x}$ . The source can be considered to be located at the point $- s ρ$ for the direction $s \in S^{1}$ with the distance $ρ \in R^{+}$ going to infinity. In order to compensate the geometric spreading of the wave field and also the phase shift, the source magnitude is assumed to scale like $\sqrt{ρ} e^{- i ω ρ}$ as $ρ$ goes to infinity. Under this setup, we have

(10)		$lim ρ \to \infty (G_{0} Π_{S}) (x, s)$
(11)	$=$	$lim ρ \to \infty (1 / \sqrt{ρ}) e^{i ω \| x - (- s ρ) \|} \sqrt{ρ} e^{- i ω ρ}$
(12)	$=$	$lim ρ \to \infty (1 / \sqrt{ρ}) e^{i ω (ρ + s \cdot x)} \sqrt{ρ} e^{- i ω ρ}$
(13)	$=$	$e^{i ω s \cdot x} .$

Similarly, one can also regard the receiver prescription as a limiting process as well. The receiver is located at point $r ρ^{'}$ for a fixed unit direction $r \in S^{1}$ with $ρ^{'} \in R^{+}$ going to infinity. Again in order to to compensate the geometric spreading and the phase shift, one scales the received signal with $\sqrt{ρ^{'}} e^{- i ω ρ^{'}}$ . As a result, we have

(14)		$lim ρ^{'} \to \infty (Π_{R} G_{0}) (r, x)$
(15)	$=$	$lim ρ^{'} \to \infty (1 / \sqrt{ρ^{'}}) e^{i ω \| r ρ^{'} - x \|} \sqrt{ρ^{'}} e^{- i ω ρ^{'}}$
(16)	$=$	$lim ρ^{'} \to \infty (1 / \sqrt{ρ^{'}}) e^{i ω (ρ^{'} - r \cdot x)} \sqrt{ρ^{'}} e^{- i ω ρ^{'}}$
(17)	$=$	$e^{- i ω r \cdot x} .$

In this limiting setting, one redefine the observation data as

(18)

d = lim ρ, ρ^{'} \to \infty (Π_{R} G_{0}) (E + E G_{0} E + \dots) (G_{0} Π_{S}) .

Now taking the two limits (10) and (14) under consideration, one arrives at the following representation of the observation data $d (r, s)$ for $r \in R$ and $s \in S$

(19)

d (r, s) = \sum x \in X \sum y \in X e^{- i ω r \cdot x} (E + E G_{0} E + \dots) (x, y) e^{i ω s \cdot y} .

3.2. Low-rank property

The intuition behind the proposed NN architecture comes from examining (19) when $E$ (or $η$ ) is small. In such a situation, we simply retain the term that is linear in $E$ . Using the fact that $E = diag (η)$ , (19) becomes

d (r, s) \approx \sum x \in X e^{i ω (s - r) \cdot x} η (x)

for $r \in R \subset S^{1}$ and $s \in S \subset S^{1}$ . This linear map takes $η (x)$ defined on a Cartesian grid $X \subset Ω$ to $d (r, s)$ defined on yet another Cartesian grid $R \times S \subset S^{1} \times S^{1}$ . Recalling that both $R$ and $S$ are of size $M$ and working with a vectorized $d \in C^{M^{2}}$ , we can write the above equation compactly as

(20)

d \approx A η,

where the element of the matrix $A \in C^{M^{2} \times N^{2}}$ at $(r, s) \in R \times S$ and $x \in X$ is given by

(21)

A (r s, x) = exp (i ω (s - r) \cdot x) .

The following theorem concerning the matrix $A$ plays a key role in the design of our NN. Let us first partition $Ω$ uniformly into $\sqrt{P_{X}} \times \sqrt{P_{X}}$ Cartesian squares of side-length equal to $1 / \sqrt{ω}$ . Here we assume that $\sqrt{P_{X}}$ is an integer. Note that, since the diameter of $Ω$ is of $O (1)$ , $\sqrt{P_{X}} \approx \sqrt{ω}$ . This naturally partitions the set of grid points $X$ into $P_{X}$ subgroups depending on which square each point belongs to. We shall denote these subgroups by $X_{0}, \dots, X_{P_{X} - 1}$ . Similarly, we also partition $S^{1} \times S^{1}$ uniformly (in the angular parameterization) into $\sqrt{P_{D}} \times \sqrt{P_{D}}$ squares $D_{0}, \dots, D_{P_{D} - 1}$ of side-length equal to $1 / \sqrt{ω}$ . $\sqrt{P_{D}}$ is also assumed to be an integer and obviously $\sqrt{P_{D}} \approx \sqrt{ω}$ . This further partitions the set $R \times S$ into $P_{D}$ subgroups depending on which square they belong to. We shall denote these subgroups by $D_{0}, \dots, D_{P_{D} - 1}$ . Fig. 2 illustrates the partition for $\sqrt{P_{X}} = \sqrt{P_{D}} = 4$ .

Theorem 1.

For any $D_{i}$ and $X_{j}$ , the submatrix

(22)

A_{i j} := [A (r s, x)]_{(r, s) \in D_{i}, x \in X_{j}}

is numerically low-rank.

Proof.

The proof of this theorem follows the same line of argument in [2, 29, 20] and below we outline the key idea. Denote the center of $D_{i}$ by $(r_{i}, s_{i})$ and the center of $X_{j}$ by $x_{j}$ . For each $(r, s) \in D_{i}$ and $x \in X_{j}$ , we write

(23)

Note that for fixed $D_{i}$ and $X_{j}$ each of the last three terms is either a constant or depends only on $x$ or $(r, s)$ . As a result, $exp (i ω (s - r) \cdot x)$ is numerically low-rank if and only if the first term $exp (i ω ((s - r) - (s_{i} - r_{i})) \cdot (x - x_{j}))$ is so. Such a low-rank property can be derived from the conditions concerning the side-lengths of $D_{i}$ and $X_{j}$ . More precisely, since $(r, s)$ resides in $D_{i}$ with center $(r_{i}, s_{i})$ , then

(24)

| (s - r) - (s_{i} - r_{i}) | \leq \frac{1}{\sqrt{ω}} .

Similarly as $x$ resides in $X_{j}$ with center $x_{j}$ , then

(25)

| x - x_{j} | \leq \frac{1}{\sqrt{ω}} .

Multiplying these two estimates results in the estimate

(26)

ω | ((s - r) - (s_{i} - r_{i})) \cdot (x - x_{j})) | \leq 1

for the phase of $exp (i ω ((s - r) - (s_{i} - r_{i})) \cdot (x - x_{j}))$ . Therefore,

(27)

exp (i ω ((s - r) - (s_{i} - r_{i})) \cdot (x - x_{j}))

for $(r, s) \in D_{i}$ or $x \in X_{j}$

is non-oscillatory and hence can be approximated effectively by applying, for example, Chebyshev interpolation in both the

(r, s)

and

x

variables. Since the degree of the Chebyshev polynomials only increases poly-logarithmically with respect to the desired accuracy,

exp (i ω ((s - r) - (s_{i} - r_{i})) \cdot (x - x_{j}))

is numerically low-rank by construction. This proves that the submatrix

A_{i j}

defined in (22) is also numerically low-rank. ∎

3.3. Matrix factorization

In this subsection, we show that Theorem 1 guarantees a low-complexity factorization of the matrix $A$ . Let the row and column indices of $A \in C^{M^{2} \times N^{2}}$ be partitioned into index sets ${D_{i}}_{i = 0}^{P_{D} - 1}$ and ${X_{j}}_{j = 0}^{P_{X} - 1}$ , respectively, as in Theorem 1. To simplify the presentation, we assume $P_{X} = P_{D} = P$ , $| X_{0} | = \dots | X_{P - 1} | = N^{2} / P$ , and $| D_{0} | = \dots = | D_{P - 1} | = M^{2} / P$ .

Since the submatrix

A_{i j} := [A (r s, x)]_{r s \in D_{i}, x \in X_{j}}

is numerically low-rank, assume that

(28)

A_{i j} \approx U_{i j} V_{i j}^{*},

where $U_{i j} \in C^{M^{2} / P \times t}$ and $V_{i j} \in C^{N^{2} / P \times t}$ . Here $t$ can be taken to be the maximum of the numerical ranks of all submatrices $A_{i j}$ . Theorem 1 implies that $t$ is a small constant.

By applying (28) to each block $A_{i j}$ , $A$ can be approximated by

(29)

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} U_{00} V_{00}^{*} & U_{01} V_{01}^{*} & \dots & U_{0 (P - 1)} V_{0 (P - 1)}^{*} U_{10} V_{10}^{*} & U_{11} V_{11}^{*} & \dots & U_{1 (P - 1)} V_{1 (P - 1)}^{*} ⋮ & ⋱ & ⋮ U_{(P - 1) 0} V_{(P - 1) 0}^{*} & U_{(P - 1) 1} V_{(P - 1) 1}^{*} & \dots & U_{(P - 1) (P - 1)} V_{(P - 1) (P - 1)}^{*} \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .

The next step is to write (29) into a factorized form. First, introduce $U_{i}$ and $V_{j}$

(30)

U_{i} = [\begin{matrix} U_{i 0}, U_{i 1}, \dots, U_{i (P - 1)} \end{matrix}] \in C^{M^{2} / P \times t P}, V_{j} = [\begin{matrix} V_{0 j}, V_{1 j}, \dots, V_{(P - 1) j} \end{matrix}] \in C^{N^{2} / P \times t P},

and define in addition

(31)

In addition, introduce

(32)

Σ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} Σ_{00} & Σ_{01} & \dots & Σ_{0 (P - 1)} Σ_{10} & Σ_{11} & \dots & Σ_{1 (P - 1)} ⋮ & ⋱ & ⋮ Σ_{(P - 1) 0} & Σ_{(P - 1) 1} & \dots & Σ_{(P - 1) (P - 1)} \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ \in C^{P^{2} t \times P^{2} t},

where the submatrix $Σ_{i j} \in C^{P t \times P t}$ itself is a $P \times P$ block matrix with blocks of size $t \times t$ . $Σ_{i j}$

is defined to be zero everywhere except being the identity matrix at the

(j, i)

-th

t \times t

block. In order to help understand the NN architecture discussed below sections, it is imperative to understand the meaning of

Σ

. Let us assume for simplicity that

t = 1

. Then for an arbitrary vector

z \in C^{P^{2}}

Σ z

essentially performs a “switch” that shuffles

z

as follows

(33)

(Σ z) (j P + i) = z (i P + j), i, j = 0, \dots, P - 1.

With the above definitions for $U$ , $V$ , and $Σ$ , the approximation in (29) can be written compactly as

(34)

A \approx U Σ V^{*} .

Notice that although $A$ has $M^{2} \times N^{2}$ entries, using the factorization (34), $A$ can be stored using $t P (M^{2} + P + N^{2})$ entries. In this paper, $P \approx max (M, N)$ and $M$ and $N$ are typically on the same order. Therefore, instead of $O (N^{4})$ , one only needs $O (N^{3})$ entries to parameterize the map $A$ approximately using (34). Such a factorization is also used in [21] for the compression of Fourier integral operators.

We would like to comment on another property that may lead to further reduction in the parameters used for approximating $A$ . Let us focus on any two submatrices $A_{i j}$ and $A_{i k}$ of $A$ . For two regions $X_{j}$ and $X_{k}$ , where the center of $X_{j}$ and $X_{k}$ are $x_{j}$ and $x_{k}$ respectively, $X_{k} = X_{j} + (x_{k} - x_{j})$ . Let $(r, s) \in D_{i}$ . For $x \in X_{j}$ and $x^{'} = x + (x_{k} - x_{j}) \in X_{k}$ , we have

(35)		$exp (i ω (s - r) \cdot x)$	$=$	$g_{1} (r, s) h ((r, s), x),$
(36)		$exp (i ω (s - r) \cdot x^{'})$	$=$	$g_{2} (r, s) h ((r, s), x),$

where

(37)

g_{1} (r, s) = exp (i ω (s - r) \cdot x_{j}), g_{2} (r, s) = exp (i ω (s - r) \cdot x_{k}), h ((r, s), x) = exp (i ω (s - r) \cdot (x - x_{j})) .

Therefore, the low-rank factorizations of $A_{i j}$ and $A_{i k}$ are solely determined by the factorization of $h (r s, x)$ . This implies that it is possible to construct low-rank factorizations for $A_{i j}$ and $A_{i k}$ :

(38)

A_{i j} \approx U_{i j} V_{i j}^{*}, A_{i k} \approx U_{i k} V_{i k}^{*},

such that $V_{i j}^{*} = V_{i k}^{*}$ . Since this is true for all possible $j, k$ , one can pick low-rank factorizations so that $V_{0} = V_{1} = \dots = V_{P - 1}$ .

As a final remark in this section, this low complexity factorization (34) for $A$ can be easily converted to one for $A^{*}$ since

(39)

A^{*} \approx V Σ^{*} U^{*},

where $U, Σ, V$ are provided in (30), (31), and (32).

3.4. Neural networks

Based on the low-rank property of $A$ in Section 3.2 and its low-complexity factorization in Section 3.3, we propose new NN architectures for representing the inverse map $d \to η$ and the forward map $η \to d$ .

3.4.1. NN for the inverse map $d \to η$ .

As pointed out earlier, $d \approx A η$ when $η$ is sufficiently small. The usual filtered back-projection algorithm [25] solves the inverse problem $d \to η$ via

(40)

η \approx (A^{*} A + ϵ I)^{- 1} A^{*} d,

where $ϵ$ is the regularization parameter. In the far field pattern problem, $(A^{*} A + ϵ I)^{- 1}$ can be understood as a deconvolution operator. To see this, a direct calculation reveals that

(41)

(A^{*} A) (x, y) = \sum r s \in R \times S e^{i ω (s - r) \cdot y} e^{- i ω (s - r) x} = \sum r s \in R \times S e^{- i ω (s - r) (x - y)}

for $x, y \in X$ . (41) shows that $A^{*} A$ is a translation-invariant convolution operator. Therefore, the operator $(A^{*} A + ϵ I)^{- 1}$ , as a regularized inverse of $A^{*} A$ , simply performs a deconvolution. In summary, the above discussion shows that in order to obtain $η$ from the scattering pattern $d$ in the regime of small $η$ , one simply needs to apply sequentially to $d$

the operator $A^{*}$ ,
a translation-invariant filter that performs the deconvolution $(A^{*} A + ϵ I)^{- 1}$ .

Although these two steps might be sufficient when $η$ is small, a nonlinear solution is needed when $η$ is not so. For this purpose, we propose a nonlinear neural network SwitchNet for the inverse map. There are two key ingredients in the design of SwitchNet.

The first key step is the inclusion of a Switch layer that sends local information globally, as depicted in Figure 4. The structure of the Switch layer is designed to mimic the matrix-vector multiplication of the operator $A^{*} \approx V Σ^{*} U^{*}$ in (29). However unlike the fixed coefficients in (29), as an NN layer, the Switch layer allows for tunable coefficients and learns the right values for the coefficients from the training data. This gives the architecture a great deal of flexibility.
The second key step is to replace the linear deconvolution in the back-projection algorithm with a few convolution (Conv) layers. This enriches the architecture with nonlinear capabilities when approximating the nonlinear inverse map.

t, P_{D}, P_{X}, N, w, α, L, d \in C^{M \times M}

η \in C^{N \times N}

d_{1} \leftarrow Vect [P_{D}] (d)

d_{2} \leftarrow Switch [t, P_{D}, P_{X}, N^{2}] (d_{1})

e_{0} \leftarrow Square [P_{X}] (d_{2})

6:for

ℓ

from

0

L - 1

e_{ℓ + 1} \leftarrow Conv [w, α] (e_{ℓ})

8:end for

η \leftarrow Conv [w, 1] (e_{L})

10:return

η

Algorithm 1 SwitchNet for the inverse map

d \to η

of far field pattern.

The pseudo-code for SwitchNet is summarized in Algorithm 1. The input $d$ is a $C^{M \times M}$ matrix, while the output $η$ is a $C^{N \times N}$ matrix. The first three steps of Algorithm 1 mimics the application of the operator $A^{*} \approx V Σ^{*} U^{*}$ . The Switch layer does most of the work, while the Vect and Square layers are simply operations that reshape the input and output data to the correct matrix form at the beginning and the end of the Switch layer. In particular, Vect groups the entries of the 2D field $d$ according to squares defined by the partition $D_{0}, \dots, D_{P_{D} - 1}$ and Square does the opposite. The remaining lines of Algorithm 1 simply apply the Conv layers with window size $w$ and channel number $α$ .

These basic building blocks of SwitchNet are detailed in the following subsection. We also take the opportunity to include the details of the pointwise multiplication PM layer that will be used in later on.

3.4.2. Layers for SwitchNet

In this section we provide the details for the layers that are used in SwitchNet. Henceforth, we assume that the entries of a tensor is enumerated in the Python convention, i.e., going through the dimensions from the last one to the first. One operation that will be used often is a

reshape, in which a tensor is changed to a different shape with the same number of entries and with the enumeration order of the entries kept unchanged.

Figure 3. An illustration of the Vect and Square layers. The detail descriptions of the layers are provided in Section 3.4.2. For the purpose of illustration we let $P = 4$ . The Vect layer vectorize a $4 \times 4$ matrix on the left hand side according to the partitioning by $2 \times 2$ blocks, to give the size 16 vector on the right hand side. The Square layer is simply the adjoint map of the Vect layer.

Vectorize layer. $z_{O} = Vect [P] (z_{I})$ with input $z_{I} \in C^{n \times n}$ . Henceforth we assume that $\sqrt{P}$ is an integer and $\sqrt{P}$ divides $n$ . This operation partitions $z_{I}$ into $\sqrt{P} \times \sqrt{P}$ square sub-blocks of equal size. Then each sub-block is vectorized, and the vectorized sub-blocks are stacked together as a vector in $C^{n^{2}}$ . Intuitively, these operations cluster the nearby entries in a sub-block together. The details of the Vect layer are given in the following:

Reshape the $z_{I}$ to a $\sqrt{P} \times \frac{n}{\sqrt{P}} \times \sqrt{P} \times \frac{n}{\sqrt{P}}$ tensor.
Swap the second and the third dimensions to get a $\sqrt{P} \times \sqrt{P} \times \frac{n}{\sqrt{P}} \times \frac{n}{\sqrt{P}}$ tensor.
Reshape the result to an $n^{2}$ vector and set it to $z_{O}$ .

Square layer. $z_{O} = Square [P] (z_{I})$ with input $z_{I} \in C^{n^{2}}$ , where $\sqrt{P}$ is an integer. The output is $z_{O} \in R^{n \times n}$ . Essentially as the adjoint operator of the Vect layer, this layer fills up each square sub-block of the matrix $z_{O}$ with a segment of entries in $z_{I}$ . The details are given as followed:

Reshape the $z_{I}$ to a $\sqrt{P} \times \sqrt{P} \times \frac{n}{\sqrt{P}} \times \frac{n}{\sqrt{P}}$ tensor.
Swap the second and the third dimensions to get a $\sqrt{P} \times \frac{n}{\sqrt{P}} \times \sqrt{P} \times \frac{n}{\sqrt{P}}$ tensor.
Reshape the result to an $n \times n$ matrix and set it to $z_{O}$ .

Figure 4. An illustration of the Switch layer where the detail description of it is provided in Section 3.4.2. For the purpose of illustration we let $n_{I} = n_{O} = 20, P_{I} = P_{O} = 4$ .

Switch layer. $z_{O} = Switch [t, P_{I}, P_{O}, n_{O}] (z_{I})$ with input $z_{I} \in C^{n_{I}}$ . It is assumed that $n_{I}$ and $n_{O}$ are integer multiples of $P_{I}$ and $P_{O}$ , respectively. This layer consists the following steps.

Apply $U^{T}$ to $z_{I}$ :

$U^{T} = ⎡ ⎢ ⎢ ⎢ ⎣ \begin{matrix} U_{0}^{T} ⋱ U_{P_{I} - 1}^{T} \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎦, U_{0}^{T}, \dots, U_{P_{I} - 1}^{T} \in C^{t P_{O} \times \frac{n_{I}}{P_{I}}}$
Reshape $z_{1}$ to be a $C^{P_{O} \times P_{I} \times t}$ tensor. Here we follow the Python convention of going through the dimensions from the last one to the first one. Then a permutation is applied to swap the first two dimensions to obtain a tensor of size $C^{P_{O} \times P_{I} \times t}$ . Finally, the result is reshaped to a vector $z_{2} \in C^{P_{I} P_{O} t}$ again going through the dimensions from the last to the first.
Apply $V$ to $z_{2}$ :

$z_{O} = V z_{2} \in R^{n_{O}},$

$V = ⎡ ⎢ ⎢ ⎢ ⎣ \begin{matrix} V_{0} ⋱ V_{P_{O} - 1} \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎦, V_{0}, \dots, V_{P_{O} - 1} \in C^{\frac{n_{O}}{P_{O}} \times t P_{I}} .$

Here the non-zero entries of $U, V$ are the trainable parameters. The Switch layer is illustrated in Figure 4.

Convolution layer. $z_{O} = Conv [w, c_{O}] (z_{I})$ with input $z_{I} = C^{n \times n \times c_{I}}$ . Here $c_{I}, c_{O}$ denote the input and output channel numbers and $w$

denotes the window size. In this paper we only use the convolution layer with stride 1 and with zero-padding:

(43)		$z_{O} (k_{1}, k_{2}, k_{3})$	$=$	$ReLU (min (n - 1, k_{1} + \frac{w - 1}{2}) \sum l_{1} = max (0, k_{1} - \frac{w - 1}{2}) min (n - 1, k_{2} + \frac{w - 1}{2}) \sum l_{2} = max (0, k_{2} - \frac{w - 1}{2}) c_{I} - 1 \sum l_{3} = 0$
(43)				$W (l_{1} - k_{1} + \frac{w - 1}{2}, l_{2} - k_{2} + \frac{w - 1}{2}, l_{3}, k_{3}) z_{I} (l_{1}, l_{2}, l_{3}) + b (k_{3}))$

with $k_{1}, k_{2} = 0, \dots, n - 1, k_{3} = 0, \dots, c_{O} - 1$ . Here $ReLU (x) = max (0, x)$ and $w$

is assumed to be odd in the presentation. Both

W \in C^{w \times w \times c_{I} \times c_{O}}

and

b \in C^{c_{O}}

are the trainable parameters.

Pointwise multiplication layer. $z_{O} = PM (z_{I})$ with input $z_{I} \in C^{n \times n \times c_{I}}$ . It is defined as

(44)

z_{O} (k_{1}, k_{2}) = W (k_{1}, k_{2}) z_{I} (k_{1}, k_{2}) + b (k_{1}, k_{2}),

$k_{1}, k_{2} = 0, \dots, n - 1$ . Both $W \in C^{n \times n}$ and $b \in C^{n \times n}$ are trainable parameters.

We remark that, among these layers, the Switch layer has the most parameters. If the input and output to the Switch layer both have size $n \times n$ , the number of parameter is $2 t P n^{2}$ where $P$ is the number of squares that partition the input field and $t$ is the rank of the low-rank approximation.

3.4.3. NN for the forward map $η \to d$ .

We move on to discuss the parameterization of the forward map $η \to d$ . The proposal is based on the simple observation that the inverse of the inverse map is the forward map.

More precisely, we simply reverse the architecture of the inverse map proposed in Algorithm 1. This results in an NN presented Algorithm 2. The basic architecture of this NN involves applying a few layers of Conv first, then followed by a Switch layer that mimics $A \approx U Σ V^{*}$ .

t, P_{D}, P_{X}, M, w, α, L, η \in R^{N \times N}

d \in R^{M \times M \times 2}

η_{0} \leftarrow η

4:for

ℓ

from

0

L - 1

η_{ℓ + 1} \leftarrow Conv [w, α] (η_{ℓ})

6:end for

d_{1} \leftarrow Conv [w, 1] (η_{L})

d_{2} \leftarrow Vect [P_{X}] (d_{1})

d_{3} \leftarrow Switch [t, P_{X}, P_{D}, M^{2}] (d_{2})

10:

d \leftarrow Square [P_{D}] (d_{3})

11:return

d

Algorithm 2 SwitchNet for the forward map

η \to d

of far field pattern.

We would also like to mention yet another possibility to parameterize the forward map $η \to d$

, via a recurrent neural network

[24]. Let

(45)

E_{e f f} = E + E G_{0} E + E G_{0} E G_{0} E + \dots =: E_{1} + E_{2} + E_{3} + \dots .

One can leverage the following recursion

(46)

E_{k + 1} = E G_{0} E_{k} k = 1, 2, \dots K

to approximate $E_{e f f}$ by treating each $E_{k}$ as an $N^{2} \times N^{2}$ image and using a recurrent neural network. At the $k$ -th level of the recurrent neural network, it takes $E_{k}$ and $E$ as inputs and outputs $E_{k + 1}$ . More specifically, in order to go from $E_{k}$ to $E_{k + 1}$ , one first apply $G_{0}$ to each column of the image $E_{k}$ , then each row of the image is reweighted by the diagonal matrix $E$ . Stopping at the $K$ -th level for a sufficiently large $K$ , $E_{e f f}$ can be approximated by

(47)

E_{e f f} \approx K + 1 \sum i = 1 E_{i} .

Once holding such an approximation to $E_{e f f}$ , we plug it into (19)

(48)

d (r, s) = \sum x \in X \sum y \in X e^{i ω r \cdot x} E_{e f f} (x, y) e^{- i ω s \cdot y} = \sum x \in X e^{i ω r \cdot x} (\sum y \in X E_{e f f} (x, y) e^{- i ω s \cdot y}) .

This shows that the map from $E_{e f f}$ to $d$ can be realized by applying a matrix product to $E_{e f f}$ first on the $y$ -dimension, then on the $x$ -dimension. If we view applying the Green’s function $G_{0}$ as applying a convolution layer in an NN, the above discussion shows that the forward map can be obtained by first applying a recurrent NN followed by a convolutional NN. The main drawback of this approach is the large memory requirement (i.e., $N^{2} \times N^{2}$ ) to store each individual $E_{k}$ . In addition, the use of a recurrent NN may lead to difficulty in training [26] due to the issue of exploding or vanishing gradient. Moreover, since the weights for parameterizing $G_{0}$ are shared over multiple layers in the recurrent NN, one might not be able to efficiently use back-propagation, which may lead to a longer training time. These are the main reasons why we decided to adopt the approach in Algorithm 2.

3.5. Numerical results

In this section, we present numerical results of SwitchNet for far field pattern at a frequency $ω \approx 60$ . The scatterer field $η (x)$ supported in $Ω = [- 0.5, 0.5]^{2}$ is assumed to be a mixture of Gaussians

(49)

n_{s} \sum i = 1 β exp (- \frac{| x - c_{i} |^{2}}{2 σ^{2}})

where $β = 0.2$ and $σ = 0.015$ . When preparing the training and testing examples, the centers ${c_{i}}_{i = 1}^{n_{s}}$ of the Gaussians are chosen to be uniformly distributed within $Ω$ . The number $n_{s}$ of the Gaussians in the mixture is set to vary between $2$ and $4$ . In the numerical experiments, the domain $Ω = [- 0.5, 0.5]^{2}$ is discretized by an $80 \times 80$ Cartesian grid $X$ . To discretize the source and receiver directions, we set both $R$ and $S$ to be a set of 80 equally spaced unit directions on $S^{1}$ . Therefore in this example, $N = M = 80$ .

In Algorithm 1, the parameters are specified as $t = 3$ (rank of the low-rank approximation), $P_{X} = 8^{2}$ , $P_{D} = 4^{2}$ , $w = 10$ (window size of the convolution layers), $α = 18$ (channel number of the convolution layers), and $L = 3$ (number of convolution layers), resulting 3100K number of parameters. The parameters for Algorithm 2 are chosen to be $t = 4$ , $P_{X} = 8^{2}$ , $P_{D} = 4^{2}$ , $w = 10$ , $α = 24$ , and $L = 3$ , with a total of 4200K parameters. Note that for both algorithms the number of parameters is significantly less than the one of a fully connected NN, which has at least $80^{4} = 40960$ K parameters.

SwitchNet is trained with the ADAM optimizer [15]

in Keras

[4]

with a step size of 0.002 and a mini-batch size of size 200. The optimization is run for 2500 epochs. Both the training and testing data sets are obtained by numerically solving the forward scattering problem with an accurate finite difference scheme with a perfectly matched layer. In the experiment, 12.5K pairs of

(η, d)

are used for training, and another 12.5K pairs are reserved for testing. The errors are reported using the mean relative errors

(50)

\frac{1}{N_{t e s t}} N_{t e s t} \sum i = 1 \frac{∥ d_{i}^{N N} - d_{i} ∥_{F}}{∥ d_{i} ∥_{F}}, \frac{1}{N_{t e s t}} N_{t e s t} \sum i = 1 \frac{∥ η_{i}^{N N} - η_{i} ∥_{F}}{∥ η_{i} ∥_{F}},

where $d_{i}^{N N}$ and $d_{i}$ denote the predicted and ground truth scattering patterns respectively for the $i$ -th testing sample, and $η_{i}^{N N}$ and $η_{i}$ denote the predicted and ground truth scatterer field respectively. Here $∥ \cdot ∥_{F}$ is the Frobenius norm.

Table 1 summarizes the test errors for Gaussian mixtures with different choices of $n_{s}$ . For the purpose of illustration, we show the predicted $d$ and $η$ by SwitchNet along with the ground truth in Figure 5 for one typical test sample.

$n_{s}$	Forward map	Inverse map
2	9.4e-03	1.2e-02
3	4.0e-02	1.4e-02
4	4.8e-02	2.4e-02

Table 1. Prediction error of SwitchNet for the maps

η \to d

and

d \to η

for far field pattern.

(a) Ground truth scattering pattern $d$ .

4. SwitchNet for seismic imaging

4.1. Problem setup

This section considers a two-dimensional model problem for seismic imaging. The scatterer $η (x)$ is again assumed to be supported in a domain $Ω$ with an $O (1)$ diameter, after appropriate rescaling. $Ω$ is discretized with a Cartesian grid $X = {x}_{x \in X}$ at the rate of at least a few point per wavelength. Compared to the source and receiver configurations in Section 3.1, the experiment setup here is simpler. One can regard both $S = {s}_{s \in S}$ and $R = {r}_{r \in R}$ to be equal to a set of uniformly sampled points along a horizontal line near the top surface of the domain. The support of $η$ is at a certain distance below the top surface so that it is well-separated from the sources and the receivers (see Figure 6 for an illustration of this configuration).

Figure 6. Illustration of a simple seismic imaging setting. The sources ( $S$ ) and receivers ( $R$ ) are located near the surface level (top) of the domain $Ω$ . The scatterer field $η (x)$ is assumed to be well-separated from the sources and the receivers.

The source and receiver operators in (9) take a particularly simple form. For the sources, the operator $(G_{0} Π_{S})$ is simply given by sampling:

(G_{0} Π_{S}) (x, s) = G_{0} (x, s) .

Similarly for the receivers, the operator $(Π_{R}^{T} G_{0})$ is given by

(Π_{R} G_{0}) (r, y) = G_{0} (r, y) .

After plugging these two formulas back into (9), one arrives at the following representation of the observation data $d (r, s)$ for $r \in R$ and $s \in S$

d (r, s) = \sum x \in X \sum y \in X G_{0} (r, x) (E + E G_{0} E + \dots) (x, y) G_{0} (y, s) .

4.2. Low-rank property

Following the approach taken in Section 3.2, we start with the linear approximation under the assumption that $η (x)$ is weak. Since $E = diag (η)$ , the first order approximation is

(51)

d (r,

SwitchNet: a neural network model for forward and inverse scattering problems

Solving Inverse Wave Scattering with Deep Learning

Sampling type method combined with deep learning for inverse scattering with one incident wave

Solving the Wide-band Inverse Scattering Problem via Equivariant Neural Networks

Reduced Order Model Approach to Inverse Scattering

Fast Fitting of Reflectivity Data of Growing Thin Films Using Neural Networks

Extracting Electron Scattering Cross Sections from Swarm Data using Deep Neural Networks

Isogeometric multilevel quadrature for forward and inverse random acoustic scattering

1. Introduction

2. Preliminary

3. SwitchNet for far field pattern

3.1. Problem setup

3.2. Low-rank property

Theorem 1.

Proof.

3.3. Matrix factorization

3.4. Neural networks

3.4.1. NN for the inverse map $d \to η$ .

3.4.2. Layers for SwitchNet

3.4.3. NN for the forward map $η \to d$ .

3.5. Numerical results

4. SwitchNet for seismic imaging

4.1. Problem setup

4.2. Low-rank property


	$U^{T} = ⎡ ⎢ ⎢ ⎢ ⎣ \begin{matrix} U_{0}^{T} ⋱ U_{P_{I} - 1}^{T} \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎦, U_{0}^{T}, \dots, U_{P_{I} - 1}^{T} \in C^{t P_{O} \times \frac{n_{I}}{P_{I}}}$

	$z_{O} = V z_{2} \in R^{n_{O}},$
	$V = ⎡ ⎢ ⎢ ⎢ ⎣ \begin{matrix} V_{0} ⋱ V_{P_{O} - 1} \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎦, V_{0}, \dots, V_{P_{O} - 1} \in C^{\frac{n_{O}}{P_{O}} \times t P_{I}} .$

SwitchNet: a neural network model for forward and inverse scattering problems

Related Research

Solving Inverse Wave Scattering with Deep Learning

Sampling type method combined with deep learning for inverse scattering with one incident wave

Solving the Wide-band Inverse Scattering Problem via Equivariant Neural Networks

Reduced Order Model Approach to Inverse Scattering

Fast Fitting of Reflectivity Data of Growing Thin Films Using Neural Networks

Extracting Electron Scattering Cross Sections from Swarm Data using Deep Neural Networks

Isogeometric multilevel quadrature for forward and inverse random acoustic scattering

1. Introduction

2. Preliminary

3. SwitchNet for far field pattern

3.1. Problem setup

3.2. Low-rank property

Theorem 1.

Proof.

3.3. Matrix factorization

3.4. Neural networks

3.4.1. NN for the inverse map d→η.

3.4.2. Layers for SwitchNet

3.4.3. NN for the forward map η→d.

3.5. Numerical results

4. SwitchNet for seismic imaging

4.1. Problem setup

4.2. Low-rank property

3.4.1. NN for the inverse map $d \to η$ .

3.4.3. NN for the forward map $η \to d$ .