标量向量化和矩阵乘法 [英] Vectorization and matrix multiplication by scalars

问题描述

我是python/numpy的新手. 我需要进行以下计算: 对于离散时间t的数组，对于$ 2 \乘以2 $的矩阵$ A $

I am new to python/numpy. I need to do the following calculation: for an array of discrete times t, calculate $e^{At}$ for a $2\times 2$ matrix $A$

我做了什么:

def calculate(t_,x_0,v_0,omega_0,c):

    # define A
    a_11,a_12, a_21, a_22=0,1,-omega_0^2,-c
    A =np.matrix([[a_11,a_12], [a_21, a_22]]) 
    print A
    # use vectorization 
    temps = np.array(t_)
    A_ = np.array([A  for k in  range (1,n+1,1)])
    temps*A_
    x_=scipy.linalg.expm(temps*A)
    v_=A*scipy.linalg.expm(temps*A)
    return x_,v_

n=10
omega_0=1
c=1
x_0=1
v_0=1
t_ = [float(5*k*np.pi/n)  for k in  range (1,n+1,1)]
x_, v_ = calculate(t_,x_0,v_0,omega_0,c)

但是，当我将A_(包含n乘以A的数组)和temps(包含要计算exp(At)的时间)相乘时会出现此错误:

However, I get this error when multiplying A_ (array containing n times A ) and temps (containg the times for which I want to calculate exp(At) :

ValueError:操作数不能与形状(10，)(10,2,2)一起广播

ValueError: operands could not be broadcast together with shapes (10,) (10,2,2)

据我理解向量化，A_中的每个元素将与临时索引相同的元素相乘；但我认为我做错了. 任何帮助/评论都非常感激

As I understand vectorization, each element in A_ would be multiplied by element at the same index from temps; but I think i don't get it right. Any help/ comments much appreciated

推荐答案

t_的纯numpy计算为(创建数组而不是列表):

A pure numpy calculation of t_ is (creates an array instead of a list):

In [254]: t = 5*np.arange(1,n+1)*np.pi/n
In [255]: t
Out[255]: 
array([ 1.57079633,  3.14159265,  4.71238898,  6.28318531,  7.85398163,
        9.42477796, 10.99557429, 12.56637061, 14.13716694, 15.70796327])

In [256]: a_11,a_12, a_21, a_22=0,1,-omega_0^2,-c
In [257]: a_11
Out[257]: 0
In [258]: A = np.array([[a_11,a_12], [a_21, a_22]]) 
In [259]: A
Out[259]: 
array([[ 0,  1],
       [-3, -1]])
In [260]: t.shape
Out[260]: (10,)
In [261]: A.shape
Out[261]: (2, 2)
In [262]: A_ = np.array([A  for k in  range (1,n+1,1)])
In [263]: A_.shape
Out[263]: (10, 2, 2)

A_是np.ndarray.我也将A设置为np.ndarray；您的是np.matrix，但是您的A_仍然是np.ndarray. np.matrix只能是2d，而A_是3d.

A_ is np.ndarray. I made A a np.ndarray as well; yours is np.matrix, but your A_ will still be np.ndarray. np.matrix can only be 2d, where as A_ is 3d.

所以t * A将是数组元素的乘法，因此出现广播错误(10,) (10,2,2).

So t * A will be array elementwise multiplication, hence the broadcasting error, (10,) (10,2,2).

要正确执行元素乘法，您需要类似

To do that elementwise multiplication right you need something like

In [264]: result = t[:,None,None]*A[None,:,:]
In [265]: result.shape
Out[265]: (10, 2, 2)

---

但是，如果您想将(10，)与(10,2,2)进行矩阵乘法，则einsum可以轻松做到这一点:

---

But if you want matrix multiplication of the (10,) with (10,2,2), then einsum does it easily:

In [266]: result1 = np.einsum('i,ijk', t, A_)
In [267]: result1
Out[267]: 
array([[   0.        ,   86.39379797],
       [-259.18139392,  -86.39379797]])

np.dot无法做到这一点，因为它的规则是倒数第二个倒数第二个". tensordot可以，但是我更喜欢einsum.

np.dot can't do it because its rule is 'last with 2nd to last'. tensordot can, but I'm more comfortable with einsum.

但是einsum表达式使(对我而言)很明显，我可以通过在第一轴上求和来从元素*中得到相同的东西:

But that einsum expression makes it obvious (to me) that I can get the same thing from the elementwise *, by summing on the 1st axis:

In [268]: (t[:,None,None]*A[None,:,:]).sum(axis=0)
Out[268]: 
array([[   0.        ,   86.39379797],
       [-259.18139392,  -86.39379797]])

或(t[:,None,None]*A[None,:,:]).cumsum(axis=0)每次获得2x2.

这篇关于标量向量化和矩阵乘法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！