paper.tex

国外免费地震资料处理软件包

%\def\SEPCLASSLIB{../../../../sepclasslib}
\def\CAKEDIR{.}

\title{Adjoint operators}
\author{Jon Claerbout}
\label{paper:conj}
\maketitle

A great many of the calculations
we do in science and engineering
are really matrix multiplication in disguise.
The first goal of this chapter is to unmask the disguise
by showing many examples.
Second, we see how the 
\bx{adjoint} operator (matrix transpose)
back-projects information from data to the underlying model.

\par
Geophysical modeling calculations
generally use linear operators that predict data from models.
Our usual task is to find the inverse of these calculations;
i.e., to find models (or make maps) from the data.
Logically, the adjoint is the first step
and a part of all subsequent steps in this \bx{inversion} process.
Surprisingly, in practice the adjoint sometimes does a better job
than the inverse!
This is because the adjoint operator tolerates imperfections
in the data and does not demand that the data provide full information.

\par
Using the methods of this chapter,
you will find that
once you grasp the relationship between operators in general
and their adjoints,
you can obtain the adjoint just
as soon as you have learned how to code
the modeling operator.

\par
If you will permit me a poet's license with words,
I will offer you the following table
of \bx{operator}s and their \bx{adjoint}s:

\begin{tabular}{p{2em}lp{1em}l}
&\bx{matrix multiply}		&&conjugate-transpose matrix multiply \\
&convolve 			&&crosscorrelate	\\
&truncate			&&zero pad	\\
&replicate, scatter, spray	&&sum or stack	\\
&spray into neighborhood	&&sum in bins	\\
&derivative (slope) 		&&negative derivative	\\
&causal integration 		&&anticausal integration	\\
&add functions			&&do integrals	\\
&assignment statements		&&added terms	\\
&plane-wave superposition	&&slant stack / beam form	\\
&superpose on a curve		&&sum along a curve \\
&stretch			&&squeeze	\\
&upward continue		&&downward continue \\
&hyperbolic modeling	 	&&normal moveout and CDP stack	\\
&diffraction modeling	 	&&imaging by migration	\\
&ray tracing			&&\bx{tomography}
\end{tabular}
\par
The left column above is often called ``\bx{modeling},''
and the adjoint operators on the right are often
used in ``data \bx{processing}.''

\par
The adjoint operator is sometimes called
the ``\bx{back projection}'' operator
because information propagated in one direction (earth to data) is projected
backward (data to earth model).
For complex-valued operators,
the transpose goes together with a complex conjugate.
In \bx{Fourier analysis}, taking the complex conjugate
of $\exp(i\omega t)$ reverses the sense of time.
With more poetic license, I say that adjoint operators
{\it undo}
the time and phase shifts of modeling operators.
The inverse operator does this too,
but it also divides out the color.
For example, when linear interpolation is done,
then high frequencies are smoothed out,
so inverse interpolation must restore them.
You can imagine the possibilities for noise amplification.
That is why adjoints are safer than inverses.

\par
Later in this chapter we relate
adjoint operators to inverse operators.
Although inverse operators are more well known than adjoint operators,
the inverse is built upon the adjoint
so the adjoint is a logical place to start.
Also, computing the inverse is a complicated process
fraught with pitfalls whereas
the computation of the adjoint is easy.
It's a natural companion to the operator itself.

%Often, the adjoint is all you need for imaging
%and when it isn't
%The adjoint is often the best place to end too,
%because inverses of giant operators (and these are giants)
%are notoriously ill-behaved.

%\par
%When the adjoint operator is
%{\it not}
%an adequate approximation to the inverse,
%then you apply the techniques of fitting and optimization
%that are not within the scope of this book.
%The computed adjoint should be, and generally can be,
%exact (within machine precision).

\par
Much later in this chapter is a formal definition of adjoint operator.
Throughout the chapter we handle an adjoint operator as a matrix transpose,
but we hardly ever write down any matrices or their transposes.
Instead, we always prepare two subroutines,
one that performs 
$\bold y =\bold A \bold x$
and another that performs
$\tilde{\bold x} =\bold A' \bold y$.
So we need a test that the two subroutines
really embody the essential aspects of matrix transposition.
Although the test is an elegant and useful test
and is itself a fundamental definition,
curiously,
that definition
does not help construct
adjoint operators,
so we postpone a formal definition of adjoint
until after we have seen many examples.

\section{FAMILIAR OPERATORS}

\par
\sx{matrix multiply}
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 \bx{pseudocode} does matrix multiplication
$\bold y=\bold B\bold x$ and multiplication by the transpose
$\tilde \bold x = \bold B' \bold y$:

