16
$\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$
12
$\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$
  • $\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 '12 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$ – Geoff Oxberry Jun 13 '12 at 15:12
  • $\begingroup$ Yes, I mean that. I will look at the permutation, thanks. $\endgroup$ – Alexander Jun 13 '12 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$ – Geoff Oxberry Jun 13 '12 at 15:20
8
$\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$
  • $\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$ – Geoff Oxberry Jun 13 '12 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 '12 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 '12 at 7:10
  • $\begingroup$ I edited my answer a little. $N$ is computed by the recipe of Geoff Oxberry. $\endgroup$ – Arnold Neumaier Jun 14 '12 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 '12 at 14:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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