paper.tex

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

\def\cakedir{.}
\def\figdir{./Fig}
\lefthead{Fomel}
\righthead{Conjugate directions}
\footer{SEP--92}

\lstset{language=c,numbers=left,numberstyle=\tiny,showstringspaces=false}

\title{Least-square inversion with inexact adjoints. \\
Method of conjugate directions: A tutorial}
%\keywords{inversion, algorithm, modeling, linear }

\email{[e-mail]}
\author{Sergey Fomel}

\maketitle

\begin{abstract}
  This tutorial describes the classic method of conjugate directions:
  the generalization of the conjugate-gradient method in iterative
  least-square inversion. I derive the algebraic equations of the
  conjugate-direction method from general optimization principles. The
  derivation explains the ``magic'' properties of conjugate gradients.
  It also justifies the use of conjugate directions in cases when
  these properties are distorted either by computational errors or by
  inexact adjoint operators. The extra cost comes from storing a
  larger number of previous search directions in the computer memory.
  A simple program and two examples illustrate the method.
\end{abstract}


\section{Introduction}
%%%%

This paper describes the method of conjugate directions for solving
linear operator equations in Hilbert space. This method is usually
described in the numerous textbooks on unconstrained optimization as
an introduction to the much more popular method of conjugate
gradients. See, for example, {\em Practical optimization} by
\cite{gill} and its bibliography. The famous conjugate-gradient solver
possesses specific properties, well-known from the original works of
\cite{hestenes} and \cite{fletcher}. For linear operators and exact
computations, it guarantees finding the solution after, at most, $n$
iterative steps, where $n$ is the number of dimensions in the solution
space. The method of conjugate gradients doesn't require explicit
computation of the objective function and explicit inversion of the
Hessian matrix.  This makes it particularly attractive for large-scale
inverse problems, such as those of seismic data processing and
interpretation.  However, it does require explicit computation of the
adjoint operator. \cite{Claerbout.blackwell.92,iee} shows dozens of
successful examples of the conjugate gradient application with
numerically precise adjoint operators.
\par
The motivation for this tutorial is to explore the possibility of
using different types of preconditioning operators in the place of
adjoints in iterative least-square inversion. For some linear or
linearized operators, implementing the exact adjoint may pose a
difficult problem. For others, one may prefer different
preconditioners because of their smoothness
\cite[]{Claerbout.sep.89.201,Crawley.sep.89.207}, simplicity
\cite[]{kleinman}, or asymptotic properties \cite[]{herman}. In those
cases, we could apply the natural generalization of the conjugate
gradient method, which is the method of conjugate directions. The cost
difference between those two methods is in the volume of memory
storage. In the days when the conjugate gradient method was invented,
this difference looked too large to even consider a practical
application of conjugate directions. With the evident increase of
computer power over the last 30 years, we can afford to do it now.
\par
I derive the main equations used in the conjugate-direction method
from very general optimization criteria, with minimum restrictions
implied. The textbook algebra is illustrated with a simple program and
two simple examples.

