Умножение матриц с использованием разделяй и властвуй в Python

Question

Я написал программу на python3, чтобы найти произведение 2 n*n матриц (, где n - степень 2 ).

Почему приведенный ниже код не работает и показывает IndexError: invalid index to scalar variable?

import numpy as np

def product(x, y, k):
    def fsum(p, q, m):
        r = [[p[i, j] + q[i, j] for j in range(m)] for i in range(m)]
        return r

    if k == 1:
        return x[0][0] * y[0][0]
    else:
        A = x[0:(k // 2), 0:(k // 2)]
        B = x[0:(k // 2), (k // 2):k]
        C = x[(k // 2):k, 0:(k // 2)]
        D = x[(k // 2):k, (k // 2):k]

        E = y[0:(k // 2), 0:(k // 2)]
        F = y[0:(k // 2), (k // 2):k]
        G = y[(k // 2):k, 0:(k // 2)]
        H = y[(k // 2):k, (k // 2):k]

        C00 = fsum(product(A, E, k // 2), product(B, G, k // 2), k // 2)
        C01 = fsum(product(A, F, k // 2), product(B, H, k // 2), k // 2)
        C10 = fsum(product(C, E, k // 2), product(D, G, k // 2), k // 2)
        C11 = fsum(product(C, F, k // 2), product(D, H, k // 2), k // 2)

        return np.array([[C00, C01], [C10, C11]])

n = int(input('Enter index(power of 2): '))
print('Input 1st matrix')
a = np.array([[int(_) for _ in input().split()] for x in range(n)])
print('Input 2nd matrix')
b = np.array([[int(_) for _ in input().split()] for x in range(n)])
print(product(a, b, n))

Answer 1

Пример прогона:

Enter index(power of 2): 2
Input 1st matrix
1 2
3 4
[[1 2]
 [3 4]]
Input 2nd matrix
3 4
5 6
[[3 4]
 [5 6]]
Traceback (most recent call last):
  File "stack51061196.py", line 35, in <module>
    print(product(a, b, n))
  File "stack51061196.py", line 21, in product
    C00 = fsum(product(A, E, k // 2), product(B, G, k // 2), k // 2)
  File "stack51061196.py", line 5, in fsum
    r = [[p[i, j] + q[i, j] for j in range(m)] for i in range(m)]
  File "stack51061196.py", line 5, in <listcomp>
    r = [[p[i, j] + q[i, j] for j in range(m)] for i in range(m)]
  File "stack51061196.py", line 5, in <listcomp>
    r = [[p[i, j] + q[i, j] for j in range(m)] for i in range(m)]
IndexError: invalid index to scalar variable.

Добавление print('p,q', p, q, type(p), type(q)) к fsum Я вижу

p,q 1 6 <class 'numpy.int64'> <class 'numpy.int64'>

То есть p - это np.int64 объект, а не массив. Он уже проиндексирован и не может идти дальше.

In [193]: x = np.array([1])[0]
In [194]: x
Out[194]: 1
In [195]: type(x)
Out[195]: numpy.int64
In [196]: x[0,0]
IndexError: invalid index to scalar variable.

Другой диагностический отпечаток

print('product AE',k,type(product(A, E, k // 2)))

показывает

product AE 2 <class 'numpy.int64'>

Таким образом, когда k равен 2, product возвращает scalar variable, int64 объект, а не массив. Передача этого в fsum приводит к ошибке.

Именно эта ветвь product вызывает проблему - x[0][0] (почему бы не x[0,0]?) Является элементом x:

if k == 1:
    return x[0][0] * y[0][0]

На данный момент x имеет (1,1) форму, так что вы можете просто написать

if k == 1:
     return x * y

С этим изменением я получаю:

0840:~/mypy$ python3 stack51061196.py 
a,b [[1 2]
 [3 4]] [[1 2]
 [3 4]]
A [[1]] <class 'numpy.ndarray'>
product AE 2 <class 'numpy.ndarray'>
p,q [[1]] [[6]] <class 'numpy.ndarray'> <class 'numpy.ndarray'>
p,q [[2]] [[8]] <class 'numpy.ndarray'> <class 'numpy.ndarray'>
p,q [[3]] [[12]] <class 'numpy.ndarray'> <class 'numpy.ndarray'>
p,q [[6]] [[16]] <class 'numpy.ndarray'> <class 'numpy.ndarray'>
[[[[ 7]]

  [[10]]]


 [[[15]]

  [[22]]]]

, который помимо размеров соответствует:

In [197]: a=np.array([[1,2],[3,4]])
In [198]: np.dot(a,a)
Out[198]: 
array([[ 7, 10],
       [15, 22]])

Ваша форма (2, 2, 1, 1), которую можно удалить с помощью squeeze, но вам действительно нужно уточнить итерацию, чтобы получить нужную форму без нее.

Он также получает правильные числа для массива 4x4, но форма теперь (2, 2, 2, 2, 1, 1).

Answer 2

Вот один из способов сделать это:

def product(X, Y):
    k = len(X)
    if k == 1:
        return X * Y
    (A, B), (C, D) = skimage.util.view_as_blocks(X, block_shape=(k // 2, k // 2))
    (E, F), (G, H) = skimage.util.view_as_blocks(Y, block_shape=(k // 2, k // 2))

    out = np.zeros((k, k))
    (I, J), (K, L) = skimage.util.view_as_blocks(out, block_shape=(k // 2, k // 2))
    I[:] = product(A, E) + product(B, G)
    J[:] = product(A, F) + product(B, H)
    K[:] = product(C, E) + product(D, G)
    L[:] = product(C, F) + product(D, H)
    return out

Излишне говорить, что это ужасно медленно

A = numpy.random.rand(32, 32)

np.allclose(product(A, A), A @ A)
# True
%timeit product(A, A)
# 526 ms ± 11.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit A @ A
# 4.36 µs ± 59.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)