Numerical Solution of Bilinear Inverse Problems

Frank Natterer

November 1995

1 Introduction

This note is concerned with a quite down-to-earth approach to the numerical solution of inverse problems. These problems often have a bilinear structure. We assume that our inverse problem is governed by a family of bilinear equations

A_{j}u_{j} + B_{j}(f, u_{j}) = a_{j} + b_{j}f , j = 1, …, p  (1.1)

where A_{j} : U_{j} → Y_{j} is linear, B_{j} : F × U_{j} → Y_{j} bilinear, a_{j} ∈ Y_{j}, b_{j} : F → Y_{j} is linear. Here, U_{j} is the Hilbert space of states, F the Hilbert space of parameters, and Y_{j} another Hilbert space. The subscript j refers to a projection direction or source position in transmission tomography or to a electrode pair in impedance tomography. Considering only one of the equations (1.1) at a time is an essential feature of our method. We assume that (1.1) is uniquely solvable for u_{j} if f is given. We consider the problem of determining f from the additional equations

C_{j}u_{j} + D_{j}(f, u_{j}) = c_{j} + d_{j}f , j = 1, …, p .  (1.2)

C_{j}, D_{j}, c_{j} and d_{j} are exactly as A_{j}, B_{j}, a_{j}, b_{j} in (1.1) with Y_{j} replaced by a third Hilbert space Z_{j}, the space of observations. Let u_{j}(f) be the solution of (1.1) for a certain f ∈ F. Then our problem amounts to solving the equations

R_{j}(f) = 0 ,  j = 1, …, p ,  (1.3)
R_{j}(f) = C_{j}u_{j}(f) + D_{j}(f, u_{j}(f)) − c_{j} − d_{j}f

for f.

A standard method for solving the nonlinear problem (1.3) is the Newton or Gauss-Newton method. These methods necessitate the computation of the operators (R′_{j}(f)^{∗}R′_{k}(f))_{j,k=1,...,p}^{−1}, where R′_{j}(f) is the Fréchet derivative of R_{j}(f). For problems of practical relevance this is totally out of question. An obvious way to circumvent this difficulty would be to compute a minimal norm solution to just one of the linearized equations

R_{j}(f) + R′_{j}(f)h_{j} = 0 ,

i.e.

h_{j} = −R′_{j}(f)^{∗}C_{j}^{−1}R_{j}(f) , C_{j} = R′_{j}(f)R′_{j}(f)^{∗}

and to update f according to f ← f + ωh_{j} with a relaxation factor ω. However, the computational burden for C_{j} is still by far too high. Thus we go one step further, replacing C_{j} by an easy to evaluate operator, e.g. by the identity operator. Then, one step of our algorithm, transforming the approximation f^{0} into a new approximation f^{1}, reads:

(i) f_{0} = f^{0}
(ii) For j = 1, 2, …, p

(1.4)

f_{j} = f_{j−1} − ωh_{j} , h_{j} = R′_{j}(f_{j−1})^{∗}C_{j}^{−1}R_{j}(f_{j−1})

(iii) f^{1} = f_{p} .  (1.5)

If R_{j} is a linear operator and C_{j} is simply the diagonal of R_{j}R_{j}^{∗}, this iteration is the well-known algebraic reconstruction technique of computerized tomography, or the (block) Kaczmarz method, for solving underdetermined linear systems. We shall prove that (1.4) converges subject to natural assumptions for 0 < ω < 2 provided that C_{j} is positive definite and C_{j} ≥ R′_{j}(f)R′_{j}(f)^{∗}. In the general nonlinear case, we do not have any convergence results. But it is easy to see what the method actually does. Suppose the problem (1.3) has a solution f - not necessarily unique. Then, in F,

∥f_{j} − f∥^{2} = ∥f_{j−1} − f − ωh_{j}∥^{2} = ∥f_{j−1} − f∥^{2} − 2ω(f_{j−1} − f, h_{j}) + ω^{2}∥h_{j}∥^{2} .

The factor of ω is

(f_{j} − f, h_{j}) = (C_{j}^{−1}R′_{j}(f_{j−1})(f_{j−1} − f) , R_{j}(f_{j−1}))
= (C_{j}^{−1}R_{j}(f_{j−1}) , R_{j}(f_{j−1})) + ⋯

where the dots stand for higher order terms in f_{j−1} − f. Thus, if C_{j} is positiv definit and ω sufficiently small, we can expect f_{j} to be closer to f than f_{j−1}, provided R_{j}(f_{j−1}) ≠ 0. The decrease of the error is best possible if

ω = \frac{(C_{j}^{−1}R_{j}(f_{j−1}), R_{j}(f_{j−1}))}{∥R′_{j}(f_{j−1})^{∗}C_{j}^{−1}R_{j}(f_{j−1})∥^{2}}

However, the emphasis of this paper is not so much on theoretical questions

such as convergence. Rather we describe the method to some detail for a variety of inverse problems for partial differential equations. We also present numerical examples which demonstrate that the method is quite easy to implement and, compared to other methods, very efficient. Its success depends of course on the degree of ill-posedness of the inverse problem which varies very much within the class of problems treated here. Also, the method is quite plausible, making use of vivid phenomena such as backpropagation. The method is closely related to the method of adjoint fields which is used in the engineering literature [1].

We do not claim any originality of our approach. Similar methods have been used in thousends of papers. A typical example is [2]. An implementation for ultrasound tomography is given in [4].

2 Linerarization

We linearize (1.1)-(1.2) about a solution u_{j} = u_{j}(f) of (1.1). Replacing u_{j} by u_{j} + w_{j}, f by f + h, (1.1) reads

A_{j}(u_{j} + w) + B_{j}(f + h, u_{j} + w) = a_{j} + b_{j}(f + h) .

