Hi Pawel,
thanks for sharing your code for finding nullspace of large matrices. This appears very useful, and I am considering using it for my research in medical image reconstruction. There is just a small problem, with which I was hoping you could help us.
We are interested in "almost-rank-deficient" matrices rather than strictly rank-deficient. So we want to be able to supply a tolerance parameter which says, "consider all eigenvectors with eigenvalues smaller than this tolerance to belong to the nullspace". My question is: is it possible to add this feature to your code?
I would really appreciate your help!
sincerely,
Ashish Raj
Overall, a really useful and functional addition - thank you!
I agree with John D'Errico that some more detailed comments/help would have been useful. Also, there's some typos in the comments header. That's really the only reason this is not a 5.
However, since the script is not meant to find an orthogonal basis, I don't think see why that should be a negative point. As I understand it anyway, orthogonal bases tend to be dense anyway - if you really need an orthogonal basis of the null space you can use the lower part of Q^{T} from a QR factorization.
I would like to see it upgraded using Tim Davis' suggestion (although dmperm now uses the code from CSparse.)
Suggestion for the author: You may want to try using dmperm first. It will split the system into its structurallly over-determined, well-determined, and under-determined parts. That function doesn't make the "coarse decomposition" clear, however. It has contains a bug when used for large matrices. Use cs_dmperm in CSparse instead.
This fills a hole in Matlab. For sparse matrices, it should be much faster than null(full(...)). Why did I rate it as only a 4?
The help seems reasonable, and if you know what you are doing, it should work.
I don't terribly like the interface though, returning a cell array. You need to remember which element number to use to extract the arrays that you need. Far more logical, easy to use, and self commenting, would be to return a structure with appropriately named fields.
Next, if what you really need is an orthogonal basis for the nullspace, then spspaces will not serve this purpose.
Finally, I'd have preferred additional comments (examples) in the help so that a less experienced user could figure out how to use this code. As it is, that person is likely to be completely lost if all they know is they need to find a sparse basis for the nullspace of the columns of a sparse matrix.
