如何用 python 中的较小矩阵项向量化填充较大矩阵
How to vectorize populating larger matrix with items of smaller matrix in python
我有一些小的对称矩阵，它们是较大对称矩阵的低维表示。我有一个向量，它是显示高维矩阵的哪些单元格应链接到低维矩阵中的哪些单元格的关键。
我想通过用低维矩阵中的相应值填充较大矩阵来重新创建这些较大矩阵。我相信应该有一个矢量化的方法来解决这个问题，但到目前为止我能想到的只是一个简单的嵌套 for 循环，这对于这些矩阵（10k+ 行和列）来说太慢了。
在这个玩具示例中，键是 vec1，低维矩阵是 source_mat，高维矩阵是 target_mat。我需要创建 target_mat，其中每个单元格都根据键从 source_mat 中填充相应的值。
import pandas as pd
import numpy as np
import random
vec1=[]
for x in range (0, 100):
vec1.append(random.randint(0, 19)) #creating the key
vec1=pd.DataFrame(vec1)
sizevec1=vec1.shape[0]
matshape=(sizevec1,sizevec1)
target_mat=np.zeros(matshape) #key and target have same shape
target_mat=pd.DataFrame(target_mat)
temp=np.random.random((20,20))
source_mat=temp*temp.T
for row in range(0,target_mat.shape[0]):
for column in range(0,target_mat.shape[1]):
print 'row is ', row
print 'column is', column
target_mat.iloc[row,column] = source_mat.item(int(vec1.iloc[row]), int(vec1.iloc[column]))
我设法想出了一个解决方案，它提供了相当大的加速，特别是对于较大的矩阵。这依赖于遍历较小的矩阵并用它的匹配元素填充大矩阵。
我用 vec1 作为具有 1000 个元素的向量尝试了这个解决方案，发现比以前的方法有 100 倍的加速。
import random
import time
import pandas as pd
import numpy as np
vec1=[]
for x in range (0, 1000):
vec1.append(random.randint(0, 19))
vec1=pd.DataFrame(vec1)
sizevec1=vec1.shape[0]
matshape=(sizevec1,sizevec1)
target_mat=np.zeros(matshape)
target_mat=pd.DataFrame(target_mat)
temp=np.random.random((20,20))
source_mat=temp*temp.T
###Slow Method###
matrixtime = time.time()
for row in range(0,target_mat.shape[0]):
for column in range(0,target_mat.shape[1]):
#print 'row is ', row
#print 'column is', column
target_mat.iloc[row,column] = source_mat.item(int(vec1.iloc[row]), int(vec1.iloc[column]))
print((time.time() - matrixtime))
target_mat_slow=target_mat
###FasterMethod###
target_mat=np.zeros(matshape)
target_mat=pd.DataFrame(target_mat)
matrixtime = time.time()
for row in range(0,source_mat.shape[0]):
for column in range(0, source_mat.shape[1]):
rowmatch = np.array(vec1==row)
rowmatch = rowmatch*1
colmatch = np.array(vec1==column)
colmatch = colmatch*1
match_matrix=rowmatch*colmatch.T
target_mat=target_mat+(match_matrix*source_mat[row,column])
print((time.time() - matrixtime))
target_mat_fast=target_mat
#Test Equivalence
target_mat_slow==target_mat_fast
这比您的 "fast" 答案快 3 倍。
import random
import time
import numpy as np
vec1=[]
for x in range (0, 1000):
vec1.append(random.randint(0, 19))
vec1=np.array(vec1)
sizevec1=vec1.shape[0]
matshape=(sizevec1,sizevec1)
target_mat=np.zeros(matshape)
temp=np.random.random((20,20))
source_mat=temp*temp.T
###FasterMethod###
target_mat=np.zeros(matshape)
def matrixops(vec1, source_mat, target_mat):
matrixtime = time.time()
for row in range(0,source_mat.shape[0]):
for column in range(0, source_mat.shape[1]):
rowmatch = np.array(vec1==row)
rowmatch = rowmatch*1
colmatch = np.array(vec1==column)
colmatch = colmatch*1
match_matrix=rowmatch*colmatch.T
target_mat=target_mat+(match_matrix*source_mat[row,column])
print((time.time() - matrixtime))
if __name__ == "__main__":
matrixops(vec1, source_mat, target_mat)
你的快版时间：4.246443033218384
本次版本时间：1.4500105381011963
正如我的评论所说，Cython 版本一点也不快。让它更快的唯一方法是采用依赖于 Python GIL 的行并转换为 C++ 风格的操作（就像我对 == 部分所做的那样，编写一个类似 C++ 的循环与 NumPy 函数相同，但 MemoryViews 不支持。由于我花了很多时间，所以在此发布以供参考：
cimport numpy
from numpy import array, multiply, asarray, ndarray, zeros, dtype, int
cimport cython
from cython cimport view
from cython.parallel cimport prange #this is your OpenMP portion
from openmp cimport omp_get_max_threads #only used for getting the max # of threads on the machine
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cpdef matrixops(int[::1] vec1, double[:,::1] source_mat, double[:,::1] target_mat):
cdef int[::1] match_matrix =zeros(vec1.shape[0], dtype=int)
cdef int[::1] rowmatch =zeros(vec1.shape[0], dtype=int)
cdef int[::1] colmatch =zeros(vec1.shape[0], dtype=int)
cdef int maxthreads = omp_get_max_threads()
cdef int row, column, i
# here's where you'd substitute
# for row in prange(source_mat.shape[0], nogil=True, num_threads=maxthreads, schedule='static'): # to use all cores
for row in range(0,source_mat.shape[0]):
for column in range(0, source_mat.shape[1]):
#this is how to avoid the GIL
for i in range(vec1.shape[0]):
rowmatch[i]=(row==vec1[i])
for i in range(vec1.shape[0]):
colmatch[i]=(column==vec1[i])
# this part has to be modified to not call Python GIL functions like was done above
match_matrix=multiply(rowmatch,colmatch.T)
target_mat=target_mat+(multiply(match_matrix,source_mat[row,column]))
上面就是你的 .PYX 文件。如果幸运的话，您通常会在 4 核上看到 3 倍的加速。抱歉，我没有成功，但是比你的 100 倍快的解决方案快 3 倍，使用直接 Python 库仍然不错。
下面是对代码的两个单独更新，它们带来了相当显着的加速。
首先- 找出向量化的解决方案，所以现在计算一步完成。即使在第二次更改之后，这也是最快的方法-
其次 - 将所有 pandas 数据帧更改为 numpy 数组。此更改对 for 循环代码的影响最大 - 现在运行速度提高了几个数量级。
下面的代码计算所有 3 种方法，'slow'、'fast' 和 'Xu Mackenzie'，以想出矢量化解决方案的朋友命名；-P
#初始化变量
import time
import random
import pandas as pd
import numpy as np
n=13000
k=2000
i=0
vec1=[]
for x in range(0, n):
vec1.append(random.randint(0, k-1))
temp=np.random.random((k,k))
#vec1=pd.DataFrame(vec1)
vec1=np.array(vec1)
#vec=pd.DataFrame(np.arange(0,300))
#vec2=pd.concat([vec,vec1], axis=1)
#sizevec1=vec1.shape[0]
sizevec1=len(vec1)
matshape=(sizevec1,sizevec1)
target_mat=np.zeros(matshape)
#target_mat=pd.DataFrame(target_mat)
source_mat=temp*temp.T
transform_mat=np.zeros((len(source_mat),len(target_mat)))
慢解
matrixtime = time.time()
for row in range(0,target_mat.shape[0]):
#print 'row is ', row
for column in range(0,target_mat.shape[1]):
#print 'column is', column
target_mat[row,column] = source_mat.item(int(vec1[row]), int(vec1[column]))
print((time.time() - matrixtime))
target_mat_slow=target_mat
target_mat=np.zeros(matshape)
XU MACKENZIE 解决方案
matrixtime = time.time()
for i in range(0,len(target_mat)):
transform_mat[vec1[i],i]=1
temp=np.dot(source_mat,transform_mat)
target_mat=np.dot(temp.T,transform_mat)
target_mat_XM=target_mat
target_mat=np.zeros(matshape)
XM_time= time.time() - matrixtime
print((time.time() - matrixtime))
上一个 'fast' 解决方案
matrixtime = time.time()
for row in range(0,source_mat.shape[0]):
print 'row is ', row
#for column in range(0, source_mat.shape[1]):
for column in range(0, row):
rowmatch = np.array([vec1==row])
rowmatch = rowmatch*1
colmatch = np.array([vec1==column])
colmatch = colmatch*1
match_matrix=rowmatch*colmatch.T
target_mat=target_mat+(match_matrix*source_mat[row,column])
print((time.time() - matrixtime))
target_mat_fast=target_mat
target_mat=np.zeros(matshape)
等效测试
target_mat_slow==target_mat_fast
target_mat_fast==target_mat_XM
我有一些小的对称矩阵，它们是较大对称矩阵的低维表示。我有一个向量，它是显示高维矩阵的哪些单元格应链接到低维矩阵中的哪些单元格的关键。
我想通过用低维矩阵中的相应值填充较大矩阵来重新创建这些较大矩阵。我相信应该有一个矢量化的方法来解决这个问题，但到目前为止我能想到的只是一个简单的嵌套 for 循环，这对于这些矩阵（10k+ 行和列）来说太慢了。
在这个玩具示例中，键是 vec1，低维矩阵是 source_mat，高维矩阵是 target_mat。我需要创建 target_mat，其中每个单元格都根据键从 source_mat 中填充相应的值。
import pandas as pd
import numpy as np
import random
vec1=[]
for x in range (0, 100):
vec1.append(random.randint(0, 19)) #creating the key
vec1=pd.DataFrame(vec1)
sizevec1=vec1.shape[0]
matshape=(sizevec1,sizevec1)
target_mat=np.zeros(matshape) #key and target have same shape
target_mat=pd.DataFrame(target_mat)
temp=np.random.random((20,20))
source_mat=temp*temp.T
for row in range(0,target_mat.shape[0]):
for column in range(0,target_mat.shape[1]):
print 'row is ', row
print 'column is', column
target_mat.iloc[row,column] = source_mat.item(int(vec1.iloc[row]), int(vec1.iloc[column]))
我设法想出了一个解决方案，它提供了相当大的加速，特别是对于较大的矩阵。这依赖于遍历较小的矩阵并用它的匹配元素填充大矩阵。
我用 vec1 作为具有 1000 个元素的向量尝试了这个解决方案，发现比以前的方法有 100 倍的加速。
import random
import time
import pandas as pd
import numpy as np
vec1=[]
for x in range (0, 1000):
vec1.append(random.randint(0, 19))
vec1=pd.DataFrame(vec1)
sizevec1=vec1.shape[0]
matshape=(sizevec1,sizevec1)
target_mat=np.zeros(matshape)
target_mat=pd.DataFrame(target_mat)
temp=np.random.random((20,20))
source_mat=temp*temp.T
###Slow Method###
matrixtime = time.time()
for row in range(0,target_mat.shape[0]):
for column in range(0,target_mat.shape[1]):
#print 'row is ', row
#print 'column is', column
target_mat.iloc[row,column] = source_mat.item(int(vec1.iloc[row]), int(vec1.iloc[column]))
print((time.time() - matrixtime))
target_mat_slow=target_mat
###FasterMethod###
target_mat=np.zeros(matshape)
target_mat=pd.DataFrame(target_mat)
matrixtime = time.time()
for row in range(0,source_mat.shape[0]):
for column in range(0, source_mat.shape[1]):
rowmatch = np.array(vec1==row)
rowmatch = rowmatch*1
colmatch = np.array(vec1==column)
colmatch = colmatch*1
match_matrix=rowmatch*colmatch.T
target_mat=target_mat+(match_matrix*source_mat[row,column])
print((time.time() - matrixtime))
target_mat_fast=target_mat
#Test Equivalence
target_mat_slow==target_mat_fast
这比您的 "fast" 答案快 3 倍。
import random
import time
import numpy as np
vec1=[]
for x in range (0, 1000):
vec1.append(random.randint(0, 19))
vec1=np.array(vec1)
sizevec1=vec1.shape[0]
matshape=(sizevec1,sizevec1)
target_mat=np.zeros(matshape)
temp=np.random.random((20,20))
source_mat=temp*temp.T
###FasterMethod###
target_mat=np.zeros(matshape)
def matrixops(vec1, source_mat, target_mat):
matrixtime = time.time()
for row in range(0,source_mat.shape[0]):
for column in range(0, source_mat.shape[1]):
rowmatch = np.array(vec1==row)
rowmatch = rowmatch*1
colmatch = np.array(vec1==column)
colmatch = colmatch*1
match_matrix=rowmatch*colmatch.T
target_mat=target_mat+(match_matrix*source_mat[row,column])
print((time.time() - matrixtime))
if __name__ == "__main__":
matrixops(vec1, source_mat, target_mat)
你的快版时间：4.246443033218384 本次版本时间：1.4500105381011963
正如我的评论所说，Cython 版本一点也不快。让它更快的唯一方法是采用依赖于 Python GIL 的行并转换为 C++ 风格的操作（就像我对 == 部分所做的那样，编写一个类似 C++ 的循环与 NumPy 函数相同，但 MemoryViews 不支持。由于我花了很多时间，所以在此发布以供参考：
cimport numpy
from numpy import array, multiply, asarray, ndarray, zeros, dtype, int
cimport cython
from cython cimport view
from cython.parallel cimport prange #this is your OpenMP portion
from openmp cimport omp_get_max_threads #only used for getting the max # of threads on the machine
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cpdef matrixops(int[::1] vec1, double[:,::1] source_mat, double[:,::1] target_mat):
cdef int[::1] match_matrix =zeros(vec1.shape[0], dtype=int)
cdef int[::1] rowmatch =zeros(vec1.shape[0], dtype=int)
cdef int[::1] colmatch =zeros(vec1.shape[0], dtype=int)
cdef int maxthreads = omp_get_max_threads()
cdef int row, column, i
# here's where you'd substitute
# for row in prange(source_mat.shape[0], nogil=True, num_threads=maxthreads, schedule='static'): # to use all cores
for row in range(0,source_mat.shape[0]):
for column in range(0, source_mat.shape[1]):
#this is how to avoid the GIL
for i in range(vec1.shape[0]):
rowmatch[i]=(row==vec1[i])
for i in range(vec1.shape[0]):
colmatch[i]=(column==vec1[i])
# this part has to be modified to not call Python GIL functions like was done above
match_matrix=multiply(rowmatch,colmatch.T)
target_mat=target_mat+(multiply(match_matrix,source_mat[row,column]))
上面就是你的 .PYX 文件。如果幸运的话，您通常会在 4 核上看到 3 倍的加速。抱歉，我没有成功，但是比你的 100 倍快的解决方案快 3 倍，使用直接 Python 库仍然不错。
下面是对代码的两个单独更新，它们带来了相当显着的加速。
首先- 找出向量化的解决方案，所以现在计算一步完成。即使在第二次更改之后，这也是最快的方法-
其次 - 将所有 pandas 数据帧更改为 numpy 数组。此更改对 for 循环代码的影响最大 - 现在运行速度提高了几个数量级。
下面的代码计算所有 3 种方法，'slow'、'fast' 和 'Xu Mackenzie'，以想出矢量化解决方案的朋友命名；-P
#初始化变量
import time
import random
import pandas as pd
import numpy as np
n=13000
k=2000
i=0
vec1=[]
for x in range(0, n):
vec1.append(random.randint(0, k-1))
temp=np.random.random((k,k))
#vec1=pd.DataFrame(vec1)
vec1=np.array(vec1)
#vec=pd.DataFrame(np.arange(0,300))
#vec2=pd.concat([vec,vec1], axis=1)
#sizevec1=vec1.shape[0]
sizevec1=len(vec1)
matshape=(sizevec1,sizevec1)
target_mat=np.zeros(matshape)
#target_mat=pd.DataFrame(target_mat)
source_mat=temp*temp.T
transform_mat=np.zeros((len(source_mat),len(target_mat)))
慢解
matrixtime = time.time()
for row in range(0,target_mat.shape[0]):
#print 'row is ', row
for column in range(0,target_mat.shape[1]):
#print 'column is', column
target_mat[row,column] = source_mat.item(int(vec1[row]), int(vec1[column]))
print((time.time() - matrixtime))
target_mat_slow=target_mat
target_mat=np.zeros(matshape)
XU MACKENZIE 解决方案
matrixtime = time.time()
for i in range(0,len(target_mat)):
transform_mat[vec1[i],i]=1
temp=np.dot(source_mat,transform_mat)
target_mat=np.dot(temp.T,transform_mat)
target_mat_XM=target_mat
target_mat=np.zeros(matshape)
XM_time= time.time() - matrixtime
print((time.time() - matrixtime))
上一个 'fast' 解决方案
matrixtime = time.time()
for row in range(0,source_mat.shape[0]):
print 'row is ', row
#for column in range(0, source_mat.shape[1]):
for column in range(0, row):
rowmatch = np.array([vec1==row])
rowmatch = rowmatch*1
colmatch = np.array([vec1==column])
colmatch = colmatch*1
match_matrix=rowmatch*colmatch.T
target_mat=target_mat+(match_matrix*source_mat[row,column])
print((time.time() - matrixtime))
target_mat_fast=target_mat
target_mat=np.zeros(matshape)
等效测试
target_mat_slow==target_mat_fast
target_mat_fast==target_mat_XM