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

[Fande Kong]

On Wed, Oct 12, 2016 at 9: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?
>

I think it is just a notation. Let us denote Bb as y, that is, y = Bb. N =
{x_0, x_1, ..., x_r}, where x_i is a vector. Applying  I-N to y means  y =
y - \sum (y, x_i) x_i.  (y, x_i) is the inner product of vectors y and
x_i.  Look at this code for more details,
http://www.mcs.anl.gov/petsc/petsc-current/src/mat/interface/matnull.c.html#MatNullSpaceRemove

Fande,



>
>>     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
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20161012/a9b1ab46/attachment.html>