\par
\vbox{\begin{tabbing}
indent \= atechars \= atechars \=  atechars \= \kill
\>if operator itself	\\
\>	\>then erase y	\\
\>if adjoint	\\
\>	\>then erase x	\\
\>do iy = 1, ny \{	\\
\>do ix = 1, nx \{	\\
\>	\>if operator itself	\\
\>	\>	\>y(iy) = y(iy) + b(iy,ix) $\times$ x(ix)	\\
\>	\>if adjoint	\\
\>	\>	\>x(ix) = x(ix) + b(iy,ix) $\times$ y(iy)	\\
\>	\>\}\}
\end{tabbing} }
\par

\noindent
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 operation
and its adjoint.

\par
A formal subroutine
%\footnote{
%	The programming language used in this book is Ratfor,
%	a dialect of Fortran.
%	For more details, see Appendix A.}
for \bx{matrix multiply} and its adjoint is found below.
The first step is a subroutine, \texttt{adjnull()},
for optionally erasing the output.
With the option {\tt add=true}, results accumulate like {\tt y=y+B*x}.
\moddex{adjnull}{erase output}{26}{48}{filt/lib}
The subroutine \texttt{matmult()}
for matrix multiply and its adjoint
exhibits the style that we will use repeatedly.
\opdex{matmult}{matrix multiply}{33}{45}{user/fomels}

% NEW
\par
Sometimes a matrix operator reduces to a simple row or a column.

\par\noindent
A {\bf row} \quad\ \ is a summation operation.
\par\noindent
A {\bf column}       is an impulse response.
\vspace{.2in}
\par\noindent
If the inner loop of a matrix multiply ranges within a
\par\noindent
{\bf row,} \quad\ \ the operator is called {\it sum} or {\it pull}.
\par\noindent
{\bf column,}       the operator is called {\it spray} or {\it push}.

\par
A basic aspect of adjointness is that the
adjoint of a row matrix operator is a column matrix operator.
For example,
the row operator $[a,b]$
\begin{equation}
y \eq
\left[ \ a \ b \ \right] 
\left[
\begin{array}{l}
	x_1 \\
	x_2
\end{array}
\right] 
\eq
a x_1 + b x_2
\end{equation}
has an adjoint that is two assignments:
\begin{equation}
	\left[
	\begin{array}{l}
		\hat x_1 \\
		\hat x_2
	\end{array}
	\right]
	\eq
	\left[
	\begin{array}{l}
		a \\
		b
	\end{array}
	\right]
	\ y
\end{equation}
\par
\boxit{
The adjoint of a sum of $N$ terms
is a collection of $N$ assignments.
}

\subsection{Adjoint derivative}
Given a sam\-pled sig\-nal,
its time \bx{derivative} can be esti\-mated
by con\-vo\-lu\-tion with the fil\-ter $(1,-1)/\Delta t$,
expressed as the matrix-multiply below:
\begin{equation}
\left[ \begin{array}{c}
	y_1 \\
	y_2 \\
	y_3 \\
	y_4 \\
	y_5 \\
	y_6
	\end{array} \right]
\eq
\left[ \begin{array}{cccccc}
	-1& 1& .& .& .& . \\
	 .&-1& 1& .& .& . \\
	 .& .&-1& 1& .& . \\
	 .& .& .&-1& 1& . \\
	 .& .& .& .&-1& 1 \\
	 .& .& .& .& .& 0
	\end{array} \right] \ 
\left[ \begin{array}{c}
	x_1 \\
	x_2 \\
	x_3 \\
	x_4 \\
	x_5 \\
	x_6
	\end{array} \right]
 \label{eqn:ruff1}
