On the solution of the linear systems which evolve from large scale inverse problems

Eldad Haber

Dept of Computer Science, the University of British Columbia Vancouver Canada

In this talk we consider an inverse problem where the forward problem is governed by a (discrete) Partial Differential Equation of the form

(1) A(m) u = f

where A(m) is a matrix which depends on the "model" m, u is some physical field and f is the source and boundary conditions. Such problems include the DC resistivity (or impedance tomography), inverse Maxwell's equations, hydrology and many other problems. The goal of the inverse problem is to recover the model m based on the measurements of the field u.

In the inverse problem we consider recovering m from measurements of u in a discrete number of locations. The problem is usually solved by minimizing an objective function of the form

(2) min phi = ||A(m)^{-1}f - b ||^2 + beta || Wm||^2

where b is the measured data beta is a penalty parameter and W has the a priori information.

The optimization problem leads to the solution of the Gauss-Newton (or quasi-Newton) step at each iteration for the perturbation dm.

(3) Ldm = (J(m)'J(m) + beta W'W) dm = -r(m)

Where the matrix J is a dense sensitivity matrix

(4) J(m) = db(m)/dm

and r(m) is the gradient of the objective function in (2).

The solution of this system is complicated because it involves large dense matrices. Furthermore, preconditioning the iteration is not simple due to the different nature of the sensitivity matrix J (which can be viewed as an integral operator) and the weighting matrix W (which is usually a difference matrix).

In this talk we present a new approach for the solution of such systems. First, we show that this system is equivalent to the solution of a sparse KKT matrix. Second, we use the KKT matrix in order to construct a sparse approximation approximate to J based on the reduced Hessian method applied to the sparse KKT system. Thirdly, we show that using the sparse approximation we can use standard preconditioners in order to solve the approximate system.

We demonstrate the efficiency of the method on two model problems from electromagnetics.