**3. Numerical examples**

Here we discuss the numerical computation of the Born approximation first proposed in [21] for backscattering and fixed angle data. We assume that the scatterer *h*0(*x*) is supported in the rectangle *R*<sup>0</sup> <sup>=</sup> -1, <sup>1</sup> <sup>×</sup> -1, <sup>1</sup> ⊂ℝ<sup>2</sup> . We consider the following examples:

Example 1: (linear)

$$\mu\left(\boldsymbol{\omega}, \left|\boldsymbol{\mu}\right|\right) = \chi\_{\boldsymbol{\Omega}}\left(\boldsymbol{\omega}\right) \left|\boldsymbol{\mu}\right| \frac{\left|\boldsymbol{\mu}\right|^{2}}{1 + \left|\boldsymbol{\mu}\right|^{2}} \text{ (saturation)}$$

Example 3: (linear, saturation)

$$\begin{array}{c} h\left(\mathbf{x}, \left|u\right|\right) = 0.75 \,\chi\_{\Omega\_{\parallel}}\left(\mathbf{x}\right) \frac{\left|u\right|^{2}}{1 + \left|u\right|^{2}} + \chi\_{\Omega\_{\parallel}}\left(\mathbf{x}\right) \left|u\right|\\ \text{(saturation, linear)}. \end{array}$$

Here Ω, Ω1, Ω<sup>2</sup> are the shifted ellipses rotated in angle Θ theta (counter-clockwise) detailed in **Table 1** and shown in **Figure 1**.


**Table 1.** Geometries of the ellipses Ω, Ω1 and Ω<sup>2</sup>

**Figure 1.** Geometries of the ellipses Ω, Ω1 and Ω2.

Consider the Born approximation for full data (11) in the form

