在 Python 中返回向量列表的向量化函数

Question

我正在尝试在 Python 中编写一个函数，该函数将参数列表 a 作为输入并返回依赖于这些参数的向量列表 u。然后我希望能够使用标准向量运算（例如标量积）对u 进行操作。在下面的示例中（灵感来自this post 的第二个答案），该函数返回两个初始向量v 和w 的组合：

import numpy as np
v = np.array([1,0,0])
w = np.array([0,1,0])
a = np.linspace(0,np.pi,4)
def u(a):
   u = np.cos(a)*v + np.sin(a)*w
   return u  
uvec = np.vectorize(u, otypes=[np.ndarray])    
ufin = np.array(uvec(a).tolist())
print "ufin =", ufin
print "ufin.v =", np.dot(ufin,v)

这会返回：

ufin = [[  1.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  5.00000000e-01   8.66025404e-01   0.00000000e+00]
 [ -5.00000000e-01   8.66025404e-01   0.00000000e+00]
 [ -1.00000000e+00   1.22464680e-16   0.00000000e+00]]
ufin.v = [ 1.   0.5 -0.5 -1. ]

这是我想要获得的，因为ufin 的行为就像一个向量。你能告诉我是否还有其他更简单的方法来实现这一点？我需要编写一个代码，其中必须定义大量向量和向量操作，并希望它尽可能紧凑。

提前谢谢你！

---

编辑：

我根据this post 的最后一个答案找到了另一个（显然更紧凑）解决方案。这个想法是将参数列表重塑为列数组，以便函数的输出自动（无需向量化）返回向量列表作为二维数组。这是通过这种方式完成的：

import numpy as np
from numpy.core.umath_tests import inner1d

v = np.array([1,0,0])
w = np.array([0,1,0])
a = np.linspace(0,np.pi,4).reshape((4,1))
b = np.linspace(0,np.pi/2,4).reshape((4,1))
def u(a):
   u = np.cos(a)*v + np.sin(a)*w
   return u 
print "u(a) =",u(a)
print "u(b) =",u(b) 
print "u(a).v =",np.dot(u(a),v)
print "u(a)^v =",np.cross(u(a),v)
# print "u(a).u(b) =",np.dot(u(a),u(b)) # does not work
print "u(a).u(b) =",inner1d(u(a),u(b)) # works
print "u(a)^u(b) =",np.cross(u(a),u(b)) 

这会返回：

u(a) = [[  1.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  5.00000000e-01   8.66025404e-01   0.00000000e+00]
 [ -5.00000000e-01   8.66025404e-01   0.00000000e+00]
 [ -1.00000000e+00   1.22464680e-16   0.00000000e+00]]
u(b) = [[  1.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  8.66025404e-01   5.00000000e-01   0.00000000e+00]
 [  5.00000000e-01   8.66025404e-01   0.00000000e+00]
 [  6.12323400e-17   1.00000000e+00   0.00000000e+00]]
u(a).v = [ 1.   0.5 -0.5 -1. ]
u(a)^v = [[  0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00  -8.66025404e-01]
 [  0.00000000e+00   0.00000000e+00  -8.66025404e-01]
 [  0.00000000e+00   0.00000000e+00  -1.22464680e-16]]
u(a).u(b) = [  1.00000000e+00   8.66025404e-01   5.00000000e-01   6.12323400e-17]
u(a)^u(b) = [[ 0.         0.         0.       ]
 [ 0.         0.        -0.5      ]
 [ 0.         0.        -0.8660254]
 [ 0.         0.        -1.       ]]

对于涉及初始（u 和 v）和向量输出列表（u(a) 和 u(b)）的操作以及涉及向量的两个输出列表的操作，这都是正确的行为。唯一需要注意的是（对于向量输出列表之间的操作）是必须使用inner1d 函数而不是标准np.dot，因为后者被解释为由于两个矩阵不一致而无法完成的矩阵乘积尺寸。

Answer 1

如果您不反对获取二维数组而不是对象数组（无论如何，对于大多数用途来说，前者比后者更方便），那么 外积 对此很有用问题：

