I have a large sparse matrix A which is column rank-defficient. Typical size of A is over 100000x100000. In my computation, I need the matrix W whose columns span the null space of A. But I do not know how to fastly compute all the columns of W.
If A is small-scale, I know there are several numerical methods based on matrix factorization, such as LU, QR, SVD. But for large-scale matrices, I can not find an efficient iterative method to do this.
Could you please give me some help?
Haibo Li, I'm sorry, but I may not be able to help you much with this problem because I don't have your data (You can upload it here or perhaps share with me). However, I can tell you that it is possible to compute the null space of a large sparse matrix with Python's NumPy and SciPy libraries. The size of your matrix, 100,000x100,000, is very large. You will need a computer with a lot of RAM and processing power to compute the null space of this matrix. For any kind of assistance you can discuss with me for a prompt response on my WA; https://wa.me/+923440907874
This paper: http://cedric.cnam.fr/~bentzc/INITREC/Files/DW10.pdf might be a good starting point.
This: https://math.mit.edu/~gs/learningfromdata/lafe019 might, also, be useful.
I suggest inverse iterations for the matrix B = alpha I + A' A, where alpha is a small positive number (e.g. 1e-6 times the square of the average value of the nonzero elements of A, see also the Note below), I is the identity matrix and A' is the transpose of A.
1. Take a random matrix U_0 of size n rows by m columns, where n is the number of columns in A and m is your initial guess for the size of the null space of A. Compute the SVD for U_0, i.e. U_0 = U Sigma V', and take the columns of W_0 = U as initial approximations to vectors in the null space of A.
2: Solve m linear systems to find the columns of a matrix U_1 such that B U_1 = W_0. As this step is going to be repeated many times, it is a good idea to factorize B, which should not be too difficult since B is sparse and symmetric positive definite (consider e.g. scipy.sparse.linalg.splu).
3: Compute the SVD U_1 = U Sigma V' and take the columns of W_1 = U as the next approximations to vectors in the null space of A.
4: Repeat steps 2 and 3 (B U_2 = W_1, U_2 = W_1 Sigma V' etc. etc.) until the largest singular values in Sigma get very close to 1/alpha.
5: If the number of computed columns of W_i corresponding to singular values (that are very close to) 1/alpha is less than m, and other singular values in Sigma are noticeably smaller than 1/alpha, you are done - the computed columns of W_i corresponding to the singular values 1/alpha are the sought orthonormal basis vectors of the null space of A. Otherwise, repeat steps 1, 2, 3, 2, 3, ... to compute more vectors in the null space but this time orthogonalize the columns of U_* to ALL already computed approximations to the vectors in the null space of A before the SVD (you may need to repeat this procedure again and again until you get all the vectors you need).
Note: the smaller alpha the faster the convergence of W_i, so try alpha as small as the factorization of B allows (it will fail if alpha is too small).
HTH
Evgueni
- Badges
- Science topic
- Similar topics
- Indexes