$$\int\_{\mathbb{R}^2} e^{ik(\theta' - \theta, y)} f(\mathbf{y}) d\mathbf{y} = A(k, \theta', \theta), \ f = q\_B. \tag{19}$$

To discretize the unknown function *f* we divide the rectangle *R*0 into *N* = *n* × *n* disjoint subrectangles *rj* of equal size, i.e.

$$R\_0 = \bigcup\_{j=1}^{N} \tau\_j.$$

Then, we represent *f* on *R*0 in piecewise constant form:

$$f(\mathbf{y}) = \sum\_{j=1}^{N} f\_j \chi\_{r\_j}(\mathbf{y}).$$

where *fj* 's are the unknown values on *rj* 's. Substitution into (19) yields

**Semi axis Θ Centre** Ω 1/2, 1/5 17*π*/18 (−0.3, − 0.4) Ω<sup>1</sup> 1/2, 1/4 *π*/3 (0.5, 0) Ω<sup>2</sup> 1/3, 1/5 2*π*/3 (−0.5, 0)

**Table 1.** Geometries of the ellipses Ω, Ω1 and Ω<sup>2</sup>

130 Applied Linear Algebra in Action

**Figure 1.** Geometries of the ellipses Ω, Ω1 and Ω2.

of equal size, i.e.

Then, we represent *f* on *R*0 in piecewise constant form:

subrectangles *rj*

Consider the Born approximation for full data (11) in the form

To discretize the unknown function *f* we divide the rectangle *R*0 into *N* = *n* × *n* disjoint

$$\sum\_{j=1}^{N} f\_j \int\_{r\_j} e^{i k \left(\theta' - \theta, \mathbf{y}\right)} d\mathbf{y} = A(k, \theta', \theta).$$

Evaluating this at some points *<sup>k</sup>* <sup>=</sup>*ks*, *<sup>θ</sup>* ′ =*θ<sup>t</sup>* ′ and *θ* = *θp* leads us to

$$\sum\_{j=1}^{N} f\_j \int\_{r\_j} e^{ik\_s(\theta\_t' - \theta\_p, y)} dy = A(k\_s, \theta\_t', \theta\_p), \qquad s = 1, \dots, N\_1, t = 1, \dots, N\_2, p = 1, \dots, N\_3.$$

If we denote *M* = *N*1*N*2*N*3, we may form one linearized index *l* = *l*(*s*, *t*, *p*), *l* = 1, …, *M* and write the latter equation as

$$\sum\_{j=1}^{N} f\_j E(j, l) = g\_l, \ l = 1, \dots, M,\tag{20}$$

where *E*( *j*, *l*)=*∫ rj e iks* (*θt* ′ -*θp*,*y*) d*y* are easily evaluated and *gl* = *A*(*ks*, *θ<sup>t</sup>* ′ , *θp*) needs to be computed. The computation of *gl* is carried out using numerical integration, see [21] for details. In matrix form (20) is clearly *Ef* = *g*.

For backscattering data and fixed angle data the system (20) is modified accordingly. We note that the system *Ef* = *g* does not depend on the scatterer but only on data type and measurement setup.

The fixed energy case differs from the first three data sets. In fixed energy inversion we approximate the scattering transform as

$$
\mathcal{T}\_h(\xi) \approx \mathcal{T}\_{h,1}(\xi).
$$

We choose *M* = *m* × *m* points *ξ* uniformly from the rectangle [−*s*, *s*] × [−*s*, *s*]. The function *Th*,1(*ξ*) is evaluated by numerical integration, see [19,21] for details. Then the inverse Born approxi‐ mation (17) is computed similarly to (19).

We use the following parameter values:

(19)

$$k\_s = s \text{ } , \theta\_p = \,^\rho\theta\_p' = \left(\cos\frac{p2\pi(N\_2 - 1)}{N\_2}, \sin\frac{p2\pi(N\_2 - 1)}{N\_2}\right)$$

We use *N* = 2500 and *N*1 = 12 for each data set. For full data we use *N*2 = *N*3 = 6. For backscattering and fixed angle data, we use *N*2 = 24. For fixed energy scattering we use *m* = 40 and *s* = 6.

In all four cases we obtain the linear system *Ef* = *g* whose coefficient matrix *E* is of size *M* × *N*. The data *g* is corrupted with zero mean Gaussian noise with standard deviation *σ* = 0.01 max |*g*|.

The size *M* as well as the ranks *r*(*E*) and (approximate) condition numbers log10*κ*(*E*) measuring the ill-posedness of the linear system *EF* = *g* are shown in **Table 2**.


**Table 2.** Matrix sizes, ranks and condition numbers

As the linear system is rank-deficient and ill-posed, we use regularization method to solve it. More precisely, we apply the total variation method (TV) which is defined as

$$f = \arg\min\_{\mathbf{v}} \{ \left\| \left| Ez - g \right| \right\|\_{2}^{2} + \delta \left\| \left| Lz \right| \right\|\_{1} \} \,,$$

where the matrix *L* implements the differences between neighbouring elements in horizontal and vertical directions (for details, see [21]). As in [21] we formulate this minimization problem as a quadratic problem in standard form for more efficient solution. As the regularization parameter we use *δ* = 2 ⋅ 10− 3 for *DE* and *δ* = 10− 3 otherwise.

All computations are carried out in a UNIX system with 512 GB of memory and 40 logical CPU cores, each running at 2.8 GHz. The software platform is MATLAB R2015b. We have used 12 workers in parallel in computing the right-hand side *g*. We note that a desktop PC with 8 logical cores running at 3.4 GHz and 16 GB of memory is also capable of our computations, but with 7 workers it is considerably slower in computing *g*.

The computational times to both form and solve the linear system are shown in **Table 3**. We point out that the right-hand side *g* contains synthetic data and actual physical measurements may take longer or shorter time.


**Table 3.** Computational times

We use *N* = 2500 and *N*1 = 12 for each data set. For full data we use *N*2 = *N*3 = 6. For backscattering and fixed angle data, we use *N*2 = 24. For fixed energy scattering we use *m* = 40 and *s* = 6.

In all four cases we obtain the linear system *Ef* = *g* whose coefficient matrix *E* is of size *M* × *N*. The data *g* is corrupted with zero mean Gaussian noise with standard deviation

The size *M* as well as the ranks *r*(*E*) and (approximate) condition numbers log10*κ*(*E*) measuring

As the linear system is rank-deficient and ill-posed, we use regularization method to solve it.

where the matrix *L* implements the differences between neighbouring elements in horizontal and vertical directions (for details, see [21]). As in [21] we formulate this minimization problem as a quadratic problem in standard form for more efficient solution. As the regularization

All computations are carried out in a UNIX system with 512 GB of memory and 40 logical CPU cores, each running at 2.8 GHz. The software platform is MATLAB R2015b. We have used 12 workers in parallel in computing the right-hand side *g*. We note that a desktop PC with 8 logical cores running at 3.4 GHz and 16 GB of memory is also capable of our computations, but with

The computational times to both form and solve the linear system are shown in **Table 3**. We point out that the right-hand side *g* contains synthetic data and actual physical measurements

the ill-posedness of the linear system *EF* = *g* are shown in **Table 2**.

**Table 2.** Matrix sizes, ranks and condition numbers

**Data** *M r***(***E***) log10***κ***(***E***)** Full 1728 360 63 Backscattering 288 281 15 Fixed angle 288 241 17 Fixed energy 1600 268 15

More precisely, we apply the total variation method (TV) which is defined as

parameter we use *δ* = 2 ⋅ 10− 3 for *DE* and *δ* = 10− 3 otherwise.

7 workers it is considerably slower in computing *g*.

may take longer or shorter time.

*σ* = 0.01 max |*g*|.

132 Applied Linear Algebra in Action

The contour plots of scatterers *h*0(*x*) and their TV reconstructions via Born approximation from full data, backscattering data and fixed angle data (with *θ*0 = (1, 0)) are shown in **Figures 2**–**5** for all examples. For fixed energy scattering we only consider the linear Example 1, since otherwise we do not have direct comparison to a scatterer. The TV reconstruction is shown in **Figure 6**. In each figure solid white line indicates the true geometry of the scatterer.

We see that the location of the scatterer is located quite nicely in all cases. The shape of the scatterer is best seen from full data and backscattering data. By computing the Born approxi‐ mation from full angle data, we close an open problem from [5].

**Figure 2.** Scatterer *h*0(*x*) and its TV reconstruction via Born approximation, Example 1. (a) scatterer; (b) full data; (c) backscattering data; (d) fixed angle data.

**Figure 3.** Scatterer *h*0(*x*) and its TV reconstruction via Born approximation, Example 2. (a) scatterer; (b) full data; (c) backscattering data; (d) fixed angle data

**Figure 4.** Scatterer *h*0(*x*) and its TV reconstruction via Born approximation, Example 3. (a) scatterer; (b) full data; (c) backscattering data; (d) fixed angle data.

**Figure 5.** Scatterer *h*0(*x*) and its TV reconstruction via Born approximation, Example 4. (a) scatterer; (b) full data; (c) backscattering data; (d) fixed angle data.

**Figure 6.** Scatterer *h*0(*x*) and its TV reconstruction via Born approximation, Example 1. (a) scatterer; (b) fixed energy.