\end{equation}
Technically the output should be {\tt n-1} points long,
but I appended a zero row,
a small loss of logical purity,
so that the size of the output vector will match that of the input.
This is a convenience for plotting
and for simplifying the assembly of other operators building on this one.
\par
The \bx{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.
\par
A complicated way to think about the adjoint of equation
(\ref{eqn:ruff1}) is to note that it is the negative of the derivative
and that something must be done about the ends.
A simpler way to think about it
is to apply the idea that the adjoint of a sum of $N$ terms
is a collection of $N$ assignments.
This is done in subroutine \texttt{igrad1()},
which implements equation~(\ref{eqn:ruff1})
and its adjoint.
\moddex{igrad1}{first difference}{25}{40}{filt/proc}
\par
\noindent
Notice that the do loop in the code
covers all the outputs for the operator itself,
and that in the adjoint operation it gathers all the inputs.
This is natural because in switching from operator
to adjoint, the outputs switch to inputs.
\par
As you look at the code,
think about matrix elements being $+1$ or $-1$ and
think about the forward operator
``pulling'' a sum into {\tt yy(i)}, and
think about the adjoint operator
``pushing'' or ``spraying'' the impulse {\tt yy(i)} back into {\tt xx()}.

\par
You might notice that you can simplify the program
by merging the ``erase output'' activity with the calculation itself.
We will not do this optimization however because in many applications
we do not want to include the ``erase output'' activity.
This often happens when we build complicated operators from simpler ones.


\subsection{Zero padding is the transpose of truncation}
Surrounding a dataset by zeros
(\bx{zero pad}ding)
is adjoint to throwing away the extended data
(\bx{truncation}).
Let us see why this is so.
Set a signal in a vector $\bold x$, and
then to make a longer vector $\bold y$,
add some zeros at the end of $\bold x$.
This zero padding can be regarded as the matrix multiplication
\begin{equation}
\bold y\eq
 \left[ 
  \begin{array}{c}
   \bold I \\ 
   \bold 0
  \end{array}
 \right] 
 \ 
 \bold x
\end{equation}
The matrix is simply an identity matrix $\bold I$
above a zero matrix $\bold 0$.
To find the transpose to zero-padding, we now transpose the matrix
and do another matrix multiply:
\begin{equation}
\tilde {\bold x} \eq
 \left[ 
  \begin{array}{cc}
   \bold I & \bold 0
  \end{array}
 \right] 
\ 
\bold y
\end{equation}
So the transpose operation to zero padding data
is simply {\it truncating} the data back to its original length.
Subroutine \texttt{zpad1()} below
pads zeros on both ends of its input.
Subroutines for two- and three-dimensional padding are in the
library named {\tt zpad2()} and {\tt zpad3()}.
\opdex{zpad1}{zero pad 1-D}{21}{32}{user/gee}


\subsection{Adjoints of products are reverse-ordered products of adjoints}
Here we examine an example of the general idea that
adjoints of products are reverse-ordered products of adjoints.
For this example we use the Fourier transformation.
No details of \bx{Fourier transformation} are given here
and we merely use it as an example of a square matrix $\bold F$.
We denote the complex-conjugate transpose (or \bx{adjoint}) matrix
with a prime,
i.e.,~$\bold F'$.
The adjoint arises naturally whenever we consider energy.
The statement that Fourier transforms conserve energy is 
$\bold y'\bold y=\bold x'\bold x$ where $\bold y= \bold F \bold x$.
Substituting gives $\bold F'\, \bold F = \bold I$, which shows that
the inverse matrix to Fourier transform
happens to be the complex conjugate of the transpose of $\bold F$.
\par
With Fourier transforms,
\bx{zero pad}ding and \bx{truncation} are especially prevalent.
Most subroutines transform a dataset of length of $2^n$,
whereas dataset lengths are often of length $m \times 100$.
The practical approach is therefore to pad given data with zeros.
Padding followed by Fourier transformation $\bold F$
can be expressed in matrix algebra as
\begin{equation}
{\rm Program} \eq
\bold F \ 
 \left[ 
  \begin{array}{c}
   \bold I \\ 
   \bold 0
  \end{array}
 \right] 
\end{equation}
According to matrix algebra, the transpose of a product,
say $\bold A \bold B = \bold C$,
is the product $\bold C' = \bold B' \bold A'$ in reverse order.
So the adjoint subroutine is given by
\begin{equation}
{\rm Program'} \eq
 \left[ 
  \begin{array}{cc}
   \bold I & \bold 0
  \end{array}
 \right] 
\
\bold F'
\end{equation}
Thus the adjoint subroutine
{\it truncates} the data {\it after} the inverse Fourier transform.
This concrete example illustrates that common sense often represents
the mathematical abstraction
that adjoints of products are reverse-ordered products of adjoints.
It is also nice to see a formal mathematical notation
for a practical necessity.
Making an approximation need not lead to collapse of all precise analysis.


\subsection{Nearest-neighbor coordinates}
\sx{nearest neighbor coordinates}
In describing physical processes,
we often either specify models as values given on a uniform mesh
or we record data on a uniform mesh.
Typically we have
a function $f$ of time $t$ or depth $z$
and we represent it by {\tt f(iz)}
corresponding to $f(z_i)$ for $i=1,2,3,\ldots, n_z$
where $z_i = z_0+ (i-1)\Delta z$.
We sometimes need to handle depth as
an integer counting variable $i$
and we sometimes need to handle it as
a floating-point variable $z$.
Conversion from the counting variable to the floating-point variable
is exact and is often seen in a computer idiom
such as either of \begin{verbatim}
            do iz= 1, nz {   z = z0 + (iz-1) * dz 
            do i3= 1, n3 {  x3 = o3 + (i3-1) * d3
\end{verbatim}
%{\tt
%   \begin{tabbing}  indent \= \kill
%	\> do iz= 1, nz \{   z = z0 + (iz-1) * dz  \\
%	\> do i3= 1, n3 \{  x3 = o3 + (i3-1) * d3