next up previous [pdf]

Next: Transient convolution Up: FAMILIAR OPERATORS Previous: FAMILIAR OPERATORS

Adjoint derivative

In numerical analysis, we represent the derivative of a time function by a finite difference. This subtracts neighboring time points and divides by the sample interval $\Delta t$. Finite difference amounts to convolution with the filter $(1,-1)/\Delta t$. Omitting the $\Delta t$, we express this concept as:
\begin{displaymath}
\left[ \begin{array}{c}
y_1 \\
y_2 \\
y_3 \\
y_4 \\
...
...x_2 \\
x_3 \\
x_4 \\
x_5 \\
x_6
\end{array} \right]
\end{displaymath} (3)

The filter is seen in any column in the middle of the matrix, namely $(1,-1)$. In the transposed matrix, the filter-impulse response is time-reversed to $(-1,1)$. Therefore, mathematically, we can say that the adjoint of the time derivative operation is the negative time derivative. Likewise, in the Fourier domain, the complex conjugate of $-i\omega$ is $i\omega$. We can also speak of the adjoint of the boundary conditions: we might say that the adjoint of ``no boundary condition'' is a ``specified value'' boundary condition. The last row in equation (3) is optional. It may seem unnatural to append a null row, but it can be a small convenience (when plotting) to have the input and output be the same size.

Equation (3) is implemented by the code in module igrad1 that does the operator itself (the forward operator) and its adjoint.

api/c/igrad1.c
void  sf_igrad1_lop(bool adj, bool add, 
		    int nx, int ny, float *xx, float *yy)
/*< linear operator >*/
{
    int i;

    sf_adjnull(adj,add,nx,ny,xx,yy);
    for (i=0; i < nx-1; i++) {
        if (adj) {
	    xx[i+1] += yy[i];
	    xx[i]   -= yy[i];
	} else {
	    yy[i] += xx[i+1] - xx[i];
	}
    }
}

The adjoint code may seem strange. It might seem more natural to code the adjoint to be the negative of the operator itself; and then, make the special adjustments for the boundaries. The code given, however, is correct and requires no adjustments at the ends. To see why, notice for each value of i, the operator itself handles one row of equation (3), while for each i, the adjoint handles one column. That is why coding the adjoint in this way does not require any special work on the ends. The present method of coding reminds us that the adjoint of a sum of $N$ terms is a collection of $N$ assignments. Think of the meaning of $y_i=y_i+a_{i,j}x_j$ for any particular $i$ and $j$. The adjoint simply accumulates that same value of $a_{i,j}$ going the other direction $x_j=x_j+a_{i,j}y_i$.

where we see that each component of the matrix is handled both by the operator and the adjoint. Think about the forward operator ``pulling'' a sum into yy(i), and think about the adjoint operator ``pushing'' or ``spraying'' the impulse yy(i) back into xx().

Figure 1 illustrates the use of module igrad1 for each north-south line of a topographic map. We observe that the gradient gives an impression of illumination from a low sun angle.

stangrad
stangrad
Figure 1.
Topography near Stanford (top) southward slope (bottom).
[pdf] [png] [scons]

To apply igrad1 along the 1-axis for each point on the 2-axis of a two-dimensional map, we use the loop

for (iy=0; iy < ny; iy++) 
     igrad1_lop(adj, add, nx, nx, map[iy], ruf[iy]);


next up previous [pdf]

Next: Transient convolution Up: FAMILIAR OPERATORS Previous: FAMILIAR OPERATORS

2014-09-27