[petsc-users] Best way to compute the null space of a sparse maximum-rank rectangular matrix

[Jose E. Roman]

El 08/03/2014, a las 20:42, Luca Argenti escribió:

> Dear all,
> 
>       I need to evaluate the null space Span{v_i} of a sparse rectangular matrix A \in C^{ N x (N+k)}, 
> 
>                    A v_i = 0,            (1)
> 
> where 
> 
> i) N is typically very big ( N ~ 500’000 ) and k, in comparison, is very small  ( k ~ 500 ).
> ii) The sub-matrix A_ij  i,j<= N is hermitian and non-singular. Equation (1), therefore, has exactly k solutions. 
> iii)The matrix is sparse, with a fill typically <= 3%, and its columns/rows can be reordered in such a way that 
>     a very large block, A_ij with i,j > n, & i,j<=N, n << N, is band-diagonal. 
> iv) A has a dominant diagonal. 
> v) For large values of i,j, the number of non-zero diagonals in the central band drop by about an order of magnitude.
> vi) Finally, this problem must be solved for several (thousands) closely-spaced values of an external parameter 
>      Q on which A depends continuously, A = A(Q). Most of the time, therefore, the null space at Q_{i+1} is arguably
>      very close to the null-space at Q_i .
> 
> My feeling is that this problem is very well defined, and that a parallel sparse iterative method should be
> able to solve it with no issues or unnecessary operations. Yet, probably because I am not an expert 
> of either PETSc or SLEPc, the two libraries I have considered so far, all the possible solutions that I found 
> seem to provide much more information than needed (thus, consuming much more resources than
> warranted). For example: is it really necessary to make a sparse LU factorization for the *whole* matrix? 
> In practice, one is looking for the null eigenspace of  A^h A. However, SLEPc suggests that this operation is 
> much more expensive than for a sparse A matrix alone (is it so? Shouldn’t Lanczos be implementable at just 
> twice the cost?), or maybe I misinterpreted the user guide.
> 
> Your suggestions will be greatly appreciated. Thank you so much for your help!
> 
> Cheers,
> 
>       Luca
> 

The nullspace of A^h A is equal to the right singular space of A corresponding to the zero singular value. It should be possible to compute this with SLEPc's SVD. Computing a large number of zeros may be problematic, so I cannot say in advance if the method will succeed. If you can generate a small matrix with these properties, send it to my personal address (not the list) and I will give it a try.

Jose