title: Faster matrix power than numpy?
body: <p>I needed to compute Q^N for quite a lot of various values of N <em>(between 1 to 10000)</em> and Numpy was a bit too slow.</p>
  
  <p>I've asked on <a href="https://math.stackexchange.com/questions/499544/expressing-the-values-of-a-matrix-at-pow-n">math.stackexchange.com</a> if I could avoid to compute Q^N for my specific need and someone answered me that computing Q^N should be quite fast using the <code>P D^N P^-1</code> method.</p>
  
  <p>So basically, instead of doing:</p>
  
  <pre><code>import numpy as np
  from numpy import linalg as LA
  ...
  LA.matrix_power(m, N)
  </code></pre>
  
  <p>I've tried :</p>
  
  <pre><code>diag, P = LA.eig(m)
  DN = np.diag(diag**N)
  P1 = LA.inv(P)
  
  P*DN*P1
  </code></pre>
  
  <p>And I obtain the same matrix as result <em>(tried on a single example)</em></p>
  
  <p>On a more complex matrix, Q:</p>
  
  <pre><code>% timeit.Timer('Q**10000', setup=setup).repeat(2, 100)
  [5.87254786491394, 5.863131046295166]
  
  % timeit.Timer('diag, P = linalg.eig(Q); DN=np.diag(diag**10000);P1=linalg.inv(P); P*DN*P1', setup=setup).repeat(2, 100)
  [2.0032401084899902, 2.018735885620117]
  </code></pre>
  
  <p>And regarding my initial problem, the second method allows me to compute <code>P, diag and P1</code>only once and use it thousands of times. It's 8 times faster using this method.</p>
  
  <p><strong>My questions are:</strong></p>
  
  <ul>
  <li>In which case it is not possible to use this last method to compute Q^N?</li>
  <li>Is it fine to use it in my case (matrix Q as given <a href="https://math.stackexchange.com/questions/499544/expressing-the-values-of-a-matrix-at-pow-n">here</a>)?</li>
  <li>Is there in numpy a function that already does it?</li>
  </ul>
  
  <p><strong>Edit:</strong></p>
  
  <ul>
  <li>It appears that for another matrix, P is not invertible. So I am adding a new question: how can I change the matrix P so it becomes invertible but the resulting matrix is not too altered? I mean, it's ok if the values are close to the real result, by close I mean ~0.0001.</li>
  </ul>