Modified LSQR

I can resolve this difficulty at the end of the calculation by computing the data resolution using the formula

_data = _r = _r(_r)^-1_r^T,   instead of (lsqrdatares). Formula (dataresLSQR) implies that I actually multiply the matrices and $\matX_r$ to find the true resolution, instead of taking the inappropriate shortcut indicated in (lsqrdatares). It is not obvious whether (dataresLSQR) is symmetric or not, although I expect that it is. Another approach is to use the Moore-Penrose uniqueness condition for the generalized inverse to strip off the parts of $\matX_r$ that lie in the null space of , i.e.,

^_r_r,   where again the matrix multiplies appearing in the formula must actually be executed in order to achieve the required stripping effect. I call (LSQRMPinverse) the Modified LSQR inverse. Once $\matM^\dagger$ has been constructed this way, (solutionvecs) is replaced by

= ^= _r_r   and I can use the standard formula $\calR_{data} = \matM\matM^\dagger$ to compute the fully symmetric data resolution matrix, if this is required.

This iterative procedure may be terminated early (i.e., for k < r) by setting b_2k+1 = 0. Then, an approximate LSQR inverse is given by

_k = _k(_k)^-1_k^T.   It is easy to check that $\matX_k\matM\matX_k \ne \matX_k$ and that $\matM\matX_k\matM \ne \matM$ if k < r. The effective model resolution matrix is given easily for the LSQR method by

_model = _k_k^T = _i=1^k ^(2i) (^(2i))^T.   The same problems that were discussed above for $\calR_{data}$ apply to $\calE_{data}$ as well, so the formula for effective data resolution in this situation is

_data = _k = _k(_k)^-1_k^T.   Effective resolution $\calE_{model}$ is clearly symmetric for any choice of k. I expect $\calE_{data}$ is also symmetric, but it is again not immediately obvious from the formula. Similarly, the Modified LSQR for partial iteration then uses $\matX_k\matM\matX_k$ as the inverse, rather than $\matX_k$ itself.

Not surprisingly, when and are consistent, the LSQR algorithm has the same qualitative characteristics as the Lanczos method from which it derives. The main advantage of this method is that it produces the generalized inverse of directly, rather than producing $\left(\matM^T\matM\right)^\dagger$.Thus, the dynamic range problems due to poor conditioning of are aggravated when the Lanczos method is applied directly to the normal equations since the eigenvalues are squared in $\matM^T\matM$; in contrast, no such aggravation occurs in the LSQR algorithm. Both algorithms are qualitatively similar to truncated SVD, but they are easier to implement and more economical than many SVD routines. Of course, both methods can be turned into SVD routines simply by including an additional step to find the eigenvalues and eigenvectors of the tridiagonal matrix $\matT_r = \matB_r^T\matB_r$. However, such a step is not necessary either to solve the inversion problem or to compute the resolution matrices.

One major difference between LSQR and the method of Lanczos is that storage is greater for LSQR by a factor of (m+n)/n, since both $\matG_k$and $\matH_k$ must be stored, especially when the resolution matrices must be computed. This difference can be significant if the data dimension m is much greater than the model dimension n as it is for overdetermined problems.