3
$\begingroup$

I need to determine ranks of the powers of a triangular $0$-diagonal integer matrix of a high dimension. Doing this I noted that the time Mathematica needs to determine MatrixRank grows initially linearly with the power of the matrix.

Example:

n=100;M=Table[If[j>i,RandomInteger[9],0],{i,n},{j,n}];m=IdentityMatrix[n];
r=Table[{AbsoluteTiming[m=m.M][[1]],AbsoluteTiming[MatrixRank[m]][[1]]},{i,50}]

{{0.000301,0.005631},{0.000245,0.014556},{0.000123,0.027389},{0.000119,0.041173},{0.000135,0.061031},{0.000174,0.079355},{0.000145,0.091159},{0.000687,0.104788},{0.000745,0.117567},{0.006144,0.139681},{0.004766,0.206832},{0.004758,0.222089},{0.004867,0.240039},{0.004661,0.234097},{0.003819,0.258007},{0.006516,0.256},{0.007784,0.278658},{0.006997,0.291819},{0.00692,0.301719},{0.00704,0.316412},{0.006978,0.339455},{0.006344,0.338896},{0.010009,0.447298},{0.007158,0.370554},{0.006502,0.431556},{0.005816,0.49495},{0.008424,0.585843},{0.003758,0.481113},{0.012146,0.539855},{0.006745,0.470209},{0.004193,0.441336},{0.004241,0.597987},{0.004239,0.447619},{0.004298,0.452553},{0.003841,0.50051},{0.003694,0.466988},{0.003739,0.454393},{0.003692,0.46096},{0.003239,0.77577},{0.00803,0.677106},{0.008887,0.550943},{0.012085,0.629503},{0.011328,0.700753},{0.004624,0.037448},{0.003418,0.030416},{0.003388,0.035711},{0.003187,0.028575},{0.003415,0.026698},{0.003462,0.026366},{0.00313,0.025726}}

 ListPlot[{Transpose[r][[2]],10 Transpose[r][[1]]}]

enter image description here

From the figure one can see that whereas the time required for the matrix multiplication remains almost constant for the matrix power larger than 10, the time for computing the matrix rank increases linearly with the matrix power till the power of about the half of the matrix dimension and then abruptly drops down. What is the reason for this behavior and how can the initial linear behavior be overcome?

Improve this question
$\endgroup$
8
  • $\begingroup$ One of the problems here is that the successive multiplication will increase the entries of m to an extend that they cannot be handled by machine integers (64 bit signed integers). That is, at some point. parts of the arithmetic (if not all of it) has to be emulated in software. $\endgroup$
    – Henrik Schumacher
    Mar 10 '20 at 9:02
  • $\begingroup$ You can try to do calculations with packed arrays of real number, e.g., with setting M = Developer`ToPackedArray@N@M. However, this will inevitiably lead to rounding errors that might be substantial for your application. Successive powers of M might have higher and higher condition number.At latest once you are above, say 1/ $MachineEpsilon, total chaos breaks loose. $\endgroup$
    – Henrik Schumacher
    Mar 10 '20 at 9:05
  • 1
    $\begingroup$ Oh, in that case, a JordanDecomposition might be more efficient. Sorry, I have not time at the moment to fill in the details.... $\endgroup$
    – Henrik Schumacher
    Mar 10 '20 at 14:08
  • 4
    $\begingroup$ (1) Depending on matrix size you might try numeric methods such as the singular values decomposition (best to do this using precision greater than machine precision). (2) No need to iterate beyond the number of zeros on the main diagonal. $\endgroup$
    – Daniel Lichtblau
    Mar 10 '20 at 14:52
  • 1
    $\begingroup$ Another thing worth checking is the factorization of the characteristic polynomial. If it is irreducible (over the rationals) then each multiplication will drop the rank by 1. $\endgroup$
    – Daniel Lichtblau
    Mar 11 '20 at 18:25
4
$\begingroup$

Mathematica's function MatrixRank can take a Modulus argument when used with integer matrices. In my experiments (using a matrix m computed using your code) I noted:

  1. It is very much faster than an exact computation.
  2. It gives the correct answer most of the time when using a prime modulus that is not too small.

For example:

{#, MatrixRank[m, Modulus -> #]} & /@ Array[Prime, 100]
(* {{2, 45}, {3, 57}, {5, 81}, {7, 82}, {11, 86}, {13, 86}, {17,
   86}, {19, 86}, {23, 86}, {29, 86}, {31, 86}, {37, 86}, {41, 
  86}, {43, 86}, {47, 86}, {53, 86}, {59, 86}, {61, 86}, {67, 
  86}, {71, 86}, {73, 86}, {79, 86}, {83, 86}, {89, 86}, {97, 
  86}, {101, 86}, {103, 86}, {107, 86}, {109, 86}, {113, 86}, {127, 
  86}, {131, 86}, {137, 86}, {139, 86}, {149, 86}, {151, 86}, {157, 
  86}, {163, 86}, {167, 86}, {173, 86}, {179, 86}, {181, 86}, {191, 
  86}, {193, 86}, {197, 86}, {199, 86}, {211, 86}, {223, 86}, {227, 
  86}, {229, 86}, {233, 86}, {239, 86}, {241, 86}, {251, 86}, {257, 
  86}, {263, 86}, {269, 86}, {271, 86}, {277, 86}, {281, 86}, {283, 
  86}, {293, 86}, {307, 86}, {311, 86}, {313, 86}, {317, 86}, {331, 
  86}, {337, 86}, {347, 86}, {349, 86}, {353, 86}, {359, 86}, {367, 
  86}, {373, 86}, {379, 86}, {383, 86}, {389, 86}, {397, 86}, {401, 
  86}, {409, 86}, {419, 86}, {421, 86}, {431, 86}, {433, 86}, {439, 
  86}, {443, 86}, {449, 86}, {457, 86}, {461, 86}, {463, 86}, {467, 
  86}, {479, 86}, {487, 86}, {491, 86}, {499, 86}, {503, 86}, {509, 
  86}, {521, 86}, {523, 86}, {541, 86}} *)

I don't know whether there are any proven results on this or whether an approximate answer is good enough for you.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ This is a great idea! However I need the exact result. Therefore it would be nice to understand what value of the prime number guarantees the correct result. $\endgroup$
    – drer
    Mar 11 '20 at 13:55
  • 1
    $\begingroup$ It is effectively impossible to precompute a prime that is guaranteed to work. There are methods to estimate/bound probability of failure though. $\endgroup$
    – Daniel Lichtblau
    Mar 11 '20 at 18:22

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.