Rapidly determining whether or not a dense matrix is of low rank

Question

In a software project that I'm working on, certain computations are vastly easier for dense low-rank matrices. Some problem instances involve dense low-rank matrices, but they're given to me in full, rather than as factors, so I'll have to check the rank and factor the matrix if I want to take advantage of the low-rank structure.

The matrices in question are typically fully or nearly fully dense, with n ranging from one hundred up to a few thousand. If a matrix has low rank (say less than 5 to 10), then computing the SVD and using it form a low-rank factorization is worth the effort. However, if the matrix is not of low rank, then the effort would be wasted.

Thus I'd like to find a fast and reasonably reliable way of determining whether or not the rank is low before investing the effort to do a full SVD factorization. If at any point it becomes clear that the rank is above the cutoff, the process can stop immediately. If the procedure mistakenly declares the matrix to be of low rank when it isn't, this isn't a huge issue, since I'd still be doing a full SVD to confirm the low rank and find a low-rank factorization.

Options that I've considered include a rank revealing LU or QR factorization followed by a full SVD as the check. Are there other approaches that I should consider?

Answer 1

There is a neat trick I have recently learned from this paper. You start doing rank-revealing QR, and stop after the first $k$ Householder reflections, when you have a matrix of the form $$ \begin{bmatrix} R_1 & R_{12}\\ 0 & R_{22} \end{bmatrix}, $$ with $R_1$ triangular of size $k\times k$, and $R_{22}$ typically not triangular (since we stopped after the first $k$ iterations of our main loop). At this point, you check if $\|R_{22}\| \leq \varepsilon$: if it holds, then $A$ is at distance at most $\varepsilon$ from a matrix of rank $\leq k$; otherwise it shouldn't be (barring numerical errors).

This procedure costs $O(n^2k)$ for a dense $n\times n$ matrix.

Answer 2

The problem, of course, is that computing the true rank (e.g., via a QR decomposition) is not really any cheaper than computing a low-rank representation of the matrix.

The best you can probably do is to use a randomized algorithm to find low-rank approximations. These can, at least in theory, be significantly faster than working on the entire matrix because, in essence, they only compute decompositions for projections of the matrix onto random subspaces.

Whether that's worth it for a matrix of size $100\times 100$ may be a good question, but if your problems really become large, I would suspect that it pays off.

Answer 3

Another approach worth trying is to use Adaptive Cross Approximation (ACA). It is a pretty popular algorithm that has many implementations available online. For the reference, you can see the original paper:

• S. A. Goreinov, E. E. Tyrtyshnikov, N. L. Zamarashkin, "A theory of pseudoskeleton approximations," Linear Algebra Appl., vol. 261, no. 1–3, pp. 1–21, Aug. 1997.

ACA and its variations (say, ACA+, hybrid cross approximation HCA) can be used in different scenarios. You, already having the whole dense matrix computed is one of the favorable, as you will be able to calculate residuals exactly if needed.

If heuristical residuals (see the algorithm) suffice, I believe your complexity will be $\mathcal O(Nr)$, where $N$ is the size of the square matrix and $r(\epsilon)$ is the rank. Note, that rank $r$ is a function of the prescribed truncation tolerance $\epsilon$. While the exact and guaranteed error bounds will require $\mathcal O(N^2r)$.

Answer 4

For the simple case where the matrix $A$ is symmetric positive-definite, calculate its say 20 biggest eigenvalues, and see if they're $\to 0$, or compare norms. ARPACK is fast for this; more important, it needs only a function $x \to A \, x$ . So for general $A$, look at the eigenvalues of $A^T A$ (as a LinOp, without instantiating it.)

scipy.sparse.linalg.svds does this: LinOp$( A^T A ) \to$ Arpack, for $A$ of any size:

from scipy.sparse.linalg import svds
sing = svds( A, k=20, tol=1e-4, return_singular_vectors=False )  # v0=random
# runtimes on random-normal n x n:
# n = 100, 1k, 2k
#       5, 130, 770 ms