Neglecting higher order terms and using that u_{j}, f solve (1.1) we obtain

A_{j}w_{j} + B_{j}(f, w_{j}) = b_{j}h − B_{j}(h, u_{j}) .  (2.1)

Similarly, neglecting higher order terms in (1.2) yields

R_{j}(f + h) = R_{j}(f) + C_{j}w_{j} + D_{j}(f, w_{j}) + D_{j}(h, u_{j}) − d_{j}h .

Thus the operator R′_{j}(f) : F → Z_{j} is given by

R′_{j}(f)h = C_{j}w_{j} + D_{j}(f, w_{j}) + D_{j}(h, u_{j}) − d_{j}h  (2.2)

where w_{j} is the solution of (2.1).

In the application to concrete problems, the bulk of the analytical work consists in determining the adjoint operator (R′_{j}(f))^{∗} : Z_{j} → F. While an explicit formula for this operator in the functional analytic framework is easy to derive, this formula is not very useful in practice. Thus we derive (R′_{j}(f))^{∗} for specific cases in the next section.

3 Special inverse problems

Now we derive the operator (R′_{j}(f))^{∗} for special inverse problems of partial differential equations. The derivation is in no way rigouros. We neglect subtleties such as the exact domains of definition. The essential tools are Gauss’ integral theorem and Green’s formula.

1. Impedance tomography.

Let Ω ⊂ IR^{n}, n = 2, 3 be sufficiently regular, and let σ be a function in Ω, the conductivity of the material in Ω. We want to determine σ from measurements of currents f_{j} and voltages g_{j}, j = 1, …, p, on ∂Ω. The currents are usually applied at pairs of electrodes on ∂Ω.

The mathematical model is as follows. The potential u_{j} satisfies

div(σ∇u_{j}) = 0 in Ω , σ \frac{∂u_{j}}{∂ν} = f_{j} on ∂Ω .  (3.1)

We normalize u_{j} by requiring that the mean value on ∂Ω is zero. The measurements provide us with the values of

u_{j} = g_{j} on ∂Ω .  (3.2)

We have to determine σ from (3.1)-(3.2) for j = 1, …, p. We assume σ to be known on ∂Ω.

It is clear that (3.1) can be written in the form of (1.1) with U_{j} = L_{2}(Ω), Y_{j} = L_{2}(Ω) × L_{2}(Ω), F = L_{2}(Ω), and σ playing the role of f. Also, (3.2) is a special case of (1.2) with Z_{j} = L_{2,σ}(∂Ω) which is simply L_{2}(∂Ω) with the (known) scalar product ∫_{∂Ω} σuvds, and we have

R_{j}(σ) = u_{j}|_{∂Ω} − g_{j}

where u_{j} is the solution of (3.1). Hence, (2.2) reads

R′_{j}(σ)η = w_{j}|_{∂Ω} ,
div(σΔw_{j}) = −div(η∇u_{j}) in Ω , σ \frac{∂w_{j}}{∂ν} = 0 on ∂Ω .  (3.3)

Theorem 3.1 The operator (R′_{j}(σ))^{∗} : L_{2,σ}(∂Ω) → L_{2}(Ω) is given by

(R′_{j}(σ))^{∗}g = ∇u_{j} · ∇z  (3.4)

where z is the solution of

div(σ∇z) = 0 in Ω , σ \frac{∂z}{∂ν} = −g on ∂Ω .  (3.5)

Proof: For any z, w we have

∫_{Ω} (div(σ∇w)z − w div(σ∇z))dx = ∫_{∂Ω} σ (\frac{∂w}{∂ν}z − w\frac{∂z}{∂ν} ) ds .

Choosing w = w_{j} from (3.3) and z from (3.5) we obtain

− ∫_{Ω} div(η∇u_{j})zdx = − ∫_{∂Ω} σ(R′_{j}(σ)η)gds .

Since η = 0 on ∂Ω an integration by parts yields

∫_{Ω} η∇u_{j} · ∇zdx = ∫_{∂Ω} σ(R′_{j}(σ)η)gds

or

(η, ∇u_{j} · ∇z)_{L_{2}(Ω)} = (R′_{j}(σ)η, g)_{L_{2,σ}(∂Ω)} .

This holds for each η vanishing on ∂Ω, verifying (3.4).

□

In the light of the theorem, the iteration (1.4) proceeds as follows. For j = 1, 2, …, p

Solve div (σ_{j−1}∇u_{j}) = 0 in Ω , σ_{j−1} \frac{∂u_{j}}{∂ν} = f_{j} on ∂Ω.

Solve div (σ_{j−1}∇z_{j}) = 0 in Ω , σ_{j−1} \frac{∂z_{j}}{∂ν} = u_{j} − g_{j} on ∂Ω.

Put σ_{j} = σ_{j−1} + ω∇u_{j} · ∇z_{j} .

We see that each sweep of the iteration requires the solution of 2p boundary value problems of nearly identical shape. This is not only easy to program but also very cheap computationally, compared to other methods.

2. A parabolic inverse problem.

Consider the initial-boundary value problem

\frac{∂u_{j}}{∂t} − div(σ∇u_{j}) = q_{j} in Ω × [0, T]

(3.6)

u_{j}(x, 0) = 0 for x ∈ Ω , \frac{∂u_{j}}{∂ν} (x, t) = 0 for x ∈ ∂Ω .

The source terms q_{j} are assumed to be known, while the diffusion coefficient σ(x) has to be determined from the knowledge of u_{j} on ∂Ω:

u_{j}(x, t) = g_{j}(x, t) , x ∈ ∂Ω , t ∈ [0, T] .  (3.7)

Again we assume σ to be known on ∂Ω.

This problem is of the form (1.1)-(1.2). We have

