Change numpy array inside recursive function-CodePudding

Question 4.1-1 from CLRS Introduction to Algorithms

I have code that multiply n*n matrix, where n is exact power of 2, using recursion. I need to modify this code, so that n doesn't need to be exact power of 2. Obvious way to do it is adding zeroes to matrix so it became k*k, where k is power of 2.
I try two ways to modify matrix, but it seems that when I try to add zeroes to C, it's only add it to local version of C inside function, but global version remains the same. I don't see ways to using return statement since I have recursion.
Code works when n is power of 2 but raise error otherwise.

def matrix_multiply_recursive_aux(A, B, C, n, A_row, A_col, B_row, B_col, C_row, C_col):
    """Compute C = C   (A * B), where A, B, and C are n x n matrices.
       Other parameters:
        A_row, A_col: starting row and column for the n x n submatrix of A
        B_row, B_col: starting row and column for the n x n submatrix of B
        C_row, C_col: starting row and column for the n x n submatrix of C
    """
    if n == 1:
        C[C_row, C_col]  = A[A_row, A_col] * B[B_row, B_col]  # base case
        return

# Code that I add, first version
    if np.log2(n) % 1 != 0:
        k = int(2 ** np.ceil(np.log2(n)))
        for i in range(n, k):
            C = np.append(C, [[0] * n], axis=0)
        for i in range(n, k):
            C = np.column_stack((C, [0] * k))
        n = k
# end of my first version

# second version

#     if np.log2(n) % 1 != 0:
#         C_temp = np.zeros((k, k))
#         for i in range(n):
#             for j in range(n):
#                 C_temp[i, j] = C[i, j]
#         C = C_temp
#         n = int(2 ** np.ceil(np.log2(n)))

# end of second version

    half = n // 2
    matrix_multiply_recursive_aux(A, B, C, half, A_row, A_col, B_row, B_col, C_row, C_col)
    matrix_multiply_recursive_aux(A, B, C, half, A_row, A_col, B_row, B_col   half, C_row, C_col   half)
    matrix_multiply_recursive_aux(A, B, C, half, A_row   half, A_col, B_row, B_col, C_row   half, C_col)
    matrix_multiply_recursive_aux(A, B, C, half, A_row   half, A_col, B_row, B_col   half, C_row   half, C_col   half)
    matrix_multiply_recursive_aux(A, B, C, half, A_row, A_col   half, B_row   half, B_col, C_row, C_col)
    matrix_multiply_recursive_aux(A, B, C, half, A_row, A_col   half, B_row   half, B_col   half, C_row, C_col   half)
    matrix_multiply_recursive_aux(A, B, C, half, A_row   half, A_col   half, B_row   half, B_col, C_row   half, C_col)
    matrix_multiply_recursive_aux(A, B, C, half, A_row   half, A_col   half, B_row   half, B_col   half, C_row   half, C_col   half)

When I run code I get error.

A = np.array([[1, 2, 0],
              [0, 3, 1],
              [1, 0, 2]])
B = np.array([[1, 3, 1],
              [1, 0, 0],
              [1, 1, 2]])
n = 3
C = np.array([[0, 0, 0], 
              [0, 0, 0], 
              [0, 0, 0]])
