2

我正在寻找一种A_{ij}从给定向量b_i和另一个矩阵创建新矩阵的快速方法C_{ij}。新矩阵的组件应具有以下形式

A_{ij} = b_i * C_{ij}.

到目前为止我正在使用dot(diag(b), C)，但是点积自然有很多与零的乘法，这是非常低效的。有没有更好的办法？

4

5 回答 5

2

使用*，具有适当广播的元素乘积：

>>> 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为列向量。

于 2012-12-03T13:19:20.620 回答
1

您还可以使用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 是一个非常有用的函数。使用索引表示法更容易制定许多简单的线性代数运算，并且增加的执行速度只是锦上添花。

于 2013-03-08T03:30:45.030 回答
0

将矩阵转换为数组，执行（广播）数组乘法，然后转换回：

np.matrix(b * np.asarray(C))
于 2012-12-03T13:15:06.393 回答
0

做

b[:,np.newaxis] * C

或者

b[:,None] * C
于 2012-12-03T22:35:12.893 回答
0

结合已经给出的答案，并尝试一种新的（但效率低下的方法），我给出了以下时间安排。

  • einsum: 0.54 s (Thycydides)
  • 广播: 0.44 s (larsmans)
  • 点对角线：15.6 s (henrikr, TS)
  • 稀疏直径矩阵：16.8 s（自己尝试，惨败:-p）

如果我每十次重复 einsum 和广播（更准确的时间），我发现：

  • einsum：5.2 秒
  • 广播：4.5秒

这是它的测试代码：

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())
于 2015-02-20T06:27:24.303 回答