Basic operators and adjoints

Programming linear operators

The operation $y_i = \sum_j b_{ij} x_j$ is the multiplication of a matrix $\bold B$ by a vector $\bold x$. The adjoint operation is $\tilde x_j = \sum_i b_{ij} y_i$. The operation adjoint to multiplication by a matrix is multiplication by the transposed matrix (unless the matrix has complex elements, in which case we need the complex-conjugated transpose). The following pseudocode does matrix multiplication $\bold y=\bold B\bold x$ and multiplication by the transpose $\tilde{\bold x} = \bold{B}' \bold{y}$:

		if adjoint    


then erase x  


if operator itself    


then erase y  


do iy = 1, ny {      


do ix = 1, nx {      


if adjoint    


x(ix) = x(ix) + b(iy,ix) $\times$ y(iy)       


if operator itself    


y(iy) = y(iy) + b(iy,ix) $\times$ x(ix)       

}}

Notice that the ``bottom line'' in the program is that $x$ and $y$ are simply interchanged. The above example is a prototype of many to follow, so observe carefully the similarities and differences between the adjoint and the operator itself.

Next we restate the matrix-multiply pseudo code in real code.

The module matmult for matrix multiply and its adjoint exhibits the style that we will use repeatedly. At last count there were 53 such routines (operator with adjoint) in this book alone.

user/fomels/matmult.c

void matmult_lop (bool adj, bool add, 
		  int nx, int ny, float* x, float*y) 
/*< linear operator >*/
{
    int ix, iy;
    sf_adjnull (adj, add, nx, ny, x, y);
    for (ix = 0; ix < nx; ix++) {
	for (iy = 0; iy < ny; iy++) {
	    if (adj) x[ix] += B[iy][ix] * y[iy];
	    else     y[iy] += B[iy][ix] * x[ix];
	}
    }
}

We have a module with two entries, one named _init for the physical scientist who defines the physical problem by defining the matrix, and another named _lop for the least-squares problem solver, the computer scientist who will not be interested in how we specify $\bold B$, but who will be iteratively computing $\bold B\bold x$ and $\bold B' \bold y$ to optimize the model fitting.

To use matmult, two calls must be made, the first one

           matmult_init( bb);

is done by the physical scientist after he or she has prepared the matrix. Most later calls will be done by numerical analysts in solving code like in Chapter . These calls look like

            matmult_lop( adj, add, nx, ny, x, y);

where adj is the logical variable saying whether we desire the adjoint or the operator itself, and where add is a logical variable saying whether we want to accumulate like $\bold y \leftarrow \bold y+\bold B\bold x$ or whether we want to erase first and thus do $\bold y \leftarrow \bold B\bold x$.

We split operators into two independent processes, the first is used for geophysical set up while the second is invoked by mathematical library code (introduced in the next chapter) to find the model that best fits the data. Here is why we do so. It is important that the math code contain nothing about the geophysical particulars. This enables us to use the same math code on many different geophysical problems. This concept of ``information hiding'' arrived late in human understanding of what is desireable in a computer language. Subroutines and functions are the way that new programs use old ones. Object modules are the way that old programs (math solvers) are able to use new ones (geophysical operators).

Basic operators and adjoints