matrix_multiply_recursive_aux(A, B, C, n, 0, 0, 0, 0, 0, 0)

      7     """
      8     if n == 1:
----> 9         C[C_row, C_col]  = A[A_row, A_col] * B[B_row, B_col]  # base case
     10         return
     11 

IndexError: index 3 is out of bounds for axis 1 with size 3

When I comment recursive calls (last eight lines) and try to print only C I get it 3*3 when I expect 4*4.

matrix_multiply_recursive_aux(A, B, C, n, 0, 0, 0, 0, 0, 0)
print(C)

[[0 0 0]
 [0 0 0]
 [0 0 0]]

How can I change nonlocal version of array?

CodePudding user response：

I tested it only with your example - so I don't know if there are other problems.

I see two problems:

You resize C but you have to resize also A and B to make calculations on matrixes with the same shapes. And this was generating IndexError: index 3 is out of bounds for axis 1 with size 3 because you wanted to create C w shape 4x4 using A and B with shapes 3x3 and it couldn't find last row/column in A and B

if np.log2(n) % 1 != 0:
    k = int(2 ** np.ceil(np.log2(n)))
    for i in range(n, k):
        C = np.append(C, [[0] * n], axis=0)
        A = np.append(A, [[0] * n], axis=0)
        B = np.append(B, [[0] * n], axis=0)
    for i in range(n, k):
        C = np.column_stack((C, [0] * k))
        A = np.column_stack((A, [0] * k))
        B = np.column_stack((B, [0] * k))
    n = k

Python sends reference to C (reference to the some place in memory) (without dupicate data) and in original version all recursions work with the same data in memory - but when you use np.append() and C = ... then it creates new matrix in memory and it assign to C reference to new place in memory - and recursions work with different values in C.

One solution is to use global C and remove C from arguments.

Other solution is to use return C and C = matrix_multiply_recursive_aux() to send reference to other recursions.

Version with return seems cleaner and better for debuging.

Version with return

In code I calculate A @ B to compare result.

def matrix_multiply_recursive_aux(A, B, C, n, A_row, A_col, B_row, B_col, C_row, C_col):
    """Compute C = C   (A * B), where A, B, and C are n x n matrices.
       Other parameters:
        A_row, A_col: starting row and column for the n x n submatrix of A
        B_row, B_col: starting row and column for the n x n submatrix of B
        C_row, C_col: starting row and column for the n x n submatrix of C
    """
    
    if n == 1:
        C[C_row, C_col]  = A[A_row, A_col] * B[B_row, B_col]  # base case
        return C


    if np.log2(n) % 1 != 0:
        k = int(2 ** np.ceil(np.log2(n)))
        for i in range(n, k):
            C = np.append(C, [[0] * n], axis=0)
            A = np.append(A, [[0] * n], axis=0)
            B = np.append(B, [[0] * n], axis=0)
        for i in range(n, k):
            C = np.column_stack((C, [0] * k))
            A = np.column_stack((A, [0] * k))
            B = np.column_stack((B, [0] * k))
        n = k

    half = n // 2
    C = matrix_multiply_recursive_aux(A, B, C, half, A_row, A_col, B_row, B_col, C_row, C_col)
    C = matrix_multiply_recursive_aux(A, B, C, half, A_row, A_col, B_row, B_col   half, C_row, C_col   half)
    C = matrix_multiply_recursive_aux(A, B, C, half, A_row   half, A_col, B_row, B_col, C_row   half, C_col)
    C = matrix_multiply_recursive_aux(A, B, C, half, A_row   half, A_col, B_row, B_col   half, C_row   half, C_col   half)
    C = matrix_multiply_recursive_aux(A, B, C, half, A_row, A_col   half, B_row   half, B_col, C_row, C_col)
    C = matrix_multiply_recursive_aux(A, B, C, half, A_row, A_col   half, B_row   half, B_col   half, C_row, C_col   half)
    C = matrix_multiply_recursive_aux(A, B, C, half, A_row   half, A_col   half, B_row   half, B_col, C_row   half, C_col)
    C = matrix_multiply_recursive_aux(A, B, C, half, A_row   half, A_col   half, B_row   half, B_col   half, C_row   half, C_col   half)
    
    return C

# ------------------------------------------

import numpy as np

A = np.array([[1, 2, 0],
              [0, 3, 1],
              [1, 0, 2]])
B = np.array([[1, 3, 1],
              [1, 0, 0],
              [1, 1, 2]])

print('--- A @ B ---')
print(A @ B) 

n = 3
C = np.array([[0, 0, 0], 
              [0, 0, 0], 
              [0, 0, 0]])

C = matrix_multiply_recursive_aux(A, B, C, n, 0, 0, 0, 0, 0, 0)

print('--- matrix_multiply_recursive_aux ---')
print(C)

Result:

--- A @ B ---
[[3 3 1]
 [4 1 2]
 [3 5 5]]
--- matrix_multiply_recursive_aux ---
[[3 3 1 0]
 [4 1 2 0]
 [3 5 5 0]
 [0 0 0 0]]

Result needs to remove added rows and colums - so other idea is to create function which adds all needed rows and columns, next runs original function with recursion, and next remove rows and columns from result.