20
$\begingroup$

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.

$\endgroup$

3 Answers 3

17
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Alexander
    Jun 13, 2012 at 15:04
  • $\begingroup$ 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? $\endgroup$ Jun 13, 2012 at 15:12
  • $\begingroup$ Yes, I mean that. I will look at the permutation, thanks. $\endgroup$
    – Alexander
    Jun 13, 2012 at 15:14
  • $\begingroup$ 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. $\endgroup$ Jun 13, 2012 at 15:20
  • $\begingroup$ 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. $\endgroup$
    – Akh
    Apr 20, 2020 at 20:13
10
$\begingroup$

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.

$\endgroup$
12
  • $\begingroup$ +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. $\endgroup$ Jun 13, 2012 at 15:34
  • $\begingroup$ Indeed, it sounds reasonable, however my matrices fit into 16 GB of RAM, so I would stay with standard matlab qr. $\endgroup$
    – Alexander
    Jun 13, 2012 at 16:15
  • $\begingroup$ 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. $\endgroup$
    – Alexander
    Jun 14, 2012 at 7:10
  • $\begingroup$ I edited my answer a little. $N$ is computed by the recipe of Geoff Oxberry. $\endgroup$ Jun 14, 2012 at 10:10
  • $\begingroup$ 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? $\endgroup$
    – Alexander
    Jun 14, 2012 at 14:30
0
$\begingroup$

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/

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.