[petsc-users] KSPLSQR convergence criterion

[Chetan Jhurani]

Replying to my previous mail, this time with self-contained
Petsc code if it helps.  As before this is to use KSPLSQR
to solve

min_x ||Ax - b||, with A = [1;2], b = [3;4],

Results:

KSPGetIterationNumber 10000
KSPGetResidualNorm 0.894427
KSPConvergedReason -3

-------------------------------------------------------------

#include <petscsys.h>
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>
#include <iostream>

int main(int argc, char* argv[])
{
	PetscErrorCode ierr;

	ierr = PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL); CHKERRQ(ierr);

	PetscScalar A_array[2] = {1, 2};
	PetscScalar rhs_array[2] = {3, 4};
	PetscInt pos[2] = {0, 1};
	Mat A;
	Vec soln, rhs;
	KSP ksp;

	ierr = MatCreateSeqDense(PETSC_COMM_SELF, 2, 1, A_array, &A); CHKERRQ(ierr);
	ierr = VecCreateSeq(PETSC_COMM_SELF, 2, &rhs); CHKERRQ(ierr);
	ierr = VecSetValues(rhs, 2, pos, rhs_array, INSERT_VALUES); CHKERRQ(ierr);
	ierr = VecAssemblyBegin(rhs); CHKERRQ(ierr);
	ierr = VecAssemblyEnd(rhs); CHKERRQ(ierr);
	ierr = VecCreateSeq(PETSC_COMM_SELF, 1, &soln); CHKERRQ(ierr);
	ierr = KSPCreate(PETSC_COMM_SELF, &ksp); CHKERRQ(ierr);
	ierr = KSPSetType(ksp, KSPLSQR); CHKERRQ(ierr);
	ierr = KSPSetOperators(ksp, A, A, SAME_PRECONDITIONER); CHKERRQ(ierr);
	ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
	ierr = KSPSolve(ksp, rhs, soln); CHKERRQ(ierr);

	PetscInt num_iters;
	ierr = KSPGetIterationNumber(ksp, &num_iters); CHKERRQ(ierr);

	PetscReal rnorm;
	ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr);

	KSPConvergedReason reason;
	ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr);

	std::cout << "KSPGetIterationNumber " << num_iters << '\n';
	std::cout << "KSPGetResidualNorm " << rnorm << '\n';
	std::cout << "KSPConvergedReason " << reason << '\n';

	ierr = KSPDestroy(ksp); CHKERRQ(ierr);
	ierr = VecDestroy(soln); CHKERRQ(ierr);
	ierr = VecDestroy(rhs); CHKERRQ(ierr);
	ierr = MatDestroy(A); CHKERRQ(ierr);

	ierr = PetscFinalize(); CHKERRQ(ierr);
}

-------------------------------------------------------------

Chetan

> -----Original Message-----
> From: Chetan Jhurani [mailto:chetan.jhurani at gmail.com]
> Sent: Friday, August 06, 2010 10:40 AM
> To: 'PETSc users list'
> Subject: KSPLSQR convergence criterion
> 
> 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
> 
> ----------------------------------------------------
>