Jon Claerbout

Mathematical physics presents us with many partial differential equations, such as the heat-flow equation, the acoustic, seismic, or electromagnetic wave equation. Such equations are frequently simulated on computers by finite-difference approximations. In applied science we often need to find the inverse operator to such a simulation to image the source field responsible for the observed data. Many inversion methods require an implementation of the adjoint (matrix transpose) of the linear simulator process.

In this short extract of a more detailed paper (http://sepwww.stanford.edu/redoc/), we demonstrate a straightforward programming scheme to implement the adjoint to a finite-difference simulation of a linear operator.

The basic idea is that the adjoint of a row vector is a column vector. So the adjoint of the equation y=ax₁+bx₂ is two equations, $\hat x_1 = ay$ and $\hat x_2 = by$.Since we often combine operators by element-wise addition, we often code them in the form $y\leftarrow y+ax_1+bx_2$.A convenient notation comes from the C computer language where such assignments are written as y += a*x1+b*x2. In this case the adjoint can be implemented as x1 += a*y and x2 += b*y.

The one-dimensional heat-flow equation can be used to represent viscosity and is the simplest place to start:  

$${\partial q \over \partial t} \quad =\quad\sigma {\partial^2 q \over \partial x^2} \ +\ s(x,t)$$ (1)

where q is temperature, $\sigma$ is heat conductivity divided by heat capacity and s is a possible source of heat. Discretizing equation 1 with the usual explicit finite-difference operator gives  

$$q_{t+1}^x \quad +=\quad q_t^x \ +\ \sigma ( q_t^{x-1}-2q_t^x+q_t^{x+1} ) \ +\ s_t^x$$ (2)

where t takes integer values, and $\sigma$ is like the $\sigma$ in the differential equation but includes $\Delta t/\Delta x^2$.Ignoring the source s_t^x, the adjoint of 2 is

$$\begin{array} {lcr} q_t^x &\quad +=\quad& q_{t+1}^x \\ q_t^{... ...+1}^x \\ q_t^{x+1} &\quad +=\quad& \sigma q_{t+1}^x \end{array}$$ (3)

A step of the heat-flow equation transforms from q_t to q_t+1. The adjoint operation (transpose matrix) transforms q_t+1 to $\hat q_t$ (which in general is different from q_t), thus stepping time backward. Temperature smooths as time runs forward. If we had the inverse matrix we could step backward in time and would find the rougher temperature q_t. Running backward in time with the adjoint matrix $\hat q_t$,we find temperature getting smoother. For heat-flow, this is the main difference between the adjoint and the inverse.

The heat-flow operator H marches one step in time. Next, using a similar method, we should find finite-difference operators for the pressure operator P and the velocity operator V. These details are found in the electronic document referred to above. Finally, we represent the viscous scalar wave equation by the product HVP, so its adjoint is simply P'V'H'.

Any time we implement an adjoint operator, we need to confirm that we compute the correct adjoint of a forward operator. Mathematics demands that y'(Ax)=(A'y)'x for a linear operator A, its adjoint A', and ``all possible functions'' x and y. In practice, it suffices to test y'(Ax)=(A'y)'x using any set of random numbers for x and y and our program for the operators A and A'. The electronic version of this document shows that for random inputs of size 100³, the quantities y'(Ax) and (A'y)'x are equal to the expected six digits of numerical accuracy. The result of this dot-product test is stored in the file dot.txt.

Since finite-difference representations of differential equations can be unstable, a further necessary test of the viscous wave equation program is the sample output below. I put random point sources in the frog pond on a $100 \times 100$ matrix.

FIG: A rectangular pond after a frog's random
hops disturbed its surface. The subroutines that compute
this figure are included in the electronic version of
this document.