U_{j} = \{u ∈ L_{2}(Ω × [0, T]) : u = 0 for t = 0 , \frac{∂u}{∂ν} = 0 on ∂Ω\} ,
F = L_{2}(Ω) ,  Z_{j} = L_{2,σ}(∂Ω × [0, T]) .

The operator R_{j} : L_{2}(Ω) → L_{2,σ}(∂Ω × [0, T]) is given by

R_{j}(σ) = u_{j} − g_{j} on ∂Ω × [0, T]

where u_{j} is the solution of (3.6). According to (2.2), the Fréchet derivative is

R′_{j}(σ)η = w_{j}|_{∂Ω×[0,T]} ,
\frac{∂w_{j}}{∂t} − div(σ∇w_{j}) = div(η∇u_{j}) in Ω × [0, T] ,  (3.8)

w_{j} = 0 for t = 0 , \frac{∂w_{j}}{∂ν} = 0 on ∂Ω × [0, T] .

Theorem 3.2 The operator (R′_{j}(σ))^{∗} : L_{2,σ}(∂Ω × [0, T]) → L_{2}(Ω) is given by

(R′_{j}(σ))^{∗}g(x) = ∫_{0}^{T} ∇u_{j} · ∇z(x, t)dt

where z is the solution of

\frac{∂z}{∂t} + div(σ∇z) = 0 in Ω × [0, T] ,

(3.9)

z = 0 for t = T , \frac{∂z}{∂ν} = −g on ∂Ω × [0, T] .

Proof: For any z, w we have

∫_{0}^{T} ∫_{Ω} \{(\frac{∂w}{∂t} − div(σ∇w))z + (\frac{∂z}{∂t} + div(σ∇z))w\} dxdt

= ∫_{Ω} ((wz)(x, T) − (wz)(x, 0))dx + ∫_{0}^{T} ∫_{∂Ω} σ(\frac{∂z}{∂ν} w − z \frac{∂w}{∂ν})dsdt .

For w_{j}, z from (3.8)-(3.9) this reads

∫_{0}^{T} ∫_{Ω} div(η∇u_{j})zdxdt = − ∫_{0}^{T} ∫_{∂Ω} σgR′_{j}(σ)ηdsdt .

An integration by parts on the left hand side turns this into

∫_{Ω} η(x) ∫_{0}^{T} ∇u_{j} · ∇zdtdx = ∫_{0}^{T} ∫_{∂Ω} σgR′_{j}(σ)ηdsdt ,

hence the theorem.

□

The algorithm (1.4) is as follows: For j = 1, 2, …, p Solve

\frac{∂u_{j}}{∂t} − div(σ_{j−1}∇u_{j}) = q_{j} in Ω × T , u_{j} = 0 for t = 0 , \frac{∂u_{j}}{∂ν} = 0 on ∂Ω

and

\frac{∂z_{j}}{∂t} +div(σ_{j−1}∇z_{j}) = 0 in Ω×T , z_{j} = 0 for t = T , \frac{∂z_{j}}{∂ν} = g_{j}−u_{j} on ∂Ω .

Put

σ_{j}(x) = σ_{j−1}(x) + ω ∫_{0}^{T} ∇u_{j}(x, t) · ∇z_{j}(x, t)dt .

Alternatively we may consider the boundary condition u_{j} = 0 on ∂Ω × [0, T] instead of ∂u_{j}/∂ν = 0 in (3.6), and ∂u_{j}/∂ν = g_{j} on ∂Ω × [0, T] as data in (3.7). Then,

R_{j}(σ) = \frac{∂u_{j}}{∂ν} − g_{j} ,

and

R′(σ)η = \frac{∂w_{j}}{∂ν} |_{∂Ω×[0,T]}

where now w_{j} satisfies (3.8) with the boundary condition w_{j} = 0 on ∂Ω × [0, T]. The adjoint is given by the same expression as above, but with the boundary condition z = g on ∂Ω × [0, T].

3. Laser tomography: Stationary case Here the governing equation is

θ · ∇u(x, θ) + a(x)u(x, θ) − b(x) ∫_{S_{n−1}} η(θ · θ′)u(x, θ′)dθ′ = q(x − s) in Ω × S^{n−1},

(3.10)

u = 0  on (∂Ω × S^{n−1})_{−} ,

where

(∂Ω × S^{n−1})_{±} = \{(x, θ) ∈ ∂Ω × S^{n−1} : ν(x) · θ \begin{matrix}≥ \\ ≤ \end{matrix} 0\}

with ν(x) the exterior normal to ∂Ω at x ∈ ∂Ω. The source term q is deltalike, and the sources s are on ∂Ω. η is assumed to be known. We want to recover a, b from knowing u on ∂Ω for sources s_{1}, …, s_{p}. Denoting by u_{j} the solution of (3.10) for s = s_{j} we have

θ · ∇u_{j}(x, θ) + a(x)u_{j}(x, θ) − b(x) ∫_{S^{n−1}} η(θ · θ′)u_{j}(x, θ′)dθ′ = q(x − s_{j}) in Ω × S^{n−1},

(3.11)

u_{j} = 0  on (∂Ω × S^{n−1})_{−} .

The information on u on ∂Ω which we want to make use of is

∫_{S^{n−1}} ν(x) · θu_{j}(x, θ)dθ = g_{j}(x) , x ∈ ∂Ω .  (3.12)

Obviously, (3.11)-(3.12) is of the type of (1.1)-(1.2) with f = (a, b) and

U_{j} = {u ∈ L_{2}(Ω × S^{n−1}) : u = 0 on (∂Ω × S^{n−1})_{−}} ,
F = L_{2}(Ω) × L_{2}(Ω) , Z_{j} = L_{2}(∂Ω) .

