4
$\begingroup$

I make some numerical experiments, involving rank of integer valued matrices of the size about $14\times 24$. As the matrix is integer valued, theoretically there should be no room for errors. However numpy does it using SVD and the results become unstable, depending on the SVD cutoff.

I wonder if there are fast stable methods for computing the rank of integer-valued matrices.

$\endgroup$
7
  • $\begingroup$ I believe you can compute the rank modulo a random, fairly small, prime, or even try a few different ones $\endgroup$ Commented May 11 at 14:18
  • $\begingroup$ can you elaborate a bit, please? $\endgroup$ Commented May 11 at 14:21
  • 1
    $\begingroup$ Possible duplicate of this question. Also, since you mention that you're using NumPy, you might want to look into using SymPy or SageMath, which offer symbolic (instead of numerical) linear algebra. Or are you specifically asking about numerical algorithms? $\endgroup$ Commented May 12 at 6:00
  • 1
    $\begingroup$ By "symbolic", I just mean that the computations are exact (see symbolic algebra system). These systems allow you to do exact linear algebra over various rings, including $\mathbb Z$ and $\mathbb Q$. The intermediate results will be stored in arbitrary precision integers or fractions, instead of being rounded to floating point. Instead of worrying about numerical stability, you now have to worry about the entries becoming too large, negatively impacting the space and time complexity. These issues are addressed in the question I linked. $\endgroup$ Commented May 12 at 23:09
  • 2
    $\begingroup$ Does this answer your question? Complexity of computing matrix rank over integers $\endgroup$ Commented May 13 at 2:18

1 Answer 1

7
$\begingroup$

Note that if you take a prime $p$ and treat the matrix $A$ as a matrix $A'$ over $\mathbb{Z} / p \mathbb{Z}$, then from the property of $\operatorname{rank}(A)$ being the largest order of a non-zero minor in $A$, and the fact that treating a matrix over $\mathbb{Z}$ as a matrix over $\mathbb{Z} / p \mathbb{Z}$ does the same to the determinant, we get that $\operatorname{rank}(A') \leq \operatorname{rank}(A)$, and if $p$ doesn't divide a particular determinant (which depends on $A$) then the ranks are equal. Therefore, we can try a few random $p$s, and take the maximum of $\operatorname{rank}(A')$.

This is obviously numerically OK, as we only need to work with integers mod $p$, and requires $O(T n^3)$ operations, where $T$ is the number of tries.

Let's estimate the number of required trials - if the maximum absolute value in an $n \times n$ matrix is $m$, by Hadamard's inequality the determinant is at most $m^n n^{\frac n2}$. From this, we can derive a trivial upper bound for its number of different prime divisors - $n (\log_2(m) + \frac{\log_2(n)}2)$, so if we sample uniformly from a sufficiently large distribution, which can be efficiently in practice by just doing a primality test on random values, we only need $O(k)$ tries to get a error rate of $2^{-k}$. This requires numbers with size $O(\log n + \log \log m)$, which is usually feasible.

$\endgroup$
4
  • $\begingroup$ Thank you! Sounds reasonable, I will try to implement it now. $\endgroup$ Commented May 11 at 14:57
  • $\begingroup$ The galois package extends numpy to implement fast linear algebra (including finding the rank) over finite fields, so it should be easy to implement this (galois also has some functions for generating random primes). $\endgroup$ Commented May 12 at 0:28
  • $\begingroup$ Thanks! @Tuomas, I will look at galois package $\endgroup$ Commented May 12 at 15:22
  • $\begingroup$ @TuomasLaakkonen, in fact, galois package did the job very nicely. Thanks for the recommendation! $\endgroup$ Commented Jun 2 at 0:40

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .