Finding Null space of a large sparse matrix in MATLAB

4
votes

Over the course of finite difference discretization of an elliptic equation and application of Neumann BCs on all sides of the 2D domain I have a large sparse matrix. I need to find null space of its transpose to enforce consistency condition on either side. For a computational domain of 50X50 and 100X100 I am able to do this with the available RAM of 32 GB in Mathematica and MATLAB with the available full matrix commands of NullSpace and null respectively.

If the computational domain is 500X250 (which is of the order of what I would have in general) the RAM required for the storage of the matrix of size (500X250)X(500X250) is 125 GB and is highly prohibitive. I use sparse matrices to store this super-matrix and I don't have have the space constraint anymore. But I can't use the 'null" command on this as it is for only full matrices. MATLAB suggests to use "SVDS" command on sparse matrices. SVDS(A) gives only the first 6 singular values and singular vectors. There is another command SVDS(A,k,sigma) which gives "k" singular values and vectors around the scalar singular value of "sigma". When I use sigma=0 so as to find the singular vector corresponding to the the "zero" value, which would be the vector from null space basis I get an error that the "sigma" is too close to the eigen value of the matrix.

My matrix is singular and hence one of its Eigen values is "zero". How do I circumvent this error? Or Is there a better way to do this with the tools at hand. I have both MATLAB and Mathematica at hand.

Thanks in advance for your help and suggestions.

Best Trinath

1
matlabmatrixsparse-matrix
I am not an expert, but please ask yourself the question whether there is another way to achieve what you want. Are all of these 3 goals (solve your problem, enforce consistency, calculate nullspace) really neccesary? On a semi-serious note, if you are working with 64 bit the matrix size should not be a problem, and you may even be able to fit it into ram. – Dennis Jaheruddin
Thank you Dennis for a prompt reply. If I use sparse routines I circumvent the space issue in MATLAB. My question is regarding using the command "SVDS" to get a vector from the null space. I too am exploring other ways to do this and constantly asking if there is a better way to solve the problem. The Intel Poisson solver does not solve the all Neumann boundary problem for elliptic equation. I am also looking at exploiting the symmetric patterns in the matrix. Thanks – Trinath

1 Answers

1
votes

I guess you can try with some sort of decomposition.

http://www.mathworks.co.uk/matlabcentral/fileexchange/11120-null-space-of-a-sparse-matrix

Have you tried this?

Or maybe this?

http://www.mathworks.co.uk/matlabcentral/newsreader/view_thread/249467

I am confident they should work, but I have not tried them myself. Another way of proceeding would be to get into QR decomposition (which would give you the permutation of the first k independent columns, if k is the rank of your matrix. Then the vectors from k+1 to n would provide a basis for your null space).

Hope this helps.

Cheers,

GL