The operator R_{j} : L_{2}(Ω) × L_{2}(Ω) → L_{2}(∂Ω) is given by

R_{j}( \begin{matrix}a \\ b \end{matrix} )(x) = ∫_{S^{n−1}} ν(x) · θu_{j}(x, θ)dθ − g_{j}(x)

where u_{j} is the solution to (3.11). The linearized operator is

R′_{j}( \begin{matrix}a \\ b \end{matrix} )( \begin{matrix}h \\ k \end{matrix} ) = ∫_{S^{n−1}} ν(x) · θw(x, θ)dθ

where w solves

θ · ∇w(x, θ) + a(x)w(x, θ) − b(x) ∫_{S^{n−1}} η(θ · θ′)w(x, θ′)dθ

= k(x) ∫_{S^{n−1}} η(θ · θ′)u_{j}(x, θ′)dθ′ − h(x)u_{j}(x, θ) in Ω × S^{n−1} , (3.13)

w = 0 on (∂Ω × S^{n−1})_{−} .

Theorem 3.3 The operator (R′_{j}( \begin{matrix}a \\ b \end{matrix} ))^{∗} : L_{2}(∂Ω) → L_{2}(Ω) × L_{2}(Ω) is given by

(R′_{j}( \begin{matrix}a \\ b \end{matrix} ))^{∗}g(x) = ∫_{S^{n−1}} ( \begin{matrix}−u_{j}(x, θ) \\ ∫_{S^{n−1}} η(θ · θ′)u_{j}(x, θ′)dθ′ \end{matrix} ) z(x, θ)dθ

where z is the solution of

− θ · ∇z(x, θ) + a(x)z(x, θ) − b(x) ∫_{S^{n−1}} η(θ · θ′)z(x, θ′)dθ′ = 0 in Ω × S^{n−1} ,

(3.14)

z = g  on (∂Ω × S^{n−1})_{+} .

Proof: For any z, w we have

∫_{Ω} ∫_{S^{n−1}} (θ · ∇w(x, θ) + a(x)w(x, θ) − b(x) ∫_{S^{n−1}} η(θ · θ′)w(x, θ′)dθ′)z(x, θ)dxdθ
= ∫_{Ω} ∫_{S^{n−1}} (−θ · ∇z(x, θ) + a(x)z(x, θ) − b(x) ∫_{S^{n−1}} η(θ′ · θ)z(x, θ′)dθ′)w(x, θ)dxdθ
+ ∫_{∂Ω} ∫_{S^{n−1}} z(x, θ)w(x, θ)ν(x) · θdsdθ .

For z, w from (3.13)-(3.14) this implies

∫_{Ω} ∫_{S^{n−1}} (k(x) ∫_{S^{n−1}} η(θ · θ′)u_{j}(x, θ′)dθ′ − h(x)u_{j}(x, θ)) z(x, θ)dxdθ

= ∫_{∂Ω} g(x)R′_{j}( \begin{matrix}a \\ b \end{matrix} )( \begin{matrix}h \\ k \end{matrix} )ds ,

hence the theorem.

□

The algorithm is as follows: For j = 1, …, p Solve

θ · ∇u_{j}(x, θ) + a_{j−1}(x)u_{j}(x, θ) − b_{j−1}(x) ∫_{S^{n−1}} η(θ · θ′)u_{j}(x, θ)dθ′ = q(x − s_{j}) in Ω × S^{n−1} ,
u_{j} = 0  on (∂Ω × S^{n−1})_{−}

and

−θ · ∇z_{j}(x, θ) + a_{j−1}(x)z_{j}(x, θ) − b_{j−1}(x) ∫_{S^{n−1}} η(θ · θ′)z_{j}(x, θ′)dθ′ = 0 in Ω × S^{n−1} ,
z_{j} = u_{j} − g_{j}  on (∂Ω × S^{n−1})_{+} .

Put

a_{j} = a_{j−1} − w ∫_{S^{n−1}} u_{j}(x, θ)z_{j}(x, θ)dθ ,
b_{j} = b_{j−1} + w ∫_{S^{n−1}} ∫_{S^{n−1}} η(θ · θ′)u_{j}(x, θ′)z_{j}(x, θ)dθ′dθ .

4. SPECT In SPECT one has to compute the source distribution f from

θ_{j} · ∇u_{j}(x) + a(x)u_{j}(x) = f(x) in Ω  (3.15)
u_{j} = 0 on ∂Ω_{j}^{−} ,
u_{j} = g_{j} on ∂Ω_{j}^{+}  (3.16)

where ∂Ω_{j}^{±} = \{x ∈ ∂Ω : ν(x) · θ_{j} \begin{matrix}≥ \\ ≤ \end{matrix} 0\}.

Here, Ω is a bounded domain in IR^{2}, and θ_{1}, …, θ_{p} are direction vectors ∈ S^{1}. If a is known, then (3.15)-(3.16) reduces to a linear problem which can be solved by inverting the attenuated Radon transform. We consider the case in which a is unknown, too.

(3.15)-(3.16) is of the form (1.1)-(1.2) again, with (a, f) playing the role of f and

U_{j} = {u ∈ L_{2}(Ω) : u = 0 on ∂Ω_{j}^{−}}
F = L_{2}(Ω) × L_{2}(Ω) , Z_{j} = L_{2}(∂Ω_{j}^{+}) .

The operator R_{j}( \begin{matrix}a \\ f \end{matrix} ) : L_{2}(Ω) × L_{2}(Ω) → L_{2}(∂Ω_{j}^{+}) is given by

R_{j}( \begin{matrix}a \\ f \end{matrix} )(x) = u_{j}(x) − g_{j}(x) , x ∈ ∂Ω_{j}^{+}

where u_{j} is the solution of (3.15). The linearized operator is

