CODEWITHSUNDEEP
pythonnumpyrecursionadjacency-matrixarray-broadcasting
mo-alowais
mo-alowais

Reputation: 23

Recursive matrix construction with numpy array issues (broadcasting?)

The following is a seemingly simple recursion for finding a hybercube matrix.

The recursion is defined as:

(Formula)

I tried to put it into code but I keep running into broadcasting issues with numpy.

I tried instead of taking the size of the previous matrix, just taking powers of 2 to reduce recursive calls (for the parameter of the identity matrix).

import numpy as np
from numpy import linalg as LA

Q1 = np.array([[0, 1],
               [1, 0]])

def Qn(n):
    if n <= 1:
        return Q1
    else:
        return np.array([[Qn(n-1), np.identity(int(np.exp2(n-1)))],
                         [np.identity(int(np.exp2(n-1))), Qn(n-1)]])
    
Q3= Qn(3)

eig_value, eig_vectors = LA.eig(Q3)
print(eig_value)

Q1 is the base case of my matrix. It should be very simple, yet I keep running into issues.

Traceback (most recent call last):
File "e:\Coding\Python\test.py", line 15, in <module>
Q3= Qn(3)
File "e:\Coding\Python\test.py", line 12, in Qn
return np.array([[Qn(n-1), np.identity(int(np.exp2(n-1)))],
ValueError: setting an array element with a sequence. The 
requested array has an inhomogeneous shape after 2 dimensions. The 
detected shape was (2, 2) + inhomogeneous part.

I get this error ^^

Upvotes: 2

Views: 64

Answers (2)

Andrej Kesely
Andrej Kesely

Reputation: 195573

Another option, using np.r_/np.c_:

import numpy as np
from numpy import linalg as LA

Q1 = np.array([[0, 1], [1, 0]])


def Qn(n):
    if n <= 1:
        return Q1
    else:
        I = np.eye(2 ** (n - 1))
        Qnm1 = Qn(n - 1)
        return np.r_[np.c_[Qnm1, I], np.c_[I, Qnm1]]


Q3 = Qn(3)

eig_value, eig_vectors = LA.eig(Q3)
print(eig_value)

Prints:

[ 3. -3. -1. -1. -1.  1.  1.  1.]

Upvotes: 1

Matt Haberland
Matt Haberland

Reputation: 3873

Use np.block instead of np.array to assemble the block matrix.

import numpy as np
from numpy import linalg as LA

Q1 = np.array([[0, 1],
               [1, 0]])

def Qn(n):
    if n <= 1:
        return Q1
    else:
        Qnm1 = Qn(n-1)
        I = np.eye(2**(n-1))
        return np.block([[Qnm1, I], [I, Qnm1]])
    
Q3 = Qn(3)

eig_value, eig_vectors = LA.eig(Q3)
print(eig_value)
# [ 3. -3. -1. -1. -1.  1.  1.  1.]

Upvotes: 2

Related Questions