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 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 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 has been constructed this way, (solutionvecs) is replaced by
= ^= _r_r and I can use the standard formula 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 b2k+1 = 0. Then, an approximate LSQR inverse is given by
_k = _k(_k)^-1_k^T. It is easy to check that and that 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 apply to as well, so the formula for effective data resolution in this situation is
_data = _k = _k(_k)^-1_k^T. Effective resolution is clearly symmetric for any choice of k. I expect is also symmetric, but it is again not immediately obvious from the formula. Similarly, the Modified LSQR for partial iteration then uses as the inverse, rather than 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 .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 ; 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 . 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 and 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.