python实现strassen算法

话不多说，先粘代码。可以直接运行，通过自行添加print语句可以进行矩阵运算结果的输出。

实现方法：

1.随机生成一定大小的方阵（矩阵）

2.为了实现strassen算法，编写矩阵相加，相减，矩阵拆分成四个子矩阵，矩阵合并

3.矩阵乘法常规算法和strassen算法

4.compare函数对不同大小的矩阵乘法进行测试并输出运行时间

5.为了验证计算正确，可以解除注释关于使用numpy的矩阵乘法部分

# coding=utf-8
"""
Created on 2019.2.9
@author: 刘祥
参考链接：https://blog.csdn.net/qq_32531519/article/details/81239763
"""
import numpy as np
from numpy import *
import time

class MatrixMultiplication(object):
	def __init__(self, size=10):
		self.size = size
		self.matrix_A = np.zeros((size,size), dtype="i").tolist()
		self.matrix_B = np.zeros((size,size), dtype="i").tolist()

	def generate_random_matrix(self, lower_bound=0, upper_bound=100):
		self.matrix_A = np.random.randint(lower_bound, upper_bound, (self.size, self.size))
		self.matrix_B = np.random.randint(lower_bound, upper_bound, (self.size, self.size))

	def direct_multiply(self):
		matrix_product = np.zeros((self.size,self.size), dtype="i").tolist()
		for row in range(0, self.size):
			for col in range(0, self.size):
				for k in range(0, self.size):
					matrix_product[row][col] += self.matrix_A[row][k]*self.matrix_B[k][col]
		return matrix_product
	
	def matrix_add(self, matrix_a, matrix_b):
		rows = len(matrix_a)
		columns = len(matrix_a[0])
		matrix_sum = np.zeros((rows,columns), dtype="i").tolist()
		for i in range(rows):
			for j in range(columns):
				matrix_sum[i][j] = matrix_a[i][j] + matrix_b[i][j]
		return matrix_sum

	def matrix_minus(self, matrix_a, matrix_b):
		rows = len(matrix_a)
		columns = len(matrix_a[0])
		matrix_difference = np.zeros((rows,columns), dtype="i").tolist()
		# matrix_difference = [[]*columns]*rows
		for i in range(rows):
			for j in range(columns):
				matrix_difference[i][j] = matrix_a[i][j] - matrix_b[i][j]
		return matrix_difference

        def matrix_divide(self, matrix):
		#返回四个子矩阵：分别为矩阵的左上，右上，左下，右下
		rows = len(matrix)
		columns = len(matrix[0])
		x_middle = rows // 2
		y_middle = columns // 2
		matrix_11 = [M[:x_middle] for M in matrix[:y_middle]]
		matrix_12 = [M[x_middle:] for M in matrix[:y_middle]]
		matrix_21 = [M[:x_middle] for M in matrix[y_middle:]]
		matrix_22 = [M[x_middle:] for M in matrix[y_middle:]]
		return matrix_11, matrix_12, matrix_21, matrix_22
	
	def matrix_merge(self, matrix_11, matrix_12, matrix_21, matrix_22):
		matrix_total = []
		rows1 = len(matrix_11)
		rows2 = len(matrix_21)
		for i in range(rows1):
			matrix_total.append(matrix_11[i]+matrix_12[i])
		for j in range(rows2):
			matrix_total.append(matrix_21[j]+matrix_22[j])
		return matrix_total
		
	def strassen(self, matrix_a, matrix_b):
		rows_a = len(matrix_a)
		cols_a = len(matrix_a[0])
		cols_b = len(matrix_b[0])
		matrix_product = np.zeros((rows_a,cols_b), dtype="i").tolist()#矩阵乘积的行列数设置
		if rows_a == 1:
			for i in range(cols_a):
				product = 0
				for j in range(cols_b):
					product += matrix_a[0][i] * matrix_b[i][j]
				matrix_product[0].append(product)
			return matrix_product
		else:
			a,b,c,d = self.matrix_divide(matrix_a)
			e,f,g,h = self.matrix_divide(matrix_b)
			P11 = self.matrix_minus(f,h)
			P21 = self.matrix_add(a,b)
			P31 = self.matrix_add(c,d)
			P41 = self.matrix_minus(g,e)
			P51 = self.matrix_add(a,d)
			P52 = self.matrix_add(e,h)
			P61 = self.matrix_minus(b,d)
			P62 = self.matrix_add(g,h)
			P71 = self.matrix_minus(a,c)
			P72 = self.matrix_add(e,f)
			
			P1 = self.strassen(a, P11)
			P2 = self.strassen(P21, h)
			P3 = self.strassen(P31, e)
			P4 = self.strassen(d, P41)
			P5 = self.strassen(P51, P52)
			P6 = self.strassen(P61, P62)
			P7 = self.strassen(P71, P72)
			
			r = self.matrix_add(self.matrix_minus(self.matrix_add(P5, P4), P2), P6)
			s = self.matrix_add(P1, P2)
			t = self.matrix_add(P3, P4)
			u = self.matrix_minus(self.matrix_minus(self.matrix_add(P5,P1),P3),P7)
			matrix_product = self.matrix_merge(r,s,t,u)
			return matrix_product
	
	
	def compare(self):
		size_list = [2**5, 2**6, 2**7,2**8, 2**9]
		for size in size_list:
			self.size = size
			self.generate_random_matrix()
			# print("矩阵A为：")
			# print(self.matrix_A)
			# print("矩阵B为：")
			# print(self.matrix_B)
			print("在数据量为%d时两种矩阵乘法用时如下："%size)
			
			start = time.clock()
			self.direct_multiply()
			end = time.clock()
			print("常规法用时:%.4f" % (end - start))
			
			start = time.clock()
			self.strassen(self.matrix_A, self.matrix_B)
			end = time.clock()
			print("strassen算法用时:%.4f" % (end - start))
		
			# print("numpy模块的矩阵运算结果为：")
			# # print(dot(mat(self.matrix_A), mat(self.matrix_B)))
			# print(mat(self.matrix_A) * mat(self.matrix_B))
		
		
		
if __name__ == '__main__':
	mat_mul = MatrixMultiplication()
	mat_mul.compare()

程序说明：

补充：在递归中递归边界很重要，也就是递归截止时应该怎么处理，

在矩阵乘法中当矩阵A的行数为1时已经不能继续划分为4个非空子矩阵了，所以得进行相应情况下的计算。

我的这部分代码有些冗余。因为strassen算法只能处理行列数为2的幂指数的矩阵，所以当递归到rows_a=1时，

cols_b也肯定为1。

                    if rows_a == 1:
			for i in range(cols_a):
				product = 0
				for j in range(cols_b):
					product += matrix_a[0][i] * matrix_b[i][j]
				matrix_product[0].append(product)
			return matrix_product

在编写程序的时候遇到一些小问题：

1.比如初始化行列数为rows，cols的二维数组（列表list）元素全为0

这里需要注意到列表的乘积代表重复,但是重复的过程其实是浅复制，具体内容请自行百度！

解决方法：

np.zeros((rows，cols), dtype="i").tolist()

dtype=“i”时生成整数数组,默认float

2.numpy模块的矩阵矢量相乘：

对于array对象，使用dot()或者matmul()函数

对于matrix对象，使用*号或dot()或matmul()函数

 

依然有很多不足，欢迎大家评论指教！