INTRODUCTION

Geophysical mapping and imaging are applications where we seek an approximate pseudo inverse of a matrix of very high order. Say, $\bold d = \bold F\bold m$ constructs theoretical data $\bold d$from model parameters $\bold m$ using a linear operator $\bold F$.Experience shows that the transpose of the simulation operator provides a useful image. Thus the suggested image is $\hat {\bold m} = \bold F' \bold d$ where $\bold d$ is the observed data, and $\bold F'$ is the transpose (Hilbert adjoint) of $\bold F$.Mathematically, this means that when the dimensionality is very high, we often approximate an inverse by a transpose $\bold F^{-1} \approx \bold F'$.

Experience shows that a ``better'' image is often created by a unitary operator. Practitioners often discover improvements in the form of scaling diagonal matrices, say $\bold F^{-1} \approx \bold D_1 \bold F' \bold D_2$ where $\bold D_1$ and $\bold D_2$ are the scaling diagonal matrices. We seek an operator $\bold G = \bold D_1 \bold F' \bold D_2$that is as unitary as possible, i.e., $\bold G' \bold G \approx (\bold G' \bold G)^2$.We believe that $\bold G$is also in some sense ideally preconditioned for linear solvers. A basic, recurring problem is a lack of straightforward theory for finding $\bold D_1$ and $\bold D_2$.