FAQ1: This question and answer deals with providing a
      starting point for LSQR (and other iterative solvers).

--------------------------------------------------------------------------
From: "Forrester H Cole (fcole@Princeton.EDU)" <fcole@Princeton.EDU>
To: saunders@stanford.edu
Subject: lsqr bnorm
Date-Sent: Tuesday, August 16, 2005 07:38:07 AM

Dear Michael,

  I'm still working with LSQR. I noticed that it seemed to be
doing too much work (~50 iterations or so) when I gave it a
system that was very nearly solved already - i.e., A*x0 - b ~= 0.
I noticed that in this case, the reported value of bnorm was way
too low, assuming bnorm represents the euclidean norm of the
right-hand side vector b. In the code, there is this (around line
475):

dvec_scale( (-1.0), input->rhs_vec );
  func->mat_vec_prod( 0, input->sol_vec, input->rhs_vec, prod );
  dvec_scale( (-1.0), input->rhs_vec );
  
    /* compute Euclidean length of u and store as BETA */
  beta = dvec_norm2( input->rhs_vec );
  
...
...
...

 bnorm = beta;

It seems like bnorm is being set not to the norm of the original
rhs vector, but to the norm of the initial residual vector. I
moved the setting of bnorm above the code that changes
input->rhs_vec. This change gives me the output I expect, namely
that if the initial residual is very small LSQR terminates after
only 1 iteration.

Is this a correct change? I haven't noticed any problems with it
so far, but I may be missing something.

Thanks,
Forrester Cole
---------------------------------------------------------------------------


Michael A. Saunders wrote:

Dear Forrester,

Thanks for explaining your question clearly.
It goes back to the reason why my own F77 and Matlab versions
of LSQR do not allow you to provide a starting point x0.
If you have such an x0, I want users to be conscious of
computing the residual
   r0 = b - A x0,
solving the system
   min ||A dx - r0||,
and then adding
   x = x0 + dx
*after* the call to LSQR.  The only way you can obtain benefit
from x0 is to recognize that the solve for dx does not need
to be so accurate.

So -- the answer is that the C code was setting bnorm = ||r0||
correctly, but you need to use larger values for atol and btol
to cause earlier termination in the solve for dx.

It's hard to automate this process, but you could base it on
the value of ||A'r0|| compared to ||A'b||.  If x0 is good,
||A'r0|| should be small, so atol and btol could be scaled up
accordingly.

Note that if x0 is *really* good, dx will be very small.
Hence we don't want LSQR to be adding corrections to x0
during the iterations, as many digits will be lost.

Good luck,
Michael
--------------------------------------------------------------

From: Forrester Cole <fcole@Princeton.EDU>
To: saunders@stanford.edu
Subject: Re: lsqr bnorm
Date-Sent: Tuesday, August 16, 2005 07:58:58 PM

Dear Michael,
  Thanks, that makes sense. Better still, it seems to work!

Forrester
--------------------------------------------------------------