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]]}]
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?
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$M = Developer`ToPackedArray@N@M
. However, this will inevitiably lead to rounding errors that might be substantial for your application. Successive powers ofM
might have higher and higher condition number.At latest once you are above, say1/ $MachineEpsilon
, total chaos breaks loose. $\endgroup$JordanDecomposition
might be more efficient. Sorry, I have not time at the moment to fill in the details.... $\endgroup$