[petsc-users] Algorithms to remove null spaces in a singular system

[Barry Smith]

> On Oct 12, 2016, at 10:18 PM, Fande Kong <fdkong.jd at gmail.com> wrote:
> 
> 
> 
> On Wed, Oct 12, 2016 at 9:12 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> > On Oct 12, 2016, at 10:00 PM, Amneet Bhalla <mail2amneet at gmail.com> wrote:
> >
> >
> >
> > On Monday, October 10, 2016, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > > On Oct 10, 2016, at 4:01 PM, Kong, Fande <fande.kong at inl.gov> wrote:
> > >
> > > Hi All,
> > >
> > > I know how to remove the null spaces from a singular system using creating a MatNullSpace and attaching it to Mat.
> > >
> > > I was really wondering what is the philosophy behind this? The exact algorithms we are using in PETSc right now?  Where we are dealing with this, preconditioner, linear solver, or nonlinear solver?
> >
> >    It is in the Krylov solver.
> >
> >    The idea is very simple. Say you have   a singular A with null space N (that all values Ny are in the null space of A. So N is tall and skinny) and you want to solve A x = b where b is in the range of A. This problem has an infinite number of solutions     Ny + x*  since A (Ny + x*) = ANy + Ax* = Ax* = b where x* is the "minimum norm solution; that is Ax* = b and x* has the smallest norm of all solutions.
> >
> >       With left preconditioning   B A x = B b GMRES, for example, normally computes the solution in the as alpha_1 Bb   + alpha_2 BABb + alpha_3 BABABAb + ....  but the B operator will likely introduce some component into the direction of the null space so as GMRES continues the "solution" computed will grow larger and larger with a large component in the null space of A. Hence we simply modify GMRES a tiny bit by building the solution from alpha_1 (I-N)Bb   + alpha_2 (I-N)BABb + alpha_3 (I-N)BABABAb + ....  that is we remove from each new direction anything in the direction of the null space. Hence the null space doesn't directly appear in the preconditioner, just in the KSP method.   If you attach a null space to the matrix, the KSP just automatically uses it to do the removal above.
> >
> > Barry, if identity matrix I is of size M x M (which is also the size of A) then are you augmenting N (size M x R; R < M) by zero colums to make I  - N possible? If so it means that only first R values of vector Bb are used for scaling zero Eigenvectors of A. Does this choice affect iteration count, meaning one can arbitrarily choose any R values of the vector Bb to scale zero eigenvectors of A?
> 
>   Yes I wasn't very precise here. I should have written it as as the projection of the vector onto the complement of the null space which I think can be written as I - N(N'N)^-1N'
> This is done with the routine MatNullSpaceRemove(). The basis of the null space you provide to MatNullSpaceCreate() should not effect the convergence of the Krylov method at all since the same null space is removed regardless of the basis.
> 
> I think we need to make sure that the basis vectors are orthogonal to each other and they are normalized. Right?

  Yes, MatNullSpaceCreate() should report an error if all the vectors are not orthonormal. 

#if defined(PETSC_USE_DEBUG)
  if (n) {
    PetscScalar *dots;
    for (i=0; i<n; i++) {
      PetscReal norm;
      ierr = VecNorm(vecs[i],NORM_2,&norm);CHKERRQ(ierr);
      if (PetscAbsReal(norm - 1.0) > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PetscObjectComm((PetscObject)vecs[i]),PETSC_ERR_ARG_WRONG,"Vector %D must have 2-norm of 1.0, it is %g",i,(double)norm);
    }
    if (has_cnst) {
      for (i=0; i<n; i++) {
        PetscScalar sum;
        ierr = VecSum(vecs[i],&sum);CHKERRQ(ierr);
        if (PetscAbsScalar(sum) > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PetscObjectComm((PetscObject)vecs[i]),PETSC_ERR_ARG_WRONG,"Vector %D must be orthogonal to constant vector, inner product is %g",i,(double)PetscAbsScalar(sum));
      }
    }
    ierr = PetscMalloc1(n-1,&dots);CHKERRQ(ierr);
    for (i=0; i<n-1; i++) {
      PetscInt j;
      ierr = VecMDot(vecs[i],n-i-1,vecs+i+1,dots);CHKERRQ(ierr);
      for (j=0;j<n-i-1;j++) {
        if (PetscAbsScalar(dots[j]) > PETSC_SQRT_MACHINE_EPSILON) SETERRQ3(PetscObjectComm((PetscObject)vecs[i]),PETSC_ERR_ARG_WRONG,"Vector %D must be orthogonal to vector %D, inner product is %g",i,i+j+1,(double)PetscAbsScalar(dots[j]));
      }
    }
    PetscFree(dots);CHKERRQ(ierr);
  }
#endif


> 
> Fande,
>  
> 
>   Barry
> 
> 
> >
> >     With right preconditioning the solution is built from alpha_1 b   + alpha_2 ABb + alpha_3 ABABb + .... and again we apply (I-N) to each term to remove any part that is in the null space of A.
> >
> >    Now consider the case   A y = b where b is NOT in the range of A. So the problem has no "true" solution, but one can find a least squares solution by rewriting b = b_par + b_perp where b_par is in the range of A and b_perp is orthogonal to the range of A and solve instead    A x = b_perp. If you provide a MatSetTransposeNullSpace() then KSP automatically uses it to remove b_perp from the right hand side before starting the KSP iterations.
> >
> >   The manual pages for MatNullSpaceAttach() and MatTranposeNullSpaceAttach() discuss this an explain how it relates to the fundamental theorem of linear algebra.
> >
> >   Note that for symmetric matrices the two null spaces are the same.
> >
> >   Barry
> >
> >
> >    A different note: This "trick" is not a "cure all" for a totally inappropriate preconditioner. For example if one uses for a preconditioner a direct (sparse or dense) solver or an ILU(k) one can end up with a very bad solver because the direct solver will likely produce a very small pivot at some point thus the triangular solver applied in the precondition can produce HUGE changes in the solution (that are not physical) and so the preconditioner basically produces garbage. On the other hand sometimes it works out ok.
> >
> >
> > >
> > >
> > > Fande Kong,
> >
> >
> >
> > --
> > --Amneet
> >
> >
> >
> >
> 
>