Fast computing of a matrix power for large integer values in C++

Question

I'm working with squared matrices that can be quite large, for instance, a M = 50 x 50 matrix.

My objective is to compute the power of the squared matrix M^t for very large t values (for example t = 4000).

I work in R and I have used the R function matrix.power from the matrixcalc R package.

I'm exploring the possibility to write a code for matrix power in C++ and then import it in R through the R package Rcpp.

One alternative would be to use the matrix multiplication approach (as in the matrix.power function in R), but looking around I have understood that there might be faster approaches to calculate a matrix power.

Do you have any experience on that? Does anyone know a library in C++ that does matrix power calculation fast and efficiently?

Consider that I'm working on a Mac0s laptop with 16 GB of RAM and 4 CPUs.

Answer 1

Since your matrices are small $(50 \times 50)$, you can probably just compute $M^t$ through repeated exponentiation where the exponents are powers of $2$.

Write $t$ in binary so that $t = 2^{k_1} + 2^{k_2} + \dots + 2^{k_\ell}$. Then $M^t = \prod_{i=1}^\ell M^{2^{k_i}}$. Moreover, for $k_i \ge 1$ you have $M^{2^{k_i}} = \left( M^{2^{k_i - 1}} \right)^2$, so you need at most $O(\log t)$ matrix multiplications.

Here is a pseudocode where "&" denotes "bitwise and" and "~" denotes "bitwise not":

Power(M, t):
  if(t & 1):        //Handle odd values of t (this saves a multiplication later)
     R = M;
     t = t & ~1;    //Clear the least significant bit of t
  else:
     R = I;

  i=1;
  B=M;                //B will always be M^i, where i is a power of 2
  while t!=0:
     i = i*2;         //Advance i to the next power of 2
     B = B*B;         //B was M^(i/2) and is now M^i

     if(t & i):       //i is of the form 2^j. Is the j-th bit of t set?
        R = R*B;      //Multiply the result with B=A^i
        t = t & ~i;   //Clear the j-th bit of t
  
  Return R;

Answer 2

I suggest using xtensor. You can compute the 4000-th matrix power of M as xt::linalg::matrix_power(M, 4000).

Obviously you should be aware that powering in any language can incur in numerical issues. Even if your matrix is 1 x 1, M^4000 could be enormous, larger than what you could store as a floating point value.