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

[Barry Smith]

  Fande,

     Could you send me (petsc-maint at mcs.anl.gov) a non symmetric matrix you have that has a different null space for A and A'. This would be one that is failing with right preconditioning. Smaller the better but whatever size you have. Run the code with -ksp_view_mat binary and send the resulting file called binaryoutput.

   I need a test matrix to update the PETSc code for this case.


   Barry

> On Oct 11, 2016, at 3:04 PM, Kong, Fande <fande.kong at inl.gov> wrote:
> 
> 
> 
> On Tue, Oct 11, 2016 at 12:18 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> > On Oct 11, 2016, at 12:01 PM, Kong, Fande <fande.kong at inl.gov> wrote:
> >
> >
> >
> > On Tue, Oct 11, 2016 at 10:39 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > > On Oct 11, 2016, at 9:33 AM, Kong, Fande <fande.kong at inl.gov> wrote:
> > >
> > > Barry, Thanks so much for your explanation. It helps me a lot.
> > >
> > > On Mon, Oct 10, 2016 at 4:00 PM, 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
> > >
> > >  Does "I" mean an identity matrix? Could you possibly send me a link for this GMRES implementation, that is, how PETSc does this in the actual code?
> >
> >    Yes.
> >
> >     It is in the helper routine KSP_PCApplyBAorAB()
> > #undef __FUNCT__
> > #define __FUNCT__ "KSP_PCApplyBAorAB"
> > PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
> > {
> >   PetscErrorCode ierr;
> >   PetscFunctionBegin;
> >   if (!ksp->transpose_solve) {
> >     ierr = PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr);
> >     ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr);
> >   } else {
> >     ierr = PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr);
> >   }
> >   PetscFunctionReturn(0);
> > }
> >
> >
> > PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y)
> > {
> >   PetscErrorCode ierr;
> >   PetscFunctionBegin;
> >   if (ksp->pc_side == PC_LEFT) {
> >     Mat          A;
> >     MatNullSpace nullsp;
> >     ierr = PCGetOperators(ksp->pc,&A,NULL);CHKERRQ(ierr);
> >     ierr = MatGetNullSpace(A,&nullsp);CHKERRQ(ierr);
> >     if (nullsp) {
> >       ierr = MatNullSpaceRemove(nullsp,y);CHKERRQ(ierr);
> >     }
> >   }
> >   PetscFunctionReturn(0);
> > }
> >
> > "ksp->pc_side == PC_LEFT" deals with the left preconditioning Krylov methods only? How about the right preconditioning ones? Are  they just magically right for the right preconditioning Krylov methods?
> 
>    This is a good question. I am working on a branch now where I will add some more comprehensive testing of the various cases and fix anything that comes up.
> 
>    Were you having trouble with ASM and bjacobi only for right preconditioning?
> 
> 
> Yes. ASM and bjacobi works fine for left preconditioning NOT for RIGHT preconditioning. bjacobi converges, but produces a wrong solution. ASM needs more iterations, however the solution is right. 
> 
> 
>  
>     Note that when A is symmetric the range of A is orthogonal to null space of A so yes I think in that case it is just "magically right" but if A is not symmetric then I don't think it is "magically right". I'll work on it.
> 
> 
>    Barry
> 
> >
> > Fande Kong,
> >
> >
> > There is no code directly in the GMRES or other methods.
> >
> > >
> > > (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.
> > >
> > >     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.
> > >
> > > What preconditioners  are appropriate? asm, bjacobi, amg? I have an example which shows  lu and ilu indeed work, but asm and bjacobi do not at all. That is why I am asking questions about algorithms. I am trying to figure out a default preconditioner for several singular systems.
> >
> >    Hmm, normally asm and bjacobi would be fine with this unless one or more of the subblocks are themselves singular (which normally won't happen). AMG can also work find sometimes.
> >
> >    Can you send a sample code?
> >
> >   Barry
> >
> > >
> > > Thanks again.
> > >
> > >
> > > Fande Kong,
> > >
> > >
> > >
> > > >
> > > >
> > > > Fande Kong,
> 
>