THE WINNERS

For applications with small numbers of unknowns, say less than dozens or hundreds, the damped inverse operator is generally the most appropriate, whereas for problems when the unknown model is a field of more than thousands of unknowns, say reflectivity or velocity as a function of two coordinates, then the conjugate operator and its various extensions are, in practice, the answer.

Between the conjugate operator and the inverse operator lies an operator (often overlooked) that is sometimes a guide to the most effective compromise. This operator is a certain unitary operator. The solution and the unitary operator are given by  

$$\bold m_{\rm unitary} \quad =\quad[( \bold B'\bold B)^{-1/2}\ \bold B'] \bold d$$ (5)

If $\bold B'\bold B$ tends to a singular matrix, then $(\bold B'\bold B)^{-1/2}\, \bold B'$ is like 0/0 whereas $(\bold B'\bold B)^{-1}\, \bold B'$ is like 0/0² so equation (5) is much safer than (2). Equation (5) is exceedingly ambiguous because of the huge multiplicity of ways of representing the square root of an operator. For an operator of N eigenvalues there are 2^N square roots corresponding to choice of sign of each. Generally, there exists an application-specific reason for selecting one of the square roots. We might seek a symmetric phaseless operator, for example, or we might seek a causal one. A clean and illustrative example of such an operator is unitary normal moveout with linear interpolation Claerbout (1986) or 1992b, page 118. This example begins from the observation that linear interpolation causes smoothing which suppresses high frequencies and noticing that it would be nice if the transformation conserved energy. The linear interpolation is represented mathematically by an NMO transformation matrix $\bold B$containing adjacent elements that sum to unity (linear interpolation) leading to $\bold B'\bold B$ being a tridiagonal matrix. To build unitary NMO, the tridiagonal matrix $\bold B'\bold B$ is factored into a bidiagonal matrix and its conjugate, say $\bold A'\bold A$ by an obvious method that is known as the Cholesky decomposition. Then (5) reduces to  

$$\bold A \bold m_{\rm factored} \quad =\quad\bold B' \bold d$$ (6)

which is easily solved for $\bold m$. by recursion. Formally, the unitary operator is $\bold A^{-1}\bold B'$.

More typically, in practice the square-root operator is not easily found, nor is it needed precisely. This leads to  

$$\bold m_{\rm precond} \quad =\quad \ll\bold B'\bold B\gg^{-1/2} \ \bold B' \bold d$$ (7)

which has the same form as (5) except that we build various simplifications into the definition of $\ll\bold B'\bold B\gg$.A notable example of this concept arises in time-series analysis where the symbols $\ll\ \gg$ denote ensemble smoothing and the factorization is not done by the Cholesky algorithm but by that of Kolmogoroff described in Claerbout (1992b) on page 230.

One of the most important linear operators in seismic exploration is velocity analysis (summing along hyperbolic trajectories of various velocities). After having defined such a transform, you can push data into velocity space and pull it back again with the operator $\bold B\bold B'$.Unless by chance or design you have defined your operator to be already unitary, you will find, when comparing the original data $\bold d$ to $\bold B\bold B'\bold d$that the processing weakens events at various times, locations, or frequencies. The remedy for this is to redefine the operator with various scale factors, such as $\sqrt{-i\omega}\sqrt{x/v}\sqrt{x/t}\sqrt{1/t}$that you can uncover either experimentally, theoretically, or by looking at Claerbout 1992a All this activity amounts to looking for a unitary representation such as equation (7). In my experience in geophysics, such efforts are not as systematic as they should be and could be. We should make a greater effort to digest the literature on preconditioning. Further, the Fortran language itself inhibits segregation of numerical-analysis concepts from application-specific concepts. Thus successful applications of advanced numerical analysis concepts are many fewer than we would like to see. This motivates a C++ study in this report Nichols et al. (1992).

Finally, I wish to cite another example of operator factorization for the inverse needed by (7). That example is the two-dimensional prediction-error operator in Abma and Claerbout, Abma and Claerbout (1992). The confusion caused by that paper at the 1992 SEP sponsor meeting is one reason I prepared this short note. The other reason, is to bring up again our lack of systematic use of preconditioning when we attack huge inverse problems. To begin with, I believe that the practical search for unitary operators is closely related to the practical search for preconditioners, though I have not formally confirmed it.