R′_{j}( \begin{matrix}a \\ f \end{matrix} )( \begin{matrix}h \\ k \end{matrix} ) = w_{j}|_{∂Ω{+ \atop j}}

where w_{j} solves

θ_{j} · ∇w_{j} + aw_{j} = k − hu_{j} in Ω , w = 0 on ∂Ω_{j}^{−} .  (3.17)

Theorem 3.4 The operator (R′_{j}( \begin{matrix}a \\ f \end{matrix} ))^{∗} : L_{2}(∂Ω_{j}^{+}) → L_{2}(Ω)×L_{2}(Ω) is given

by

(R′_{j}( \begin{matrix}a \\ f \end{matrix} ))^{∗}g = ( \begin{matrix}−u_{j}z_{j} \\ z_{j} \end{matrix} )

where z_{j} is the solution of

− θ_{j} · ∇z_{j} + az_{j} = 0 in , z_{j} = g ν(x) · θ_{j} on ∂Ω_{j}^{+} .  (3.18)

Proof: For any z, w we have

∫_{Ω} (θ_{j} · ∇w + aw)zdx = ∫_{Ω} (−θ_{j} · ∇z + az)wdx + ∫_{∂Ω} θ_{j} · ν(x)wzds .

For z, w as in (3.13)-(3.17) this means

∫_{Ω} (k − hu_{j})z_{j}dx = ∫_{∂Ω_{j}^{+}} R′_{j}( \begin{matrix}a \\ f \end{matrix} )( \begin{matrix}h \\ k \end{matrix} )ds ,

hence the theorem.

□

The algorithm is:

For j = 1, …, p

Solve

θ_{j} · ∇u_{j} + a_{j−1}u_{j} = f_{j−1} in Ω , u_{j} = 0 on ∂Ω_{j}^{−}  (3.19)

and

− θ_{j} · ∇z_{j} + a_{j−1}z_{j} = 0 in Ω , z_{j} = u_{j} − g_{j} on ∂Ω_{j}^{+} .  (3.20)

Put

\begin{aligned}a_{j} & = & a_{j−1} − ωu_{j}z_{j} , \\ f_{j} & = & f_{j−1} + ωz_{j} . \end{aligned}

The forward problem (3.19) can be solved by inverting the attenuated Radon transform, and (3.20) is the attenuated backprojection. Thus the algorithm is very similar to the ART algorithm in computerized tomography and in fact reduces to this algorithm if a is known, e.g. a = 0.

5. Ultrasound tomography: Time harmonic case

Here we have the Helmholtz equation

Δu + k^{2}(1 − f)u = 0 in IR^{n}

where u = e^{ikx·θ} + v with the scattered wave v satisfying the radiation condition

lim_{r→∞} (\frac{∂v}{∂r} − ikv) = 0 .

Here, f is a complex valued function which has to be recovered from knowing u for each direction θ ∈ S^{n−1} outside the support of f. k is a real parameter controlling the spatial resolution.

Let Ω be a ball containing supp(f) in its interior. We rewrite the radiation condition as a boundary condition on ∂Ω, to wit

\frac{∂v}{∂ν} = Bv on ∂Ω

with a linear operator B on L_{2}(∂Ω). Assume that the measurements are made for finitely many directions θ_{1}, …, θ_{p}, and let u_{j} = e^{ikx·θ_{j}} + v_{j} be the solution for direction θ_{j}. Then,

Δv_{j} + k^{2}(1 − f)v_{j} = k^{2}fe^{ikx·θ_{j}} in Ω

(3.21)

\frac{∂v_{j}}{∂ν} = Bv_{j} on ∂Ω

and f has to be determined from

v_{j} = g_{j} on ∂Ω .  (3.22)

Obviously, (3.21)-(3.22) is of the form (1.1)-(1.2), with

U_{j} = \{v ∈ L_{2}(Ω) : \frac{∂v}{∂ν} = Bv on ∂Ω\} ,
F = L_{2}(Ω) ,  Z_{j} = L_{2}(∂Ω) .

The operator R_{j} : L_{2}(Ω) → L_{2}(∂Ω) is

R_{j}f = v_{j} − g_{j} on ∂Ω

where v_{j} solves (3.21). The linearized operator is given by R′_{j}(f)h = w_{j}, where w_{j} is a solution of

Δw_{j} + k^{2}(1 − f)w_{j} = k^{2}hu_{j} in Ω  (3.23)

\frac{∂w_{j}}{∂ν} = Bw_{j} on ∂Ω .

Theorem 3.5 The operator (IR′_{j}(f))^{∗} : L_{2}(∂Ω) → L_{2}(Ω) is given by

(R′_{j}(f))^{∗}g = k^{2}\overline{u}_{j}z

where z is the solution of

Δz + k^{2}(1 − \overline{f})z = 0 in Ω , \frac{∂z}{∂ν} − B^{∗}z = −g on ∂Ω .  (3.24)

Proof: For any z, w we have

∫_{Ω} \{(Δw + k^{2}(1 − f)w)\overline{z} − w(\overline{Δz + k^{2}(1 − \overline{f})z)}\}dx

= ∫_{∂Ω} (\frac{∂w}{∂ν}\overline{z} − w\frac{∂\overline{z}}{∂ν} ) ds .

For the functions z, w_{j} from (3.24)-(3.23) this yields

∫_{Ω} k^{2}(e^{ikx·θ_{j}} + v_{j})\overline{z}dx = ∫_{∂Ω} ((Bw)\overline{z} − w_{j}(\overline{B^{∗}z − g}))ds
= ∫_{∂Ω} w_{j}\overline{g}ds
= ∫_{∂Ω} (R′_{j}(f)h)\overline{g}ds

or

