PROBLEM SETTING

A reflectivity model $\bold m$ in geophysics is a function defined in a physical volume sampled densely enough to fill a computer's random access memory. The volume typically has somewhat fewer than about (10³)³ = 10⁹ elements (voxels). The data $\bold d$ upon which the model is built is a hundred to a thousand times larger, say 10¹¹ elements and generally requires many magnetic tapes. Thus the matrix $\bold F$ relating the data to the model would have 10²⁰ elements (which is a futile way to think about the problem). What we can readily do is manufacture synthetic data from any model by running a program that amounts to a matrix multiply

$$\bold d_{\rm theoretical} \quad=\quad\bold F \bold m$$ (1)

By paying attention to some details it is straight-forward (in all the problems that I have encountered) to write another program that applies the transpose $\bold F'$ operator to the data $\bold d$ giving us a first guess at the model. This first guess is called an image $\hat {\bold m}$ of the model $\bold m$.

$$\hat {\bold m} \quad=\quad \bold F' \bold d_{\rm observed}$$ (2)

The method of least squares offers the standard solution

$$\hat {\bold m} \quad=\quad(\bold F'\bold F)^{-1} \bold F' \bold d_{\rm observed}$$ (3)

which is wholly impractical because of the size of $\bold F'\bold F$.This suggests we seek a diagonal matrix, say $\bold D_m$ as an approximation for $(\bold F'\bold F)^{-1}$.

$$\hat {\bold m} \quad=\quad\bold D_m \bold F' \bold d_{\rm observed}$$ (4)

At my request, Bill Symes offered a suggestion for $\bold D_m$.His suggestion uses the ${\bf diag}$ operator of matrix algebra which simply takes a vector and distributes its components along the diagonal of a matrix:

$$\bold D_m \quad=\quad { {\bf diag}\ \bold F' \bold d \over ... ...\bold F \bold F' \bold d } \quad=\quad {\bf diag}\ {a \over b}$$ (5)

This suggestion seems reasonable because the denominator has two more powers of $\bold F$ than does the numerator and the suggestion seems wise because it uses the known data $\bold d$.Unfortunately, the denominator b can easily vanish. To prevent such trouble we rearrange and smooth

$$\bold D_m \quad=\quad{\bf diag}\ { <a b\gt \over <b^2\gt + \epsilon^2 }$$ (6)

where $\epsilon$ is a positive number for experimentation and the numerator and denominator can be locally averaged as indicated by the ``$<\ \gt$''. (Remember that the vectors $\bold m$ and $\bold d$ are sampled representations of continuous functions of space and time.) Our experience with Symes' suggestion was disappointing, probably for the reason suggested next.

Besides the standard pseudo inverse $\hat {\bold m} = (\bold F'\bold F)^{-1} \bold F' \bold d$there is another

$$\hat {\bold m} \quad=\quad\bold F' (\bold F\bold F')^{-1} \bold d_{\rm observed}$$ (7)

that arises when a fitting problem is underdetermined. This pseudoinverse solves the problem of minimizing $\vert\vert\hat {\bold m}\vert\vert$subject to the constraint equations $\bold d = \bold F\bold m$.This pseudoinverse suggests a right side diagonal approximation

$$\hat {\bold m} \quad=\quad\bold F' \bold D_d \bold d_{\rm observed}$$ (8)

I don't expect this right-side diagonal to work any better than did the earlier left-side diagonal. Experience with special cases suggests we need both the data-space diagonal and the model-space diagonal, namely

$$\hat {\bold m} \quad=\quad\bold D_m \bold F' \bold D_d \bold d_{\rm observed}$$ (9)

Our challange for the mathematicians is to give us some guesses for the diagonal matrices $\bold D_m$ and $\bold D_d$.This problem is one of a very general nature, far more general than the several applications in which I have found it. It is like the question, ``What is an all-purpose preconditioner?''