Adjoint derivative

In numerical analysis we represent the derivative a time function by a finite difference. We do this by subtracting each two neighboring time points and then dividing by the sample interval $\Delta t$.This amounts to convolution with the filter $(1,-1)/\Delta t$.Omitting the $\Delta t$ we express this concept as:  

$$\left[ \begin{array} {c} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_... ..._1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \end{array} \right]$$ (3)

The filter impulse response 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). So, mathematically, we can say that the adjoint of the time derivative operation is the negative time derivative. This corresponds also to the fact that 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 which does the operator itself (the forward operator) and its adjoint. igrad1first difference

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's 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.

The Ratfor90 dialect of Fortran allows us to write the inner code of the igrad1 module more simply and symmetrically using the syntax of C, C++, and Java where expressions like a=a+b can be written more tersely as a+=b. With this, the heart of module igrad1 becomes

if( adj) {   xx(i+1) += yy(i)
             xx(i)   -= yy(i)
           }
else {       yy(i)   += xx(i+1)
             yy(i)   -= xx(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 stangrad90 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.

 

stangrad90
Figure 1 Topography near Stanford (top) southward slope (bottom).

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

do iy=1,ny
      stat = igrad1_lop( adj, add, map(:,iy), ruf(:,iy))

On the other hand, to see the east-west gradient, we use the loop

do ix=1,nx 
      stat = igrad1_lop( adj, add, map(ix,:), ruf(ix,:))