(h, k^{2}(e^{−ikx·θ_{j}} + \overline{v}_{j})z) = (R′_{j}(f)h, g)_{L_{2}(∂Ω)} .
□

6. Ultrasound tomography: Time resolved case Here we have the wave equation

\frac{∂^{2}u}{∂t^{2}} = c^{2}Δu + q(x − s)δ(t) in Ω × [0, T]  (3.25)

with initial and boundary conditions

u = \frac{∂u}{∂t} = 0 on Ω × \{0\} , \frac{∂u}{∂ν} = 0 on ∂Ω × [0, T]  (3.26)

with ν the exterior normal on ∂Ω. The source term q is δ like and the sources s sit on ∂Ω. The problem is to recover the sound speed c in Ω from the measurement of u on ∂Ω × [0, T], or a part thereof, for sources s_{1}, …, s_{p}. We assume c^{2} = 1 + f with f = 0 on ∂Ω. Again, this problem is obviously a special case of (1.1), (1.2) with

Z_{j} = Z = L_{2}(∂Ω × [0, T]) , F = L_{2}(Ω) .

The linearized operator R′_{j}(f) : F → Z is given by

R′_{j}(f)h = w_{j}|_{∂Ω×[0,T]} ,

where w_{j} is the solution of

\frac{∂^{2}w_{j}}{∂t^{2}} = (1 + f)Δw_{j} + hΔu_{j}  in Ω × [0, T] ,

(3.27)

w_{j} = \frac{∂w_{j}}{∂t} = 0 on Ω × \{0\} , \frac{∂w_{j}}{∂ν} = 0 on ∂Ω × [0, T]

with u_{j} the solution to (3.25), (3.26) for c^{2} = 1 + f and s = s_{j}.

Theorem 3.6 The operator (R′_{j}(f))^{∗} : L_{2}(∂Ω × [0, T]) → L_{2}(Ω) is given by

(R′_{j}(f))^{∗}g = ∫_{0}^{T} Δu_{j} · zdt

where z is the solution to

\frac{∂^{2}z}{∂t^{2}} = Δ((1 + f)z) in Ω × [0, T]  (3.28)

z = \frac{∂z}{∂t} = 0 on Ω × \{T \} , \frac{∂z}{∂ν} = g on ∂Ω × [0, T] .

Proof: For arbitrary w, z we have

∫_{Ω} ∫_{0}^{T} (\frac{∂^{2}w}{∂t^{2}} − (1 + f)Δw) zdtdx = ∫_{Ω} ∫_{0}^{T} (\frac{∂^{2}z}{∂t^{2}} − Δ((1 + f)z)) wdtdx

− ∫_{∂Ω} ∫_{0}^{T} (\frac{∂w}{∂ν}z − w\frac{∂z}{∂ν} ) dtdx + [∫_{Ω} (\frac{∂w}{∂t}z − w\frac{∂z}{∂t} ) dx]_{0}^{T} .

For w = w_{j} from (3.27) and z as in (3.28) this reads

∫_{Ω} h ∫_{0}^{T} Δu_{j} · zdtdx = ∫_{∂Ω} ∫_{0}^{T} R′_{j}(f)h · gdtdx ,

hence

(R′_{j}(f))^{∗}g = ∫_{0}^{T} Δu_{j} · zdt

□

If u is measured only on a part ∂\overline{Ω} × [0, T] of ∂Ω × [0, T], then Z has to be replaced by \overline{Z} = L_{2}(∂\overline{Ω} × [0, T]). With E the extension from ∂\overline{Ω} × [0, T] to ∂Ω × [0, T] by zero, the corresponding adjoint (\overline{R}′_{j}(f))^{∗} : \overline{Z} → F is then given by (R′_{j}(f))^{∗}E.

7. Laser tomography in the diffusion approximation:

Time harmonic case.

Here the partial differential equation is

\frac{∂v_{j}}{∂t} (x, t) = div(σ(x)∇v_{j}(x, t)) − a(x)v_{j}(x, t) + e^{iωt}q(x − s_{j}) in Ω × IR^{1}
v_{j}(x, t) = 0 , \frac{∂v_{j}}{∂ν} (x, t) = e^{iωt}g_{j}(x) on ∂Ω × IR^{1} .

The modulating frequency ω is fixed, and g_{j} is the measured data for the source s_{j} ∈ ∂Ω.

Putting v_{j}(x, t) = e^{iωt}u_{j}(x) we obtain

0 = div(σ(x)∇u_{j}(x)) − (a(x) + iω)u_{j}(x) + q(x − s_{j}) in Ω ,

(3.29)

u_{j} = 0 , \frac{∂u_{j}}{∂ν} = g_{j} on ∂Ω .

The problem is to recover σ, a from the g_{j}. It is clear that (3.29) is of the form (1.1-2). We have

R_{j} ( \begin{matrix}σ \\ a \end{matrix} ) = \frac{∂u_{j}}{∂ν} − g_{j} ,

R_{j} being considered as an operator from L_{2}(Ω) × L_{2}(Ω) into L_{2}(∂Ω). We

readily compute the operator R′_{j} ( \begin{matrix}σ \\ a \end{matrix} ) : L_{2}(Ω)×L_{2}(Ω) → L_{2}(∂Ω), obtaining

R′_{j} ( \begin{matrix}σ \\ a \end{matrix} ) ( \begin{matrix}η \\ h \end{matrix} ) = \frac{∂w_{j}}{∂ν}

where w_{j} solves

0 = div(η∇u_{j}) + div(σ∇w_{j}) − hu_{j} − (a + iω)w_{j} in Ω
w_{j} = 0 on ∂Ω .

Theorem 3.7 The operator (R′_{j} ( \begin{matrix}σ \\ a \end{matrix} ))^{∗} = L_{2}(∂Ω) → L_{2}(Ω) × L_{2}(Ω) is