>>> np.outer(np.cos(a), v) + np.outer(np.sin(a), w)
array([[  1.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  5.00000000e-01,   8.66025404e-01,   0.00000000e+00],
       [ -5.00000000e-01,   8.66025404e-01,   0.00000000e+00],
       [ -1.00000000e+00,   1.22464680e-16,   0.00000000e+00]])
>>> _ @ v
array([ 1. ,  0.5, -0.5, -1. ])

@ 中缀运算符仅适用于 python 3.5 或更高版本。在旧版本中，您可以使用np.dot。在这种情况下，任何一个的含义都是matrix multiplication。

Answer 2

所以你告诉vectorze，函数u返回一个数组，或者实际上是一个object，因为没有dtype=array。

In [475]: uvec=np.vectorize(u, otypes=[object])
In [476]: uvec(a)
Out[476]: 
array([array([ 1.,  0.,  0.]), array([ 0.5      ,  0.8660254,  0.       ]),
       array([-0.5      ,  0.8660254,  0.       ]),
       array([ -1.00000000e+00,   1.22464680e-16,   0.00000000e+00])], dtype=object)

因此，如果您想“解开”它并创建一个二维数组，则必须将其转换为数组列表，然后使用 np.array 或 (np.stack) 将它们沿新维度连接：

In [477]: uvec(a).tolist()
Out[477]: 
[array([ 1.,  0.,  0.]),
 array([ 0.5      ,  0.8660254,  0.       ]),
 array([-0.5      ,  0.8660254,  0.       ]),
 array([ -1.00000000e+00,   1.22464680e-16,   0.00000000e+00])]

但是对于像a 这样的单个一维数组的迭代，vectorize 是多余的。一个简单的列表推导也同样有效

In [483]: np.array([u(i) for i in a])
Out[483]: 
array([[  1.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  5.00000000e-01,   8.66025404e-01,   0.00000000e+00],
       [ -5.00000000e-01,   8.66025404e-01,   0.00000000e+00],
       [ -1.00000000e+00,   1.22464680e-16,   0.00000000e+00]]

vectorize 在有多个输入数组时更有用，它们可以有多个维度。它负责广播和迭代多个级别。你也可以使用np.frompyfunc(u,1,1)，它总是返回一个对象数组，而且速度稍快。

但你也可以不循环计算：

def ufunc(a,v,w):
    a = a[:,None]
    return np.cos(a)*v + np.sin(a)*w

In [498]: ufunc(a,v,w)
Out[498]: 
array([[  1.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  5.00000000e-01,   8.66025404e-01,   0.00000000e+00],
       [ -5.00000000e-01,   8.66025404e-01,   0.00000000e+00],
       [ -1.00000000e+00,   1.22464680e-16,   0.00000000e+00]])

这会将 a 转换为 (4,1) 数组，该数组针对 (3,) v 进行广播以生成 (4,3) 数组。

---

另一种方法是定义一个接受标量的函数（所有周围）

import math
def u1(a,v,w):
    return math.cos(a)*v+math.sin(a)*w
u1v = np.vectorize(u1, otypes=[float])

In [505]: u1v(a[:,None], v, w)
Out[505]: 
array([[  1.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  5.00000000e-01,   8.66025404e-01,   0.00000000e+00],
       [ -5.00000000e-01,   8.66025404e-01,   0.00000000e+00],
       [ -1.00000000e+00,   1.22464680e-16,   0.00000000e+00]])

这充分利用了vectorize 的广播能力。我虽然它会比你的代码慢，因为它必须在两个维度上循环，但它实际上测试得更快（虽然不如我的 ufunc 快）。尝试对 u 函数进行一半矢量化可能不会为您带来任何速度。

---

您还可以定义原始函数以将v 和w 作为参数，并使用vectorize excluded 参数按原样传递它们。我认为它没有性能优势。

In [521]: def u2(a,v,w):
     ...:     return np.cos(a)*v + np.sin(a)*w
     ...: 
In [522]: u2vec=np.vectorize(u2, otypes=[object],excluded=[1,2])
In [523]: u2vec(a,v,w)
Out[523]: 
array([array([ 1.,  0.,  0.]), array([ 0.5      ,  0.8660254,  0.       ]),
       array([-0.5      ,  0.8660254,  0.       ]),
       array([ -1.00000000e+00,   1.22464680e-16,   0.00000000e+00])], dtype=object)