На этом шаге мы рассмотрим реализацию этого метода.
Пусть A и B - матрицы размерностей p*q и q*r соответственно. Размер задачи зависит от этих трех размерностей p, q и r. Тривиальное начальное условие имеет место при p = q = r = 1, когда результат - просто число. Кроме того, некоторые реализации могут рассматривать размерности 0. В таких случаях результатом должна быть пустая матрица, о чём мы также коротко скажем.
Довольно простой способ декомпозиции задачи - разделить каждую матрицу на четыре блока (образующих блочную матрицу размерности 2х2). В этом случае их произведение определяется следующим образом:
⌈ ⌉ ⌈ ⌉ ⌈ ⌉ A1,1 A1,2 B1,1 B1,2 A1,1B1,1 + A1,2B2,1 A1,1B2,1 + A1,2B2,2 = (6.5) A2,1 A2,2 B2,1 B2,2 A2,1B1,1 + A2,2B2,1 A2,1B1,2 + A2,2B2,2 ⌊ ⌋ ⌊ ⌋ ⌊ ⌋
Обратите внимание, что формула похожа на умножение двух матриц размерности 2х2. Например, верхний левый блок результата (A1,1B1,1 + A1,2B2,1) можно считать произведением первой строки (блоков) А и первого столбца (блоков) B.
Декомпозиция приводит к вычислению восьми более простых произведений матриц. Таким образом, в рекурсивном условии метод вызовет себя восемь раз. Результат каждого произведения должен добавляться и помещаться в стек так, чтобы получить матрицу-результат. В примере 6.9 приводится возможная реализация.
Пример 6.9. Умножение матриц по принципу "разделяй и властвуй"
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
import numpy as np def matrix_mult(A, B): p = A.shape[0] q = A.shape[1] r = B.shape[1] if p == 0 or q == 0 or r == 0: return np.zeros((p, r)) elif p == 1 and q == 1 and r == 1: return np.matrix([[A[0, 0] * B[0, 0]]]) else: A11 = A[0:p // 2, 0:q // 2] A21 = A[p // 2:p, 0:q // 2] A12 = A[0:p // 2, q // 2:q] A22 = A[p // 2:p, q // 2:q] B11 = B[0:q // 2, 0:r // 2] B21 = B[q // 2:q, 0:r // 2] B12 = B[0:q // 2, r // 2:r] B22 = B[q // 2:q, r // 2:r] C11 = matrix_mult(A11, B11) + matrix_mult(A12, B21) C12 = matrix_mult(A11, B12) + matrix_mult(A12, B22) C21 = matrix_mult(A21, B11) + matrix_mult(A22, B21) C22 = matrix_mult(A21, B12) + matrix_mult(A22, B22) return np.vstack([np.hstack([C11, C12]), np.hstack([C21, C22])]) A = np.matrix([[2, 3, 1, -3], [4, -2, 1, 2]]) B = np.matrix([[2, 3, 1], [4, -1, -5], [0, -6, 3], [1, -1, 1]]) print(matrix_mult(A, B)) |
Рекурсивное условие сначала определяет каждую из меньших блочных матриц, добавляет их произведения и выстраивает матрицу-результат методами vstack и hstack. Одно из начальных условий вычисляет простое произведение при p = q = r = 1. Кроме того, в коде учитывается также случай пустых входных матриц, поскольку они могут появиться в рекурсивном условии при разбиении матрицы, одна из размерностей которой равна 1 (очевидно, вектор не может быть разбит на четыре вектора согласно (6.5)). Таким образом, если какая-то из размерностей равна 0, в соответствующем начальном условии создаётся пустая матрица размерности p*r, которая должным образом обрабатывается средствами языка Python.
Предыдущий метод создаёт матрицы размерности 1x1 (в начальном условии) и постепенно накапливает их в стеке для формирования конечной матрицы размерности p*r. Кроме того, отметим, что в методы не передаются размерности входных матриц.
Другой более эффективный способ состоит в передаче матриц А и B целиком в каждом вызове метода, тогда как сами перемножаемые блоки задаются как дополнительные параметры своими границами внутри этих матриц (см. пример 1.6). Кроме того, результат можно сохранить в матрице-параметре C размерности p*r (передаваемой по ссылке). В примере 6.10 приведена возможная реализация этого альтернативного решения.
Пример 6.10. Альтернативное умножение матриц по принципу «разделяй и властвуй»
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 |
import numpy as np def add_matrices_limits(A, B, C, lp, up, lr, ur): for i in range(lp, up + 1): for k in range(lr, ur + 1): C[i, k] = A[i, k] + B[i, k] def matrix_mult_limits(A, B, C, lp, up, lq, uq, lr, ur): mp = (lp + up) // 2 mq = (lq + uq) // 2 mr = (lr + ur) // 2 if lp == up and lq == uq and lr == ur: C[mp, mr] = A[mp, mq] * B[mq, mr] elif lp <= up and lq <= uq and lr <= ur: M1 = np.zeros((A.shape[0], B.shape[1])) M2 = np.zeros((A.shape[0], B.shape[1])) matrix_mult_limits(A, B, M1, lp, mp, lq, mq, lr, mr) matrix_mult_limits(A, B, M2, lp, mp, mq + 1, uq, lr, mr) add_matrices_limits(M1, M2, C, lp, mp, lr, mr) matrix_mult_limits(A, B, M1, lp, mp, lq, mq, mr + 1, ur) matrix_mult_limits(A, B, M2, lp, mp, mq + 1, uq, mr + 1, ur) add_matrices_limits(M1, M2, C, lp, mp, mr + 1, ur) matrix_mult_limits(A, B, M1, mp + 1, up, lq, mq, lr, mr) matrix_mult_limits(A, B, M2, mp + 1, up, mq + 1, uq, lr, mr) add_matrices_limits(M1, M2, C, mp + 1, up, lr, mr) matrix_mult_limits(A, B, M1, mp + 1, up, lq, mq, mr + 1, ur) matrix_mult_limits( A, B, M2, mp + 1, up, mq + 1, uq, mr + 1, ur) add_matrices_limits(M1, M2, C, mp + 1, up, mr + 1, ur) def matrix_mult_limits_vrapper(A, B): C = np.zeros((A.shape[0], B.shape[1])) matrix_mult_limits(A, B, C, 0, A.shape[0] - 1, 0, A.shape[1] - 1, 0, B.shape[1] - 1) return C |
Метод matrix_mult_limits() в каждом вызове передает матрицы A и B целиком, сохраняя результат в третьем параметре - матрице C размерности р*r. Остальные параметры - это нижние и верхние границы перемножаемых подматриц, связанные с размерностями p, q и r. Начальное условие в matrix_mult_limits() выполняется, когда обе подматрицы, скажем A[i,j] и B[j,k], - скаляры. В этом случае метод просто сохраняет их произведение в строке i и столбце k матрицы C. Метод не является функцией, так как он не возвращает результат (матрицу). Вместо этого он изменяет параметр C, где хранится результат. Наконец, итерационный метод add_matrices_limits() складывает элементы подматриц, переданных в двух первых входных параметрах, и сохраняет результат в своём третьем параметре (подматрицы задаются параметрами-границами).
На следующем шаге мы рассмотрим алгоритм Штрассена умножения матриц.
