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.


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
     R = I;

  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;
  • $\begingroup$ Sorry for my naive question, but I'm trying to implement your pseudocode in C++. M would be a matrix, whereas t is an integer? I am asking because I need to define the object types in the Power function definition. For instance, in the case of t=4000 the equivalent binary number would be: 111110100000. $\endgroup$ May 28 at 11:12
  • $\begingroup$ Yes, $M$ is a matrix and $t$ is an integer. $\endgroup$
    – Steven
    May 28 at 12:23
  • $\begingroup$ Ok. And $I$ would correspond to an identity matrix? And $*$ is a multiplication sign..correct? Sorry for this naive questions. $\endgroup$ May 28 at 13:36
  • $\begingroup$ Correct.$\phantom{}$ $\endgroup$
    – Steven
    May 28 at 13:55
  • $\begingroup$ I tried implementing your pseudocode in C++. I compared the results of the C++ code with the results of another pre-made functions in R. They are different, which means that I made some mistake. Is there a way to show you my code? It is quite short I would say. $\endgroup$ May 28 at 15:30

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.


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.