标量向量化和矩阵乘法

Question

我是 python/numpy 的新手。我需要做以下计算：对于离散时间 t 的数组，为 $2\times 2$ 矩阵 $A$

计算 $e^{At}$

我做了什么：

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) 的时间）相乘时出现此错误：

ValueError：操作数无法与形状 (10,) (10,2,2) 一起广播

据我了解向量化，A_ 中的每个元素都将乘以来自 temps 的相同索引处的元素；但我想我没有弄清楚。非常感谢任何帮助/评论

Answer 1

这就是我要做的。

import numpy as np
from scipy.linalg import expm
A = np.array([[1, 2], [3, 4]])
for t in np.linspace(0, 5*np.pi, 20):
    print(expm(t*A))

此处未尝试矢量化。 expm 函数一次应用于一个矩阵，它肯定会占用大量计算时间。无需担心将 A 乘以标量的成本。

Answer 2

t_ 的纯 numpy 计算是（创建数组而不是列表）：

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。

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

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

的东西

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

但是如果你想要 (10,) 与 (10,2,2) 的矩阵乘法，那么 einsum 很容易做到：

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

np.dot做不到，因为它的规则是'last with 2nd to last'。 tensordot 可以，但我更喜欢 einsum。

但是那个 einsum 表达式很明显（对我来说）我可以通过在第一个轴上求和从元素 * 得到同样的东西：

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。

标量向量化和矩阵乘法

Vectorization and matrix multiplication by scalars

python

numpy

array-broadcasting