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

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
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;

• 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. May 28 '21 at 11:12
• Yes, $M$ is a matrix and $t$ is an integer. May 28 '21 at 12:23
• Ok. And $I$ would correspond to an identity matrix? And $*$ is a multiplication sign..correct? Sorry for this naive questions. May 28 '21 at 13:36
• Correct.$\phantom{}$ May 28 '21 at 13:55
• 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. May 28 '21 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.