given by

(R′_{j} ( \begin{matrix}σ \\ a \end{matrix} ))^{∗} g = ( \begin{matrix}∇\overline{u_{j}} · ∇z \\ \overline{u_{j}}z \end{matrix} )

where z is the solution of

div(σ∇z) − (a − iω)z = 0 in Ω
z = g on ∂Ω

Proof: The proof is done by verification. We compute

(R′_{j} ( \begin{matrix}σ \\ a \end{matrix} ) ( \begin{matrix}η \\ h \end{matrix} ) , g) = ∫_{∂Ω} \frac{∂w_{j}}{∂ν} \overline{g} ds

= ∫_{∂Ω} w_{j} \frac{∂\overline{g}}{∂ν} ds − ∫_{Ω} \{w_{j}div(σ∇\overline{z}) − div(σ∇w_{j})\overline{z}\} dx

by way of Green’s theorem. The integral over ∂Ω vanishes, and the integral over Ω can be rewritten as

∫_{Ω} \{w_{j}(div(σ∇\overline{z}) − (a + iω)\overline{z}) − (div(σ∇w_{j}) − (a + iω)w_{j})\overline{z}\} dx

which, by the differential equations for z and w_{j}, can be written as

− ∫_{Ω} (−div(η∇\overline{u}_{j}) + h\overline{u}_{j})\overline{z}dx .

Hence

(R′_{j} ( \begin{matrix}σ \\ a \end{matrix} ) ( \begin{matrix}η \\ h \end{matrix} ) , g) = ∫_{Ω} (−div(η∇u_{j}) + hu_{j})\overline{z}dx

= ∫_{Ω} (η∇u_{j} · ∇\overline{z} + hu_{j}\overline{z})dx

= (( \begin{matrix}η \\ h \end{matrix} ) , ( \begin{matrix}∇\overline{u}_{j} · ∇z \\ \overline{u}_{j}z \end{matrix} )) .
□

8. Laser tomography in the diffusion approximation:

Time resolved case

In this case the governing equation reads

\frac{∂u_{j}}{∂t} = div(σ(x)∇u_{j}) − a(x)u_{j} + q_{j} in Ω × [0, T],

(3.30)

u_{j}(x, 0) = 0 in Ω, u_{j}(x, t) = 0 on ∂Ω × [0, T].

The source terms q_{j} are assumed to be δ-functions at t = 0 and x = s_{j} ∈ ∂Ω. Now the invers problem is to determine both the diffusion coeffient σ(x) and the absorption coefficient a(x) from the given data

\frac{∂u_{j}}{∂ν} (x, t) = g_{j}(x, t) on ∂Ω × [0, T].

As before this problem is of the form (1.1)-(1.2).

The operator R_{j} : L_{2}(Ω) × L_{2}(Ω) → L_{2,σ}(∂Ω × [0, T]) is given by

R_{j}( \begin{matrix}σ \\ a \end{matrix} ) = \frac{∂u_{j}}{∂ν} − g_{j} on ∂Ω × [0, T]

where u_{j} is the solution of (3.30). According to (2.2) we calculate the linearized operator to be

R′_{j}( \begin{matrix}σ \\ a \end{matrix} )( \begin{matrix}η \\ h \end{matrix} ) = \frac{∂w_{j}}{∂ν} |_{∂Ω×[0,T]}

where w_{j} solves the equation

\frac{∂w_{j}}{∂t} − div(σ(x)∇w_{j}) + a(x)w_{j} = div(η(x)∇u_{j}) − h(x)u_{j} in Ω × [0, T],

(3.31)

w_{j}(x, 0) = 0 in Ω,  w_{j}(x, t) = 0 on ∂Ω × [0, T].

Theorem 3.8 The operator (R′_{j}( \begin{matrix}σ \\ a \end{matrix} ))^{∗} : L_{2,σ}(∂Ω × [0, T]) → L_{2}(Ω) × L_{2}(Ω)

is given by

(R′_{j}( \begin{matrix}σ \\ a \end{matrix} ))^{∗}g(x) = ∫_{0}^{T} ( \begin{matrix}(∇u_{j} · ∇z)(x, t) \\ (u_{j}z)(x, t) \end{matrix} ) dt.

where z(x, t) is the solution of

\frac{∂z}{∂t} + div(σ(x)∇z) − a(x)z = 0 in Ω × [0, T],

(3.32)

z(x, T) = 0 for x ∈ Ω, z(x, t) = g on ∂Ω × [0, T].

Proof: For any z, w we have

∫_{0}^{T} ∫_{Ω} \{(\frac{∂w}{∂t} − div(σ∇w) + aw)z + (\frac{∂z}{∂t} + div(σ∇z) − az)w\} dx dt

= ∫_{Ω} ((wz)(x, T) − (wz)(x, 0)) dx + ∫_{0}^{T} ∫_{∂Ω} σ(\frac{∂z}{∂ν} w − z \frac{∂w}{∂ν}) ds dt.

For w_{j}, z from (3.31)-(3.32) this reduces to

∫_{0}^{T} ∫_{Ω} (div(η∇u_{j}) − hu_{j})z dx dt = − ∫_{0}^{T} ∫_{∂Ω} σgR′_{j}( \begin{matrix}σ \\ a \end{matrix} )( \begin{matrix}η \\ h \end{matrix} ) ds dt.

An integration by parts on the left hand side leads to

∫_{Ω} \{η ∫_{0}^{T}(∇u_{j} · ∇z) dt + h ∫_{0}^{T}(u_{j}z) dt\} dx = (g, R′_{j}( \begin{matrix}σ \\ a \end{matrix} )( \begin{matrix}η \\ h \end{matrix} ))_{L_{2,σ}} .

This proofs the theorem. □