\section{IN SEARCH OF THE MINIMUM} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We are looking for the solution of the linear operator equation
\begin{equation}
{\bf d = A\,m}\;,
\label{eqn:equation}
\end{equation} 
where ${\bf m}$ is the unknown model in the linear model space, ${\bf
d}$ stands for the given data, and ${\bf A}$ is the forward modeling
operator. The data vector ${\bf d}$ belongs to a Hilbert space with
a defined norm and dot product. The solution is constructed by iterative
steps in the model space, starting from an initial guess ${\bf
m}_0$. Thus, at the $n$-th iteration, the current model ${\bf m}_n$ is
found by the recursive relation
\begin{equation}
{\bf m}_n = {\bf m}_{n-1} + \alpha_n {\bf s}_n\;,
\label{eqn:mn}
\end{equation}
where ${\bf s}_n$ denotes the step direction, and $\alpha_n$ stands
for the scaling coefficient. The residual at the $n$-th iteration is
defined by
\begin{equation}
{\bf r}_n = {\bf d - A\,m}_{n}\;.
\label{eqn:residual}
\end{equation}
Substituting (\ref{eqn:mn}) into (\ref{eqn:residual}) leads to the equation
\begin{equation}
{\bf r}_n = {\bf r}_{n-1} - \alpha_n {\bf A\,s}_n\;.
\label{eqn:rn}
\end{equation}
For a given step ${\bf s}_n$, we can choose $\alpha_n$ to minimize the
squared norm of the residual
\begin{equation}
\|{\bf r}_n\|^2 = \|{\bf r}_{n-1}\|^2 - 
2\,\alpha_n \left({\bf r}_{n-1},\,{\bf A\,s}_n\right) +
\alpha_n^2\,\|{\bf A\,s}_n\|^2\;.
\label{eqn:rnorm}
\end{equation}
The parentheses denote the dot product, and 
$\|{\bf x}\|=\sqrt{({\bf x,\,x})}$ denotes the norm of $x$ in the
corresponding Hilbert space. The optimal value of $\alpha_n$ is easily
found from equation (\ref{eqn:rnorm}) to be
\begin{equation}
\alpha_n = {{\left({\bf r}_{n-1},\,{\bf A\,s}_n\right)} \over
{\|{\bf A\,s}_n\|^2}}\;.
\label{eqn:alpha}
\end{equation}
Two important conclusions immediately follow from this fact. First,
substituting the value of $\alpha_n$ from formula (\ref{eqn:alpha}) into
equation (\ref{eqn:rn}) and multiplying both sides of this equation by ${\bf
r}_n$, we can conclude that
\begin{equation}
\left({\bf r}_n,\,{\bf A\,s}_n\right) = 0\;,
\label{eqn:rasn}
\end{equation}
which means that the new residual is orthogonal to the corresponding
step in the residual space. This situation is schematically shown in
Figure \ref{fig:dirres}. Second, substituting formula (\ref{eqn:alpha}) into 
(\ref{eqn:rnorm}), we can conclude that the new residual decreases according
to
\begin{equation}
\|{\bf r}_n\|^2 = \|{\bf r}_{n-1}\|^2 - 
{{\left({\bf r}_{n-1},\,{\bf A\,s}_n\right)^2} \over
{\|{\bf A\,s}_n\|^2}}\;,
\label{eqn:pythagor}
\end{equation}
(``Pythagoras's theorem'' ), unless ${\bf r}_{n-1}$ and ${\bf A\,s}_n$
are orthogonal. These two conclusions are the basic features of
optimization by the method of steepest descent. They will help us
define an improved search direction at each iteration.

\inputdir{XFig}

\sideplot{dirres}{height=2.5in}{Geometry of the residual in the
data space (a scheme).}

\section{IN SEARCH OF THE DIRECTION} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Let's sup\-pose we have a ge\-ne\-ra\-tor that pro\-vides
parti\-cu\-lar search directions at each step. The new direction can
be the gradient of the objective function (as in the method of
steepest descent), some other operator applied on the residual from
the previous step, or, generally speaking, any arbitrary vector in the
model space. Let us denote the automatically generated direction by
${\bf c}_n$. According to formula (\ref{eqn:pythagor}), the residual
decreases as a result of choosing this direction by
\begin{equation}
\|{\bf r}_{n-1}\|^2 - \|{\bf r}_{n}\|^2 = 
{{\left({\bf r}_{n-1},\,{\bf A\,c}_n\right)^2} \over {\|{\bf
A\,c}_n\|^2}}\;.
\label{eqn:deltar}
\end{equation}
How can we improve on this result? 

\subsection{First step of the improvement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Assuming $n>1$, we can add some
amount of the previous step ${\bf s}_{n-1}$ to the chosen direction
${\bf c}_n$ to produce a new search direction ${\bf s}_n^{(n-1)}$, as
follows:
\begin{equation}
{\bf s}_n^{(n-1)} =  {\bf c}_n + \beta_n^{(n-1)}\,{\bf s}_{n-1}\;,
\label{eqn:cn}
\end{equation}
where $\beta_n^{(n-1)}$ is an adjustable scalar coefficient. According to
to the fundamental orthogonality principle (\ref{eqn:rasn}), 
\begin{equation}
\left({\bf
r}_{n-1},\,{\bf A\,s}_{n-1}\right) = 0\;.  
\label{eqn:rasn1}
\end{equation}
As follows from equation (\ref{eqn:rasn1}), the numerator on the right-hand
side of equation (\ref{eqn:deltar}) is not affected by the new choice of the
search direction:
\begin{equation}
\left({\bf r}_{n-1},\,{\bf A\,s}_n^{(n-1)}\right)^2 = \left[
\left({\bf r}_{n-1},\,{\bf A\,c}_n\right) + \beta_n^{(n-1)}\,
\left({\bf r}_{n-1},\,{\bf A\,s}_{n-1}\right)\right]^2 =
\left({\bf r}_{n-1},\,{\bf A\,c}_n\right)^2\;.
\label{eqn:numerator}
\end{equation}
However, we can use transformation (\ref{eqn:cn}) to decrease the
denominator in (\ref{eqn:deltar}), thus further decreasing the residual
${\bf r}_n$. We achieve the minimization of the denominator
\begin{equation}
\|{\bf A\,s}_n^{(n-1)}\|^2 =  \|{\bf A\,c}_n\|^2 + 
2\,\beta_n^{(n-1)}\,\left({\bf A\,c}_n,\,{\bf A\,s}_{n-1}\right) +
\left(\beta_n^{(n-1)}\right)^2\,\|{\bf A\,s}_{n-1}\|^2
\label{eqn:denominator}
\end{equation} 
by choosing the coefficient $\beta_n^{(n-1)}$ to be
\begin{equation}
\beta_n^{(n-1)} = - {{\left({\bf A\,c}_n,\,{\bf A\,s}_{n-1}\right)} \over
{\|{\bf A\,s}_{n-1}\|^2}}\;.
\label{eqn:beta}
\end{equation}
Note the analogy between (\ref{eqn:beta}) and (\ref{eqn:alpha}). Analogously to
(\ref{eqn:rasn}), equation (\ref{eqn:beta}) is equivalent to the orthogonality condition
\begin{equation}
\left({\bf A\,s}_n^{(n-1)},\,{\bf A\,s}_{n-1}\right) = 0\;.
\label{eqn:acas}
\end{equation}
Analogously to (\ref{eqn:pythagor}), applying formula (\ref{eqn:beta}) is also equivalent to defining the
minimized denominator as
\begin{equation}
\|{\bf A\,c}_n^{(n-1)}\|^2 =  \|{\bf A\,c}_n\|^2 -
{{\left({\bf A\,c}_n,\,{\bf A\,s}_{n-1}\right)^2} \over
{\|{\bf A\,s}_{n-1}\|^2}}\;.
\label{eqn:pithagor2}
\end{equation}

\subsection{Second step of the improvement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Now let us assume $n > 2$ and add some amount of the step from
the $(n-2)$-th iteration to the search direction, determining the new
direction ${\bf s}_n^{(n-2)}$, as follows:
\begin{equation}
{\bf s}_n^{(n-2)} =  {\bf s}_n^{(n-1)} + \beta_n^{(n-2)}\,{\bf s}_{n-2}\;.
\label{eqn:cn2}
\end{equation}
We can deduce that after the second change, the value of numerator in
equation (\ref{eqn:deltar}) is still the same:
\begin{equation}
\left({\bf r}_{n-1},\,{\bf A\,s}_n^{(n-2)}\right)^2 = \left[
\left({\bf r}_{n-1},\,{\bf A\,c}_n\right) + \beta_n^{(n-2)}\,
\left({\bf r}_{n-1},\,{\bf A\,s}_{n-2}\right)\right]^2 =
\left({\bf r}_{n-1},\,{\bf A\,c}_n\right)^2\;.
\label{eqn:numerator2}
\end{equation}
This remarkable fact occurs as the result of transforming the dot product $\left({\bf
r}_{n-1},\,{\bf A\,s}_{n-2}\right)$ with the help of equation
(\ref{eqn:rn}):
\begin{equation}
\left({\bf r}_{n-1},\,{\bf A\,s}_{n-2}\right) =
\left({\bf r}_{n-2},\,{\bf A\,s}_{n-2}\right) -
\alpha_{n-1}\,\left({\bf A\,s}_{n-1},\,{\bf A\,s}_{n-2}\right) = 0\;.
\label{eqn:dotprod}
\end{equation}
The first term in (\ref{eqn:dotprod}) is equal to zero according to formula
(\ref{eqn:rasn}); the second term is equal to zero according to formula
(\ref{eqn:acas}). Thus we have proved the new orthogonality equation
\begin{equation}
\left({\bf r}_{n-1},\,{\bf A\,s}_{n-2}\right) = 0\;,
\label{eqn:rasn2}
\end{equation}
which in turn leads to the numerator invariance (\ref{eqn:numerator2}). The
value of the coefficient $\beta_n^{(n-2)}$ in (\ref{eqn:cn2}) is defined
analogously to (\ref{eqn:beta}) as
\begin{equation}
\beta_n^{(n-2)} = - 
{{\left({\bf A\,s}_n^{(n-1)},\,{\bf A\,s}_{n-2}\right)} \over
{\|{\bf A\,s}_{n-2}\|^2}} = 
- {{\left({\bf A\,c}_n,\,{\bf A\,s}_{n-2}\right)} \over
{\|{\bf A\,s}_{n-2}\|^2}}\;,
\label{eqn:beta2}
\end{equation}
where we have again used equation (\ref{eqn:acas}). If ${\bf A\,s}_{n-2}$ is
not orthogonal to ${\bf A\,c}_n$, the second step of the improvement leads
to a further decrease of the denominator in (\ref{eqn:pythagor}) and,
consequently, to a further decrease of the residual.

\subsection{Induction}
%%%%%%%%%%%%%%%%%%
Continuing by induction the process of adding a linear combination of
the previous steps to the arbitrarily chosen direction ${\bf c}_n$
(known in mathematics as the {\em Gram-Schmidt orthogonalization
process}), we finally arrive at the complete definition of the new
step ${\bf s}_n$, as follows:
\begin{equation}
{\bf s}_n = {\bf s}_n^{(1)} =  
{\bf c}_{n} + \sum_{j=1}^{j=n-1}\,\beta_n^{(j)}\,{\bf s}_{j}\;.
\label{eqn:step}
\end{equation}
Here the coefficients $\beta_n^{(j)}$ are defined by equations
\begin{equation}
\beta_n^{(j)} = 
- {{\left({\bf A\,c}_n,\,{\bf A\,s}_{j}\right)} \over
{\|{\bf A\,s}_{j}\|^2}}\;,
\label{eqn:betaj}
\end{equation} 
which correspond to the orthogonality principles
\begin{equation}
\left({\bf A\,s}_n,\,{\bf A\,s}_{j}\right) = 0\;,\;\;1 \leq j \leq n-1
\label{eqn:asj}
\end{equation}
and
\begin{equation}
\left({\bf r}_{n},\,{\bf A\,s}_{j}\right) = 0\;,\;1 \leq j \leq n\;.
\label{eqn:rasnj}
\end{equation}
It is these orthogonality properties that allowed us to optimize the
search parameters one at a time instead of solving the $n$-dimensional
system of optimization equations for $\alpha_n$ and $\beta_n^{(j)}$.

\section{ALGORITHM}
The results of the preceding sections define the method of conjugate
directions to consist of the following algorithmic steps: