Model fitting by least squares

KRYLOV SUBSPACE ITERATIVE METHODS

The solution time for simultaneous linear equations grows cubically with the number of unknowns. There are three regimes for solution; which one is applicable depends on the number of unknowns $m$. For $m$ three or less, we use analytical methods. We also sometimes use analytical methods on matrices of size $4\times 4$ if the matrix contains many zeros. A 1988 my desktop workstation solved a $100 \times 100$ system in a minute. Ten years later it would do a $600\times 600$ system in about a minute. A nearby more powerful computer would do $1000\times 1000$ in a minute. Since the computing effort increases with the third power of the size, and since $4^3=64\approx 60$, an hour's work solves a four times larger matrix, namely $4000\times 4000$ on the more powerful machine. For significantly larger values of $m$, exact numerical methods must be abandoned and iterative methods must be used.

(The compute time for a rectangular matrix is slightly more pessimistic. It is the product of the number of data points $n$ times the number of model points squared $m^2$. This happens to be the cost of computing the matrix $\bold F\T\bold F$ from $\bold F$. Since the number of data points generally exceeds the number of model points $n>m$ by a substantial factor (to allow averaging of noises), it leaves us with significantly fewer than 4000 points in model space.)

A square image packed into a 4096 point vector is a $64\times 64$ array. The computer power for linear algebra to give us solutions that fit in a $k\times k$ image is thus proportional to $k^6$, which means that even though computer power grows rapidly, imaging resolution using ``exact numerical methods'' hardly grows at all from our $64\times 64$ current practical limit.

The retina in our eyes captures an image of size about $1000\times 1000$ which is a lot bigger than $64\times 64$. Life offers us many occasions where final images exceed the 4000 points of a $64\times 64$ array. To make linear algebra (and inverse theory) relevant to such applications, we investigate special techniques. A numerical technique known as the ``conjugate-direction method'' works well for all values of $m$ and is our subject here. As with most simultaneous equation solvers, an exact answer (assuming exact arithmetic) is attained in a finite number of steps. And if $n$ and $m$ are too large to allow enough iterations, the iterative methods can be interrupted at any stage, the partial result often proving useful. Whether or not a partial result actually is useful is the subject of much research; naturally, the results vary from one application to the next.

Model fitting by least squares