4 Convergence

We now prove convergence of the iterative method for linear equations

R_{j}f = g_{j} ,  j = 1, …, p  (4.1)

where R_{j} : F → Z_{j} are linear bounded operators. To avoid purely technical difficulties we assume F and Z_{j} to be of finite dimension. The iteration reads

(i) f_{0} = f^{0} ,
(ii) For j = 1, …, p

(4.2)

f_{j} = f_{j−1} + ωh_{j} , h_{j} = −R_{j}^{∗}C_{j}^{−1}(g_{j} − R_{j}f_{j−1}) ,

(iii) f^{1} = f_{p} .

We note in passing that for dim (Z_{j}) = 1 this is just the Kaczmarz method.

Theorem 4.1 Let C_{j} be positive definite and C_{j} ≥ R_{j}R_{j}^{∗}. Then, (4.2) converges. If (4.1) is consistent and

f^{0} ∈ ∑_{j=1}^{p} range(R_{j}^{∗}) ,  (4.3)

then (4.2) converges to the solution of (4.1) with least norm.

Proof: For C_{j} = R_{j}R_{j}^{∗} this is just Theorem 3.9 of [3]. Our proof is a slight extension of this result.

To begin with we write f_{j} in the form

f_{j} = f_{0} + ∑_{k=1}^{j} R_{k}^{∗}u_{k}  (4.4)

with certain vectors u_{k}. Inserting this into the recursion (4.2) yields

f_{0} + ∑_{k=1}^{j} R_{k}^{∗}u_{k} = f_{0} + ∑_{k=1}^{j−1} R_{k}^{∗}u_{k} + ωR_{j}^{∗}C_{j}^{−1}(g_{j} − R_{j}f_{0} − ∑_{k=1}^{j−1} R_{j}R_{k}^{∗}u_{k}) .

We apply R_{j}. After obvious manipulations we obtain

R_{j}R_{j}^{∗}u_{j} = ωR_{j}R_{j}^{∗}C_{j}^{−1}(g_{j} − R_{j}f_{0} − ∑_{k=1}^{j−1} R_{j}R_{k}^{∗}u_{k}) .

Hence,

C_{j}u_{j} = ω(g_{j} − R_{j}f_{0} − ∑_{k=1}^{j−1} R_{j}R_{k}^{∗}u_{k}) .  (4.5)

Now let A be the operator (R_{j}R_{k}^{∗})_{j,k=1,...,p} and R the map (R_{j})_{j=1,...,p}, i.e. A = RR^{∗}. Let A = L + L^{∗} + D where D is the (block) diagonal and L the (block) lower left part of A. Let the vectors u, g consist of the components u_{j}, g_{j}, resp. and let C be the (block) diagonal C_{1}, …, C_{p}. Then (4.5) reads

Cu = ω(g − Rf_{0} − Lu)

or

u = ω(C + ωL)^{−1}(g − Rf_{0}) .  (4.6)

Using (4.4) with j = p we see that

f^{1} = f^{0} + R^{∗}u .  (4.7)

Now assume that f^{0} ∈ ∑_{j=1}^{p} range(R_{j}^{∗}) = range(R_{j}^{∗}). Then, f^{k} = R^{∗}u^{k} for all

iterates, and (4.6)-(4.7) yield

R^{∗}u^{1} = R^{∗}u^{0} + ωR^{∗}(C + ωL)^{−1}(g − RR^{∗}u^{0})

or

u^{1} = u^{0} + ω(C + ωL)^{−1}(g − RR^{∗}u^{0}) .

Hence

u^{1} = C(ω)u^{0} + c(ω) , C(ω) = I − ω(C + ωL)^{−1}RR^{∗} .  (4.8)

We show that all the eigenvalues of C(ω) are either 1 or less than 1 in absolute value if C_{j} ≥ R_{j}R_{j}^{∗}, j = 1, …, p. Let λ be such an eigenvalue and u the corresponding eigenvector, i.e.

(I − ω(C + ωL)^{−1}RR^{∗})u = λu

or

(C + ωL − ωRR^{∗})u = λ(C + ωL)u .

Using RR^{∗} = D + L + L^{∗} yields

(C − ωD − ωL^{∗})u = λ(C + ωL)u .  (4.9)

Now let (u, Du) = 1 and (Lu, u) = a = α + iβ, γ = (Cu, u) with α, β, γ real. (4.9) immediately yields

γ − ω − ω\overline{a} = λ(γ + ωa)

hence

λ = \frac{γ − ω − ωα + iωβ}{γ + ωα + iωβ} .

Since RR^{∗} ≥ 0 we have α ≥ −1/2, and since C_{j} ≥ R_{j}R_{j}^{∗} we have γ ≥ 1. For α = −1/2 we have λ = 1. For α > −1/2 we have |γ − ω − ωα| < |γ + ωα| provided that ω < 2γ. Thus |λ| < 1 for α > −1/2 and 0 < ω < 2.

Returning to (4.8) we see that the sequence u^{k} converges except in the eigenspace of C(ω) for the eigenvalue 1. But this eigenspace is identical to the null space of R^{∗}. This does not prevent the f^{k} = R^{∗}u^{k} from converging.

□

References

[1] Bui, H.D.: Inverse Problems in the Mechanics of Materials. CRC Press 1994.

[2] Gutman, S. - Klibanov, M.: Regularized quasi-Newton method for inverse scattering problems. Math. Comput. Modelling 18, 5-31 (1993). Oxford: Pergamon.

[3] Natterer, F.: The Mathematics of Computerized Tomography. Wiley - Teubner 1986.

[4] Natterer, F. - Wübbeling, F.: A propagation-backpropagation method for uzltrasound tomography. Inverse Problems 11, 1225-1232 (1995). IOP Publishing Ltd.