13
$\begingroup$

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?

| cite | improve this question | |
$\endgroup$
8
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ This is essentially the approach that I described in the question. I think that Wolfgang Bangerth's proposed answer could do better than $O(n^{2}k)$. $\endgroup$ – Brian Borchers Jul 12 '17 at 22:39
7
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ From what I know of these algorithms, they produce a low-rank matrix that is reasonably close in norm to the given matrix. I need to know whether or not there's a (for example) rank-10 or less matrix that is very close to the given matrix (say a relative error of 1.0e-10 or better.) $\endgroup$ – Brian Borchers Jul 8 '17 at 4:03
  • $\begingroup$ Yes, but you can also do a QR decomposition of the projected (low-dimensional) matrix and if that decomposition reveals a lack of full rank, then you will also have a rank-deficient original matrix. Wasn't that the criterion you needed to do a QR decomposition on the original matrix? $\endgroup$ – Wolfgang Bangerth Jul 10 '17 at 2:55
  • $\begingroup$ I can see that the rank of the projected matrix is less than or equal to $k$ (the number of rows in the random matrix I multiply times A) and the rank of A. If it is of rank $k$, then the original matrix cannot be of rank $k-1$ or less. If it is of rank less than $k$ then I could have just been unlucky or $A$ was of rank less than $k$. Finding the rank of the $k$ by $n$ matrix can be done in $O(k^{2}n)$ time. However, if the random matrix that I multiply times $A$ is dense, mutiplication takes $O(kn^{2})$ time. Are there sparse matrices that preserve the rank with high probability? $\endgroup$ – Brian Borchers Jul 12 '17 at 22:36
  • $\begingroup$ I don't know. I agree (and meant to imply) that the algorithm can only tell you if a matrix is not of full rank. It cannot tell you if the matrix is of full rank unless you take all $k=n$ random directions. My hope simply would be that you get an answer for sufficiently small $k$ where $kn^2\ll n^3$. $\endgroup$ – Wolfgang Bangerth Jul 19 '17 at 3:03
1
$\begingroup$

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:

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)$.

| cite | improve this answer | |
$\endgroup$
0
$\begingroup$

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
| cite | improve this answer | |
$\endgroup$

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.