Credits for the image go to Charchithowitzer.
We have seen a lot of algorithms for matrix multiplication. Some are slow, like brute-force, in which we simply solve our problem with polynomial time. We also have fast algorithms using dynamic programming. Here we will use a memoization technique based on a divide and conquer approach. This article will focus on Strassen’s multiplication recursive algorithm for multiplying nxn matrices, which is a little faster than the simple brute-force method. The time complexity of this algorithm is O(n^(2.8), which is less than O(n^3).
Matrix multiplication is based on a divide and conquer-based approach. Here we divide our matrix into a smaller square matrix, solve that smaller square matrix and merge into larger results. For larger matrices this approach will continue until we recurse all the smaller sub matrices.
Suppose we have two matrices, A and B, and we want to multiply them to form a new matrix, C.
C=AB, where all A,B,C are square matrices. We will divide these larger matrices into smaller sub matrices n/2; this will go on.
Credits for the image go to Mohit Gupta.
Now from above we see:
r=ae+bg
s=af+bh
t=ce+dg
u=cf+dh
Each of the above four equations satisfies two multiplications of n/2Xn/2 matrices and addition of their n/2xn/2 products. Using these equations to define a divide and conquer strategy we can get the relation among them as:
T(N) = 8T(N/2) + O(N2)
From the above we see that simple matrix multiplication takes eight recursion calls.
T(n)=O(n^3)
Thus, this method is faster than the ordinary one.
It takes only seven recursive calls, multiplication of n/2xn/2 matrices and O(n^2) scalar additions and subtractions, giving the below recurrence relations.
T(N) = 7T(N/2) + O(N2)
Steps of Strassen’s matrix multiplication:
Divide the matrices A and B into smaller submatrices of the size n/2xn/2.
Using the formula of scalar additions and subtractions compute smaller matrices of size n/2.
Recursively compute the seven matrix products Pi=AiBi for i=1,2,…7.
Now compute the r,s,t,u submatrices by just adding the scalars obtained from above points.
Submatrix Products:
We have read many times how two matrices are multiplied. We do not exactly know why we take the row of one matrix A and column of the other matrix and multiply each by the below formula.
Pi=AiBi
=(α1ia+α2ib+α3ic)(β1ie+β2if+β2ig)
Where a,b ,β,α are the coefficients of the matrix that we see here, the product is obtained by just adding and subtracting the scalar.
Credits for the image go to Mohit Gupta.
O(N2)
Below is the code implementation using Python, divided into two parts. In the first part we split our matrices into smaller matrices and in other functions we perform Strassen’s method of operation, which we see in the above formula of scalar addition and subtractions of the scalar.
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
import numpy as np
def split(mat):
r, c = mat.shape
r2, c2 = r //2, c//2
return mat[: r2,: cl2], mat[: r2, cl2: ], mat[r2: ,: c2], mat[r2: , c2: ]
def strassen(x, y):
if len(x) == 1:
return x * y
a, b, c, d = split(x)
e, f, g, h = split(y)
p1 = strassen(a, f - h)
p2 = strassen(a + b, h)
p3 = strassen(c + d, e)
p4 = strassen(d, g - e)
p5 = strassen(a + d, e + h)
p6 = strassen(b - d, g + h)
p7 = strassen(a - c, e + f)
c11 = p5 + p4 - p2 + p6
c12 = p1 + p2
c21 = p3 + p4
c22 = p1 + p5 - p3 - p7
c = np.vstack((np.hstack((c11, c12)), np.hstack((c21, c22))))
return c
From the above we see that Strassen’s method is looking good, but it is good to remember is that it is not often used for matrix multiplication for five reasons:
Constraints used in the above method are generally very large; it takes more time in computation.
When the matrix is sparse this method works fine because sparse matrices take less time to compute.
It is not practically possible as it is computation and theoretical approach only.
It takes more space for storing sub matrices.
There is less chance of accuracy.