[petsc-users] KSPLSQR convergence criterion

[Barry Smith]

   Ok, here is my understanding of the situation:  LSQR actually solves the normal equations A'*Ax = A'x. When Ax = b has a solution then Matlab and PETSc use || b - Ax|| < tol to determine convergence. Matlab also checks || A'*b - A'*Ax || < tol for convergence.

   In petsc-dev http://www.mcs.anl.gov/petsc/petsc-as/developers/index.html I have added an additional test based on the residual norm of the normal equations || A'*b - A'*Ax || < tol    I also added -ksp_monitor_lsqr that prints the residual norm and the residual norm of the normal equations. For example

[barry-smiths-macbook-pro:~/Downloads] barrysmith% ./ex1 -ksp_max_it 10 -ksp_converged_reason -ksp_monitor_lsqr
  0 KSP Residual norm 5.000000000000e+00
  1 KSP Residual norm 8.944271909999e-01 Residual norm normal equations 1.807443084090e-15
Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1

Please let me know if you have any troubles.

  Barry

Note that I use a different scaling of tol for convergence test than Matlab (they try to approximate the norm of the preconditioned matrix and include that inside the bound, in PETSc we never do this).



On Aug 6, 2010, at 11:40 AM, Chetan Jhurani wrote:

> Hello!
> 
> Could someone help me figure out the correct way to
> to use KSPLSQR?  Maybe I'm missing something simple.
> 
> I'm trying to use it to solve an over-determined set
> of equations.  It seems like the KSP iterations keep on
> going because the residual is not small enough.  But
> since the system is over-determined, in general, the
> exact residual Ax - b will not be zero.
> 
> For example, to solve the system
> 
> min_x ||Ax - b||, with A = [1;2], b = [3;4],
> 
> x must be inv(A'A)A'b = 2.2.  Relative residual is 0.17888.
> See the Matlab code below too, which gives the correct
> value in 1 iteration.
> 
> KSPLSQR, however, stagnates after 1 iteration.  At that
> point, the relative residual and true residual norm
> match the correct values exactly.
> 
> Is there an assumption in KSPLSQR that the system is always
> invertible or under-determined?  This is petsc-3.1p2 if
> it matters.
> 
> Thanks,
> 
> Chetan
> 
> ----------------------------------------------------
> KSPLSQR:
> 
>  0 KSP preconditioned resid norm 5.000000000000e+000 true resid norm 5.000000000000e+000 ||Ae||/||Ax|| 1.000000000000e+000
>  1 KSP preconditioned resid norm 8.944271909999e-001 true resid norm 8.944271909999e-001 ||Ae||/||Ax|| 1.788854382000e-001
>  2 KSP preconditioned resid norm 8.944271909999e-001 true resid norm 8.944271909999e-001 ||Ae||/||Ax|| 1.788854382000e-001
>  3 KSP preconditioned resid norm 8.944271909999e-001 true resid norm 8.944271909999e-001 ||Ae||/||Ax|| 1.788854382000e-001
>  4 KSP preconditioned resid norm 8.944271909999e-001 true resid norm 8.944271909999e-001 ||Ae||/||Ax|| 1.788854382000e-001
>  5 KSP preconditioned resid norm 8.944271909999e-001 true resid norm 8.944271909999e-001 ||Ae||/||Ax|| 1.788854382000e-001
> 
> Matlab:
> 
>>> a = [1;2];
>>> b = [3;4];
>>> [x,FLAG,RELRES,ITER] = lsqr(a,b, 1e-10, 10)
> x =
>   2.200000000000000
> FLAG =
>     1
> RELRES =
>   0.178885438199983
> ITER =
>     1
>>> norm(a*x-b)
> ans =
>   0.894427190999916
> 
> ----------------------------------------------------
> 
>