How to raise a numpy array to a power? (corresponding to repeated matrix multiplications, not elementwise)

I believe you want numpy.linalg.matrix_power

As a quick example:

import numpy as np
x = np.arange(9).reshape(3,3)
y = np.matrix(x)

a = y**3
b = np.linalg.matrix_power(x, 3)

print a
print b
assert np.all(a==b)

This yields:

In [19]: a
Out[19]: 
matrix([[ 180,  234,  288],
        [ 558,  720,  882],
        [ 936, 1206, 1476]])

In [20]: b
Out[20]: 
array([[ 180,  234,  288],
       [ 558,  720,  882],
       [ 936, 1206, 1476]])