paper.tex

来自「国外免费地震资料处理软件包」· TEX 代码 · 共 1,940 行 · 第 1/5 页

% copyright (c) 2001 Jon Claerbout

\title{Basic operators and adjoints}
\author{Jon Claerbout}
\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 images) 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 neighborhoods      &&sum within 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 curves              &&sum along a curve \\
&stretch                       &&squeeze       \\
&scalar field gradient         &&negative of vector field divergence       \\
&upward continue               &&downward continue \\
&diffraction modeling          &&imaging by migration  \\
&hyperbola modeling            &&CDP stacking  \\
&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
When the adjoint operator is
{\em  not}
an adequate approximation to the inverse,
then you apply the techniques of fitting and optimization
explained in Chapter \ref{lsq/paper:lsq}.
These techniques require iterative use of the
modeling operator and its adjoint.

\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).
Using complex-valued operators,
the transpose and complex conjugate go together;
and 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
{\em  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.
But nature determines in each application what is the best operator to use,
and whether to stop after the adjoint,
to go the whole way to the inverse,
or to stop partway.

\par 
The operators and adjoints above transform vectors to other vectors.
They also transform data planes to model planes, volumes, etc.
A mathematical operator transforms an ``abstract vector'' which
might be packed full of volumes of information like television
signals (time series) can pack together a movie, a sequence of frames.
We can always think of the operator as being a matrix
but the matrix can be truly huge (and nearly empty).
When the vectors transformed by the matrices are large like
geophysical data set sizes
then the matrix sizes are ``large squared,''
far too big for computers.
Thus although we can always think of an operator as a matrix,
in practice, we handle an operator differently.
Each practical application requires the practitioner to
prepare two computer programs.
One performs the matrix multiply
$\bold y =\bold A \bold x$
and another multiplys by the transpose
$\tilde{\bold x} =\bold A' \bold y$
(without ever having the matrix itself in memory).
It is always easy to transpose a matrix.
It is less easy to take a computer program that does
$\bold y =\bold A \bold x$
and convert it to another to do
$\tilde{\bold x} =\bold A' \bold y$.
In this chapter are many examples of increasing complexity.
At the end of the chapter we will see a test for any program pair
to see whether the operators $\bold A$ and $\bold A'$ are mutually adjoint
as they should be.
Doing the job correctly (coding adjoints without making approximations)
will reward us later when we tackle model and image estimation problems.

%\par
%We will see that computation of the \bx{adjoint}
%is a straightforward adjunct to the operator computation itself,
%and the computed adjoint should be, and generally can be,
%exact (within machine precision).
%We will see
%if the application's operator is computed in an approximate way,
%that it is natural and best that the adjoint
%be computed with adjoint approximations.
%
%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 routines,
%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 routines
%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
%%helps us not one bit in constructing
%does not help construct
%adjoint operators,
%so we postpone a formal definition of adjoint
%until after we have seen many examples.


\subsection{Programming linear operators}
\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\noindent
\vbox{\begin{tabbing}
indent \= atechars \= atechars \=  atechars \= \kill
\>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)       \\
\>      \>\}\}
\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 adjoint and the operator itself.

\par
Next we restate the matrix-multiply pseudo code in real code.
\begin{comment}
in a language called \bx{Loptran}\footnote{
        The programming language, Loptran,
        is based on a dialect of Fortran called Ratfor.
        For more details, see Appendix A.},
a language designed for
exposition and research
in model fitting and optimization in physical sciences.
\end{comment}

The module \texttt{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.
\opdex{matmult}{matrix multiply}{33}{45}{user/fomels}        
\begin{comment}
Notice that the module \texttt{matmult}
does not explicitly erase its output before it begins,
as does the psuedo code.
That is because Loptran will always erase for you
the space required for the operator's output.
Loptran also defines a logical variable {\tt adj}
for you to distinguish your computation of the adjoint
{\tt x=x+B'*y}
from the forward operation
{\tt y=y+B*x}. 
In computerese, the two lines beginning \#\% are macro expansions
that take compact bits of information which expand
into the verbose boilerplate that Fortran requires.
Loptran is Fortran with these macro expansions.
You can always see how they expand by looking at
\texttt{http://sep.stanford.edu/sep/prof/gee/Lib/}.

\par
What is new in Fortran 90, and will be a big help to us,
is that instead of a subroutine with a single entry.
\end{comment}

We have a module with two entries,
one named {\tt \_init} 
for the physical scientist who defines the physical problem by
defining the matrix, and
another named {\tt \_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.
\begin{comment}
The lines beginning with {\tt {\#}{\%}} are expanded by Loptran into
more verbose and distracting Fortran 90 code.
The second line in the module \texttt{matmult},
however,
is pure Fortran syntax saying that
{\tt bb} is a pointer to a real-valued matrix.
\end{comment}
\par
To use \texttt{matmult}, two calls must be made,
the first one
\begin{verbatim}
           matmult_init( bb);
\end{verbatim}
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 \ref{lsq/paper:lsq}.
These calls look like
\begin{verbatim}
            matmult_lop( adj, add, nx, ny, x, y);
\end{verbatim}
where {\tt adj} is the logical variable saying whether we desire
the adjoint or the operator itself,
and where {\tt 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$.
\begin{comment}
The return value {\tt stat} is an integer parameter,
mostly useless (unless you want to use it for error codes). 
\par
Operator initialization often allocates memory.
To release this memory, you can {\tt call matmult\_close()}
although in this case nothing really happens.
\end{comment}
\par
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.
%in other words, to find the zero-offset trace which sprays
%out to make the best approximations to all the other traces.
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.
\begin{comment}
This feature alone is valuable enough to warrant
upgrading from Fortran 77 to Fortran 90, and likewise from C to C++.
\end{comment}
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).


\section{FAMILIAR OPERATORS}
The simplest and most fundamental linear operators
arise when 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 {\em  sum} or {\em  pull}.
\par\noindent
{\bf column,}       the operator is called {\em  spray} or {\em  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}
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:
\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 \\