20

Given a dense matrix $$A \in R^{m \times n}, m >> n; max(m) \approx 100000 $$ what is the best way to find its null-space basis within some tolerance $\epsilon$?

Based on that basis can I then say that certain cols are linearly dependent within $\epsilon$? In other words, having null space basis computed, what columns of $A$ have to be removed in order to get nonsingular matrix?

References are appreciated.

Alexander
  • 1,111
  • 1
  • 8
  • 15

3 Answers3

17

Standard methods for determining the null space of a matrix are to use a QR decomposition or an SVD. If accuracy is paramount, the SVD is preferred; the QR decomposition is faster.

Using the SVD, if $A = U\Sigma V^{H}$, then columns of $V$ corresponding to small singular values (i.e., small diagonal entries of $\Sigma$) make up the a basis for the null space. The relevant tolerance here is what one considers a "small" singular value. MATLAB, for instance, takes small to be $\max(m,n) \cdot \varepsilon$, where $\varepsilon$ is related to machine accuracy (see here in MATLAB's documentation).

Using the QR decomposition, if $A^{T} = QR$, and the rank of $A$ is $r$, then the last $n-r$ columns of $Q$ make up the nullspace of $A$, assuming that the QR decomposition is rank revealing. To determine $r$, calculate the number of entries on the main diagonal of $R$ whose magnitude exceeds a tolerance (similar to that used in the SVD approach).

Don't use LU decomposition. In exact arithmetic, it is a viable approach, but with floating point arithmetic, the accumulation of numerical errors makes it inaccurate.

Wikipedia covers these topics here.

Geoff Oxberry
  • 30,394
  • 9
  • 64
  • 127
  • Geoff, talking in terms of QR, suppose I have the decomposition, how do I then relate null space basis and columns in original matrix? In other words, what columns should I remove from $A$ in order to get rid of null space? The point here is to work with $A$ itself and not with its decomposition. – Alexander Jun 13 '12 at 15:04
  • Routines that calculate the QR decomposition normally include an option to return a permutation vector indicating how columns are permuted to obtain the QR factorization. The last $n-r$ entries of that permutation vector would correspond to the rows of $A$ (columns of $A^{T}$) that are in the nullspace. The first $r$ entries of that vector correspond to the columns of $A^{T}$ that are linearly independent. I'm not sure what you mean by "get rid of the null space". Do you mean you want to remove columns of $A$ to obtain a nonsingular matrix? – Geoff Oxberry Jun 13 '12 at 15:12
  • Yes, I mean that. I will look at the permutation, thanks. – Alexander Jun 13 '12 at 15:14
  • That is a different question. What you would then do instead is calculate the QR decomposition (or SVD) of $A$. If you calculate the QR decomposition of $A$, you can calculate the rank of $A$ as in the answer above (no need to transpose the matrix), and then the first $r$ entries (where $r$ is the rank of $A$) of the permutation vector correspond to the independent columns of $A$. The same sort of algorithm applies to the SVD; if you can return a permutation vector along with the decomposition, that should provide the necessary information. – Geoff Oxberry Jun 13 '12 at 15:20
  • From my own experiments, the null space vectors appear to be in the columns of Q that correspond to the zeros of R, not necessarily the last n - r columns of Q. – Akh Apr 20 '20 at 20:13
  • Note that Matlab’s QR is not rank-revealing. You need something like this: https://www.mathworks.com/matlabcentral/fileexchange/18591-rank-revealing-qr-factorization – Amit Hochman Sep 23 '22 at 06:59
10

If $m\gg n$, as your question indicates, you can save some work by first picking an index set $I$ of $p\approx 5n$ (say) random rows and using the orthogonal factorization $A_{I:}^T=QR$. (The QR-factorization is the one where $Q$ is sqare and $R$ is rectangular of rank $r$, and the remaining $n-r$ columns of $R$ are zero. Using a permuted QR factorization will enhance stability; the permutation must then be accounted for in a more detailed recipe.)

Typically, this will give you a much lower dimensional subspace spanned by the columns of $N$, the last $n-r$ columns of $Q$. This subspace contains the null space of $A$. Now pick another, disjoint random index set and compute the QR factorization of $(A_{I:}N)^T$. Multiply the resulting null space on the left by $N$ to get an improved $N$ of probably even lower dimension. Iterate until the dimension of $N$ no longer decreases. Then you probably have the correct null space and can check by computing $AN$. If this is not yet negligible, do further iterations with the most significant rows.

Edit: Once you have $N$, you can find a maximal set $J$ of linearly independent columns of $A$ by an orthogonal factorization of $N^T=QR$ with pivoting. Indeed, the set $J$ of indices not chosen as pivots will have this property.

Arnold Neumaier
  • 11,318
  • 20
  • 47
  • +1 for an efficient way to determine the nullspace of a large matrix. I'll have to remember to consult this answer later on when I need it. – Geoff Oxberry Jun 13 '12 at 15:34
  • Indeed, it sounds reasonable, however my matrices fit into 16 GB of RAM, so I would stay with standard matlab qr. – Alexander Jun 13 '12 at 16:15
  • Prof. Neumaier, I've decided to test that algorithm, but I don't understand exactly what is $N$ and what does "compute the QR factorization of $(A_{I:}N)^T$" mean? Could you please explain a bit more. – Alexander Jun 14 '12 at 07:10
  • I edited my answer a little. $N$ is computed by the recipe of Geoff Oxberry. – Arnold Neumaier Jun 14 '12 at 10:10
  • Thank you. I implemented it. However, as far as I see, this algorithm doesn't allow to define me a set of linearly independent columns of $A$ (since we decompose $A_{I:}^T$ rather than $A_{I:}$), but just helps to estimate nullspace basis itself? – Alexander Jun 14 '12 at 14:30
  • I added a corresponding recipe to my answer. – Arnold Neumaier Jun 14 '12 at 14:44
  • Makes sense. Although I can't implement it with standard Matlab $[Q,R,E] = qr(N');$ All columns of E contain 1s and I can't determine pivots. – Alexander Jun 14 '12 at 16:11
  • @Alexander:[Q,R,E] = QR(N',0) gives the wanted information. – Arnold Neumaier Jun 14 '12 at 16:14
  • Is there a paper that proves that this approach convergences? I can see how with noisy data it an iteration could cause it to move in the wrong direction but on average I think it should step towards the correct solution. – lessthanoptimal Mar 26 '18 at 16:36
  • @PeterAbeles: With exact data, the process will produce a basis after finitely many steps if all index sets are chosen disjoint,With noisy data it is not clear what the precise goal is, hence it is difficult to discuss convergence issues. But in this case, choosing the index sets randomly should quickly produce a very acceptable result. – Arnold Neumaier Mar 26 '18 at 18:46
  • @Alexander, I have a similar problem for my application. Could you implement Prof. Neumaier's suggestion? For my specific application, I have a different implementation in MATLAB but the main issue is that the null space is not accurate up to roundoff. I used this implementation for different problems with a far higher accuracy. But for a problem involving some new transformation matrices, the accuracy is not good enough. For instance the norm of AN is on the order of 1e-7, well, small but still not at the roundoff as mentioned – Umut Tabak Dec 05 '23 at 09:57
  • Particularly, I implemented the method discussed in this link: https://mathoverflow.net/questions/253995/how-can-one-construct-a-sparse-null-space-basis-using-recursive-lu-decomposition/253997?noredirect=1#comment672418_253997 – Umut Tabak Dec 05 '23 at 10:12
0

Above are answers to all of your questions except the last one: "In other words, having null space basis computed, what columns of have to be removed in order to get nonsingular matrix?"

You can use dimension-reduction to reconstruct A using factors of reduced elements of SVD, but instead of using SVD, you may prefer CUR decomposition for very large A.

Calculate U,S,V from SVD(A), then reduced A = U2 * S2 * V2^T where U2 is U[mxr], r is the rank of A, S2 = S[rxr], and V2^T is V^T[rxn]. References are chapter 11 of "Mining of Massive Datasets" by Jure Leskovec, Anand Rajaraman, and Jeff Ullman http://www.mmds.org/

nichole
  • 1