Exponentiation a matrix by itself N times?

Question

I am realizing Exponentiation of a matrix using FOR:

import numpy as np
fl=2
cl=2
fl2=fl
cl2=cl
M = random.random((fl,cl))
M2 = M
Result = np.zeros((fl,cl))
Temp = np.zeros((fl,cl))
itera = 2
print('Matriz A:\n',M)
print('Matriz AxA:\n',M2)

for i in range (0,itera):
    for a in range(0,fl):
        for b in range (0,cl):  
            Result[a,b]+=M[a,b]*M[a,b]
            temp[a,b]=Result[a,b]
            Res[a,k]=M[a,b]
print('Potencia:\n',temp)
print('Matriz:\n', Result)

The error is that it does not perform well the multiplication in Result[a,b]+=M[a,b]*M[a,b] and when I save it in a temporary matrix to multiply it with the original matrix, it does not make the next jump in for i in range (0,itera):

I know I can perform the function np.matmul but I try to do it with the FOR loop

Example

Answer 1

There are some problems here:

1. you do not reset the result matrix after multiplication is done, hence you keep adding more values; and
2. you never assign the result back to m to perform a next generation of multiplications.

Naive power implementation

I think it is also better to "encapsulate" matrix multiplication in a separate function, like:

def matmul(a1, a2):
    m, ka = a1.shape
    kb, n = a2.shape
    if ka != kb:
        raise ValueError()
    res = np.zeros((m, n))
    for i in range(m):
        for j in range(n):
            d = 0.0
            for k in range(ka):
                d += a1[i,k] * a2[k,j]
            res[i, j] = d
    return res

Then we can calculate the power of this matrix with:

m2 = m
for i in range(topow-1):
    m = matmul(m, m2)

Note that we can not use m here as the only matrix. Since if we write m = matmul(m, m), then m is now m2. But that means that if we perform the multiplication a second time, we get m4 instead of m3.

This then produces the expected results:

>>> cross = np.array([[1,0,1],[0,1,0], [1,0,1]])
>>> matmul(cross, cross)
array([[2., 0., 2.],
       [0., 1., 0.],
       [2., 0., 2.]])
>>> matmul(cross, matmul(cross, cross))
array([[4., 0., 4.],
       [0., 1., 0.],
       [4., 0., 4.]])
>>> matmul(cross, matmul(cross, matmul(cross, cross)))
array([[8., 0., 8.],
       [0., 1., 0.],
       [8., 0., 8.]])

Logarithmic power multiplication

The above can calculate the Mⁿ in O(n) (linear time), but we can do better, we can calculate this matrix in logarithmic time: we do this by looking if the power is 1, if it is, we simply return the matrix, if it is not, we check if the power is even, if it is even, we multiply the matrix with itself, and calculate the power of that matrix, but with the power divided by two, so M^{2 n}=(M×M)ⁿ. If the power is odd, we do more or less the same, except that we multiply it with the original value for M: M^{2 n + 1}=M×(M×M)ⁿ. Like:

def matpow(m, p):
    if p <= 0:
        raise ValueError()
    if p == 1:
        return m
    elif p % 2 == 0:  # even
        return matpow(matmul(m, m), p // 2)
    else:             # odd
        return matmul(m, matpow(matmul(m, m), p // 2))

The above can be written more elegantly, but I leave this as an exercise :).

Note however that using numpy arrays for scalar comuputations is typically less efficient than using the matrix multiplication (and other functions) numpy offers. These are optimized, and are not interpreted, and typically outperform Python equivalents significantly. Therefore I would really advice you to use these. The numpy functions are also tested, making it less likely that there are bugs in it.

Answer 2

You're looking for np.linalg.matrix_power.

If you're using numpy, don't use a for loop, use a vectorized operation.

arr = np.arange(16).reshape((4,4))

np.linalg.matrix_power(arr, 3)

array([[ 1680,  1940,  2200,  2460],
       [ 4880,  5620,  6360,  7100],
       [ 8080,  9300, 10520, 11740],
       [11280, 12980, 14680, 16380]])

Which is the same as the explicit multiplication:

arr @ arr @ arr

---

>>> np.array_equal(arr @ arr @ arr, np.linalg.matrix_power(arr, 3))
True

---

Since you asked

If you really want a naive solution using loops, we can put together the pieces quite easily. First we need a way to actually multiple the matrices. There are options that beat n^3 complexity, this answer is not going to do that. Here is a basic matrix multiplication function:

def matmultiply(a, b):
  res = np.zeros(a.shape)
  size = a.shape[0]

  for i in range(size):
    for j in range(size):
      for k in range(size):
        res[i][j] += a[i][k] * b[k][j]

  return res

Now you need an exponential function. This function takes a matrix and a power, and raises a matrix to that power.

def loopy_matrix_power(a, n):
  res = np.identity(a.shape[0])
  while n > 0:
    if n % 2 == 0:
      a = matmultiply(a, a)
      n /= 2
    else:
      res = matmultiply(res, a)
      n -= 1

  return res

In action:

loopy_matrix_power(arr, 3)

array([[ 1680.,  1940.,  2200.,  2460.],
       [ 4880.,  5620.,  6360.,  7100.],
       [ 8080.,  9300., 10520., 11740.],
       [11280., 12980., 14680., 16380.]])