How do computer algebra packages like Sagemath implement rank of a matrix

Question

I am not sure if this is the right place to ask this question, but I believe there will be people here who do computations on computer algebra packages like Sage in their work. I have been using Sagemath to perform some matrix rank computations. It turns up a few bizarre results occasionally.

For example, I had to find the rank of a matrix ( $100 \times 150$ ) with large integer entries (entries of magnitude in the range of $1$ to $10^{15}$ ). When I wrote the code with the matrix M declared as matrix(ZZ, R, C), or as matrix(QQ, R, C), it returns a rank of around 90 (which I believe is correct), whereas if I declare the matrix over the reals as matrix(RR, R, C), it returns a rank of around 50, which I believe is too low based on some conjectures I have.

So, overall I am curious, what are the standard way(s) to implement rank computation (and does it differ based on reals, or rationals) and where can I read more about these? If these issues arise due to precision errors, how can I get around them? And more importantly, how do I know beforehand that my computation is susceptible to precision errors? (I tried looking up the source code of Sagemath a couple of times, but I was quickly lost, so I hope someone can point me to the precise documentation/source code)

I think that almost all rank calculations over $\mathbb{R}$ will be subject to error. I think that the standard advice is to compute singular values instead. If you see some very small singular values, then a clear gap, then some larger ones, then you can guess that the very small values should really be zero. Of course, you will not be able to prove anything rigorously by this method alone. It may be that Sage is doing something like this internally and applying an inappropriate cutoff. — Neil Strickland, Nov 28 '17 at 12:29
Three axioms of linear algebra: (1) Don’t compute the inverse of a matrix (unless it’s unitary) (2) Don’t multiply two matrices. (3) Never compute the rank of a matrix over the reals. — Chris Godsil, Nov 28 '17 at 13:02
The condition number of the matrix will give you some sense of whether to expect numerical difficulties https://en.wikipedia.org/wiki/Condition_number#Matrices It sounds like you have some entries that are 15 orders of magnitude larger than some other entries, and that is also usually bad news. — Timothy Chow, Nov 28 '17 at 17:11
@NeilStrickland But can I trust the rank computations that Sage returns over $\mathbb{Q}$ ? Do they employ different algorithms? I was initially coding with Matlab/Octave and read somewhere that the rank it returns is the number of singular values above a certain threshold and that was when I switched to Sage. — Nikhil, Nov 28 '17 at 17:13
@Nikhil: Multiplying two matrices is expensive, and if they are sparse there are often better approaches. I have seen students and postdocs hung up on this very point in the past year. (Of course, these axioms are partly tongue-in-cheek, but I think they are useful to keep in mind.) — Chris Godsil, Nov 28 '17 at 17:43
@Nikhil : I do not know for sure, but I would guess that over $\mathbb{Q}$ , Sage uses exact arithmetic. Whether that means you should trust the computation is a philosophical question. Should we ever trust a computer program? Should we trust computations that we do by hand? — Timothy Chow, Nov 28 '17 at 18:11
@TimothyChow Yes, trust is not the right word. I kinda assumed always that these are easy tasks (heck, they are well with in Polynomial time!) and that all these packages probably implemented some clever version of Gaussian elimination. Also $10^{15}$ didn't seem too big either. So this whole excercise kinda came as a surprise -- I agree I was too naive to expect theory to work super well in practice! :-) — Nikhil, Nov 28 '17 at 18:25
Gaussian elimination with integer matrices can produce huge intermediate results already for 100x150 matrices. That can sometimes be avoided with lattice reduction tricks. It can also be avoided in practice by working modulo a few "random" large primes; the resulting upper bound on the rank might in principle not be sharp but if you want a rigorous proof you can probably build on mod-p results. (Often such commands offer a flag that lets you choose between a heuristic answer and an answer that comes with a proof but takes longer to compute.) — Noam D. Elkies, Nov 28 '17 at 18:30
@FedericoPoloni: I think (4) might be “Don’t write your code for standard operations”. — Chris Godsil, Nov 28 '17 at 19:51
@NoamD.Elkies Can you elaborate a bit on these lattice reduction tricks? I understand the computation mod-p (which Igor alludes to in his answer? or rather the answer in the mathoverlfow question he linked), but for the record, I'm more interested in a proof that takes longer to compute, but is correct! — Nikhil, Nov 28 '17 at 21:20
Suppose $A$ is a $100 \times 150$ integer matrix with all entries of size at most $H$ (you gave $H = 10^{15}$ ). Computing modulo a few primes you find that the rank is at least $90$ , and guess that it's exactly $90$ . So, regarding $A$ as a map ${\bf Z}^{150} \to {\bf Z}^{100}$ , you need $60$ independent vectors in $\ker A$ . You expect their entries to have size about $H^{2.5}$ (the $2.5$ arises as $150 / (150-60)$ , cf. Siegel's lemma). Choose $M$ sufficiently larger than $H^{2.5}$ , and let $L \subset {\bf Z}^{100+150}$ be the lattice consisting of all vectors $(MTa,a)$ [cont'd] — Noam D. Elkies, Nov 29 '17 at 00:21
. . . with $a \in {\bf Z}^{100}$ . Then a reduced basis for this lattice should start with $60$ vectors $(0,a)$ with $Ta=0$ , because any vector $(MTa,a)$ with $a$ not in $\ker T$ has size at least $M$ .
Once you've found a basis containing $60$ such vectors, you're done because $T$ has rank at most $150-60$ . Note that even if $M$ must be taken as large as $10^{100} > H^6$ this is way better than Gaussian elimination which would have you doing arithmetic with (ratios of) integers of size $H^{90}$ . — Noam D. Elkies, Nov 29 '17 at 00:26

Igor Rivin · Accepted Answer · 2017-11-28T18:39:34.440

I don't know what algorithm Sage actually uses, but computing rank over the integers is fun and easy: Complexity of computing matrix rank over integers . It is NOT so easy if you want good running time, and for that there are a number of papers of Arne Storjohann, which show that one can do it asymptotically as fast as for real matrices (a surprising result, in view of coefficient blow-up). Storjohann has actually implemented his algorithms, and I am sure Sage uses this or something like it.

Chen, Zhuliang; Storjohann, Arne, A BLAS based C library for exact linear algebra on integer matrices, Kauers, Manuel (ed.), Proceedings of the 2005 international symposium on symbolic and algebraic computation, ISSAC’05, Beijing, China, July 24--27, 2005. New York, NY: ACM Press (ISBN 1-59593-095-7). 92-99 (2005). ZBL1360.65086.

As for the reals, the fastest way to compute rank is to compute the singular value decomposition, and throw away singular values below some cutoff (probably around $10^{-6}$ ). As pointed out by many in the comments, rank is very unstable, so unless you want the "pseudo-rank" (as above, defined by the size of the singular values), don't go there.

Can I access their source code somewhere? I could only find links to the paper, but not the source code. — Nikhil, Nov 28 '17 at 21:30

How do computer algebra packages like Sagemath implement rank of a matrix

1 Answers1