Bdtjtk

3D numpy array A contains a series (in this example, I am choosing 3) of 2D numpy array D of shape 2 x 2. The D matrix is as follows:

D = np.array([[1,2],[3,4]])

A is initialized and assigned as below:

idx = np.arange(3)

A = np.zeros((3,2,2))

A[idx,:,:] = D           # This gives A = [[[1,2],[3,4]],[[1,2],[3,4]],

                         # [[1,2],[3,4]]]

                         # In mathematical notation: A = {D, D, D}

Now, essentially what I require after the execution of the codes is:

Mathematically, A = {D^0, D^1, D^2} = {D0, D1, D2}
where D0 = [[1,0],[0,1]], D1 = [[1,2],[3,4]], D2=[[7,10],[15,22]]

Is it possible to apply power to each matrix element in A without using a for-loop? I would be doing larger matrices with more in the series.

I had defined, n = np.array([0,1,2]) # corresponding to powers 0, 1 and 2 and tried

Result = np.power(A,n) but I do not get the desired output.

Is there are an efficient way to do it?

Full code:

D = np.array([[1,2],[3,4]])

idx = np.arange(3)

A = np.zeros((3,2,2))

A[idx,:,:] = D           # This gives A = [[[1,2],[3,4]],[[1,2],[3,4]],

                         # [[1,2],[3,4]]]

                         # In mathematical notation: A = {D, D, D}

n = np.array([0,1,2])

Result = np.power(A,n)   # ------> Not the desired output.

A cumulative product exists in numpy, but not for matrices. Therefore, you need to make your own 'matcumprod' function. You can use np.dot for this, but np.matmul (or @) is specialized for matrix multiplication.

Since you state your powers always go from 0 to some_power, I suggest the following function:

def matcumprod(D, upto):

  Res = np.empty((upto, *D.shape), dtype=A.dtype)

  Res[0, :, :] = np.eye(D.shape[0])

  Res[1, :, :] = D.copy()

  for i in range(1,upto):

    Res[i, :, :] = Res[i-1,:,:] @ D



  return Res

By the way, a loop often times outperforms a built-in numpy function if the latter uses a lot of memory, so don't fret over it if your powers stay within bounds...

Alright, i spent a lot of time on this problem but could not seem to find a vectorized solution in the way you'd like. So i would like to instead first propose a basic solution, and then perhaps an optimization if you require finding continuous powers.

The function you're looking for is called numpy.linalg.matrix_power

import numpy as np



D = np.matrix([[1,2],[3,4]])

idx = np.arange(3)

A = np.zeros((3,2,2))

A[idx,:,:] = D           # This gives A = [[[1,2],[3,4]],[[1,2],[3,4]],

                         # [[1,2],[3,4]]]

                         # In mathematical notation: A = {D, D, D}

np.zeros(A.shape)

n = np.array([0,1,2])

result = [np.linalg.matrix_power(D, i) for i in n]

np.array(result)

#Output:

array([[[ 1,  0],

        [ 0,  1]],



       [[ 1,  2],

        [ 3,  4]],



       [[ 7, 10],

        [15, 22]]])

However, if you notice, you end up calculating multiple powers for the same base matrix. We could instead utilize the intermediate results and go from there, using numpy.linalg.multi_dot

def all_powers_arr_of_matrix(A): 

    result = np.zeros(A.shape)

    result[0] = np.linalg.matrix_power(A[0], 0)

    for i in range(1, A.shape[0]):

        result[i] = np.linalg.multi_dot([result[i - 1], A[i]])

    return result

    result = all_powers_arr_of_matrix(A)

#Output:

array([[[ 1.,  0.],

        [ 0.,  1.]],



       [[ 1.,  2.],

        [ 3.,  4.]],



       [[ 7., 10.],

        [15., 22.]]])

Also, we can avoid creating the matrix A entirely, saving some time.

    def all_powers_matrix(D, *rangeargs): #end exclusive

        ''' Expects 2D matrix. 

        Use as all_powers_matrix(D, end) or

        all_powers_matrix(D, start, end)

        '''

        if len(rangeargs) == 1:

            start = 0

            end = rangeargs[0]

        elif len(rangeargs) == 2:

            start = rangeargs[0]

            end = rangeargs[1]

        else:

            print("incorrect args")

            return None

        result = np.zeros((end - start, *D.shape))

        result[0] = np.linalg.matrix_power(A[0], start)

        for i in range(start + 1, end):

            result[i] = np.linalg.multi_dot([result[i - 1], D])

        return result



            return result

result = all_powers_matrix(D, 3)

#Output:

array([[[ 1.,  0.],

        [ 0.,  1.]],



       [[ 1.,  2.],

        [ 3.,  4.]],



       [[ 7., 10.],

        [15., 22.]]])

Note that you'd need to add error handling if you decide to use these functions as-is.

I don't have a full solution, but there are some things I wanted to mention which are a bit too long for the comments.

You might first look into addition chain exponentiation if you are computing big powers of big matrices. This is basically asking how many matrix multiplications are required to compute A^k for a given k. For instance A^5 = A(A^2)^2 so you need to only three matrix multiplies: A^2 and (A^2)^2 and A(A^2)^2. This might be the simplest way to gain some efficiency, but you will probably still have to use explicit loops.

Your question is also related to the problem of computing Ax, A^2x, ... , A^kx for a given A and x. This is an active area of research right now (search "matrix powers kernel"), since computing such a sequence efficiently is useful for parallel/communication avoiding Krylov subspace methods. If you're looking for a very efficient solution to your problem it might be worth looking into some of the results about this.

To calculate power of matrix D, one way could be to find the eigenvalues and right eigenvectors of it with np.linalg.eig and then raise the power of the diagonal matrix as it is easier, then after some manipulation, you can use two np.einsum to calculate A

#get eigvalues and eigvectors

eigval, eigvect = np.linalg.eig(D)



# to check how it works, you can do:

print (np.dot(eigvect*eigval,np.linalg.inv(eigvect)))

#[[1. 2.]

# [3. 4.]]

# so you get back on D



#use power as ufunc of outer with n on the eigenvalues to get all the one you want

arrp = np.power.outer( eigval, n).T



#apply_along_axis to create the diagonal matrix along the last axis

diagp = np.apply_along_axis( np.diag, axis=-1, arr=arrp)



#finally use two np.einsum to calculate with the subscript to get what you want

A = np.einsum('lij,jk -> lik',

              np.einsum('ij,kjl -> kil',eigvect,diagp), np.linalg.inv(eigvect)).round()

print (A)

print (A.shape)



#[[[ 1.  0.]

#  [-0.  1.]]

#

# [[ 1.  2.]

#  [ 3.  4.]]

#

# [[ 7. 10.]

#  [15. 22.]]]

#

#(3, 2, 2)