我正在寻找一种A_{ij}
从给定向量b_i
和另一个矩阵创建新矩阵的快速方法C_{ij}
。新矩阵的组件应具有以下形式
A_{ij} = b_i * C_{ij}.
到目前为止我正在使用dot(diag(b), C)
,但是点积自然有很多与零的乘法,这是非常低效的。有没有更好的办法?
使用*
,具有适当广播的元素乘积:
>>> b = array([1,2,3])
>>> C = arange(9).reshape(3,3)
>>> dot(diag(b), C)
array([[ 0, 1, 2],
[ 6, 8, 10],
[18, 21, 24]])
>>> atleast_2d(b).T * C
array([[ 0, 1, 2],
[ 6, 8, 10],
[18, 21, 24]])
atleast_2d(b).T
(或b.reshape(-1,1)
)将向量重塑b
为列向量。
您还可以使用einsum:
>>> b = np.array([1,2,3])
>>> C = np.arange(9).reshape(3,3)
>>> np.einsum('i,ij->ij', b, C)
array([[ 0, 1, 2],
[ 6, 8, 10],
[18, 21, 24]])
在我的机器上,这种方法比这两种方法都快
np.dot(np.diag(b), C)
和
np.atleast_2d(b).T * C
对于 100,000 个循环,einsum 需要 1.16 秒,dot 需要 1.61 秒,atleast_2d 需要 2.03 秒。
如果您熟悉Einstein summation notation,那么 einsum 是一个非常有用的函数。使用索引表示法更容易制定许多简单的线性代数运算,并且增加的执行速度只是锦上添花。
将矩阵转换为数组,执行(广播)数组乘法,然后转换回:
np.matrix(b * np.asarray(C))
做
b[:,np.newaxis] * C
或者
b[:,None] * C
结合已经给出的答案,并尝试一种新的(但效率低下的方法),我给出了以下时间安排。
如果我每十次重复 einsum 和广播(更准确的时间),我发现:
这是它的测试代码:
import timeit
from scipy.sparse.dia import dia_matrix
import numpy as np
def get_bc():
N = 10000
M = 10000
b = np.random.random(size=N)
c = np.random.random(size=[N, M])
return b, c
def multiply_dot_diag():
global b, c
return np.dot(np.diag(b), c)
def multiply_broadcast():
global b, c
return np.atleast_2d(b).T * c
def mulitply_sparse():
global b, c
N = b.size
return dia_matrix((b, 0), shape=[N, N]).todense().dot(c)
def multiply_einsum():
global b, c
return np.einsum('i,ij->ij', b, c)
if __name__ == '__main__':
global b, c
b, c = get_bc()
print (b, c)
print (timeit.timeit(multiply_einsum, number=1))
print (timeit.timeit(multiply_broadcast, number=1))
print (timeit.timeit(multiply_dot_diag, number=1))
print (timeit.timeit(mulitply_sparse, number=1))
np.testing.assert_array_equal(multiply_dot_diag(), multiply_broadcast())
np.testing.assert_array_equal(multiply_dot_diag(), mulitply_sparse())
np.testing.assert_array_equal(multiply_dot_diag(), multiply_einsum())