如何在 Python 中创建真正可调用的数组或矩阵
How to create a truly callable array or matrix in Python
我想制作一个矩阵，其所有条目都作为某个变量的函数 x。因此 B(x) 将以快速的方式理想地给出 N x N 输出。如果您愿意输入以函数作为条目的矩阵，这实际上是一项简单的任务。例如：
f1 = lambda x: 2*x
f2 = lambda x: x**2
B = lambda x : np.array([[f1(x),f2(x)],
[f2(x),f1(x)]])
这很天真，因为它无法扩展到数组很大且具有多种功能的场景。一个问题打出来就需要很长时间。人们通常会创建一个空数组并使用两个 Python for 循环来计算给定条目的特定函数，然后将输出存放在数组中。然后返回数组。
上述方法的问题在于每次调用函数时，它都是运行那些for循环。如果您想 运行 在 x 值的数据集上运行函数，这会使速度变慢。我尝试使用 Sympy 的 lambdfiy 函数创建静态可调用数组。对于 x 的评估，它似乎比纯 Python 中的 for 循环解决方案更快。但是，安装成本远远超过了这一点。有关详细信息，请参阅下面的代码。
有没有办法使用 Numpy 中的 vectorize 函数来加快速度？您能否找到比 for 循环版本更快的解决方案？
我也玩弄这个想法（或称之为梦想），其中可以评估 X 的整个数据集而不是单独评估每个 x。就像在 Numpy 中广播一样。例如
# Naive
result1 = [np.sin(x) for x in X]
# vs
results2 = np.sin(X)
总之，这太牵强了。
这是我写的代码。请尝试使用 N 的大小，看看速度下降有多迷人。澄清一下，我已经从整体上评估了我的程序，可调用数组的这个问题是瓶颈所在。
import numpy as np
from sympy import symbols,lambdify,zeros
from time import time
def get_sympy(N):
'''
Creates a callable array using Sympys lambdfiy capabilites.
This is truly a callable array.
'''
x = symbols('x')
output = zeros(N,N)
for i in range(N):
for j in range(N):
if i == j:
output[i,j] = x**2
elif i == 1:
output[i,j] = x**3
elif j == 0:
output[i,j] = x**4
else:
output[i,j] = x
return lambdify(x,output,'numpy')
def get_python(x,N):
'''
Uses Python loops to setup an array that mimic that of a callable array.
It is not truly a callable array as the loops run on each call of
this function.
'''
output = np.zeros((N,N))
f1 = lambda x: x**2
f2 = lambda x: x**3
f3 = lambda x: x**4
for i in range(N):
for j in range(N):
if i == j:
output[i,j] = f1(x)
elif i == 1:
output[i,j] = f2(x)
elif j == 0:
output[i,j] = f3(x)
else:
output[i,j] = x
return output
if __name__ == '__main__':
N = 30
X = np.random.uniform()
callable_sympy_array = get_sympy(N)
callable_python_array = lambda x: get_python(x,N)
t1 = time()
sympy_result = callable_sympy_array(X)
t2 = time()
python_result = callable_python_array(X)
t3 = time()
sympy_func = get_sympy(N)
t4 = time()
sympy_time = t2-t1
python_time = t3-t2
sympy_setup_time = t4-t3
print('\nSingle Evaluation Results:\n')
print('Sympy: ',round(sympy_time,5))
print('Python: ',round(python_time,5))
print('Sympy + Setup',round(sympy_setup_time,5))
evals = 100
print('\nResults for {0} evaluations of a {1} by {1} array:\n'.format(evals,N))
print('Sympy: ',sympy_setup_time + evals*sympy_time)
print('Python: ',python_time*evals)
快速 numpy 评估需要将内置编译的 operators/functions 应用于整个数组。任何类型的 python 级迭代都会减慢您的速度，就像在标量上评估（一般）Python 函数一样。快速的东西主要限于运算符（如 **）和 ufunc（np.sin 等）。
您的 sympy 生成的函数说明了这一点：
在 isympy 会话中：
In [65]: M = get_sympy(3)
使用 ipython 代码内省：
In [66]: M??
Signature: M(x)
Docstring:
Created with lambdify. Signature:
func(x)
Expression:
Matrix([[x**2, x, x], [x**3, x**2, x**3], [x**4, x, x**2]])
Source code:
def _lambdifygenerated(x):
return (array([[x**2, x, x], [x**3, x**2, x**3], [x**4, x, x**2]]))
Imported modules:
Source:
def _lambdifygenerated(x):
return (array([[x**2, x, x], [x**3, x**2, x**3], [x**4, x, x**2]]))
File: /<lambdifygenerated-8>
Type: function
所以这是 x 中的一个函数，使用 numpy 运算、** 运算符和数组创建。就像您输入它一样。sympy 通过在其符号代码中进行词法替换来创建它，因此您可以说它确实 'type it in'.
它可以对标量进行运算
In [67]: M(3)
Out[67]:
array([[ 9, 3, 3],
[27, 9, 27],
[81, 3, 9]])
在数组上，此处生成 (3,3,3) 结果：
In [68]: M(np.arange(1,4))
Out[68]:
array([[[ 1, 4, 9],
[ 1, 2, 3],
[ 1, 2, 3]],
[[ 1, 8, 27],
[ 1, 4, 9],
[ 1, 8, 27]],
[[ 1, 16, 81],
[ 1, 2, 3],
[ 1, 4, 9]]])
我希望很容易编写一个 sympy 表达式，该表达式在翻译时不能采用数组参数。 if 测试很难以数组兼容的形式编写，因为 Python if 表达式仅适用于标量布尔值。
您的 get_python 不会采用数组 x，主要是因为
output = np.zeros((N,N))
有固定尺寸；使用 np.zeros((N,N)+x.shape), x.dtype) 可能会解决这个问题。
无论如何，由于每次调用的 python 级迭代，它会很慢。
===
如果您尝试分配元素组会更快。例如在这种情况下：
In [76]: output = np.zeros((3,3),int)
In [77]: output[:] = 3
In [78]: output[:,0]=3**4
In [79]: output[1,:]=3**3
In [80]: output[np.arange(3),np.arange(3)]=3**2
In [81]: output
Out[81]:
array([[ 9, 3, 3],
[27, 9, 27],
[81, 3, 9]])
===
frompyfunc 是处理此类情况的便捷工具。在某些情况下，它提供比直接迭代快 2 倍的速度提升。但即使没有它也可以使代码更简洁。
例如，快速写下您的示例：
In [82]: def foo(x,i,j):
...: if i==j: return x**2
...: if i==1: return x**3
...: if j==0: return x**4
...: return x
In [83]: f = np.frompyfunc(foo, 3, 1)
In [84]: f(3,np.arange(3)[:,None], np.arange(3))
Out[84]:
array([[9, 3, 3],
[27, 9, 27],
[81, 3, 9]], dtype=object)
对于 Out[68] 示例：
In [98]: f(np.arange(1,4),np.arange(3)[:,None,None], np.arange(3)[:,None]).shape
Out[98]: (3, 3, 3)
In [99]: timeit f(np.arange(1,4),np.arange(3)[:,None,None], np.arange(3)[:,None]).shape
23 µs ± 471 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [100]: timeit M(np.arange(1,4))
21.7 µs ± 440 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
在标量 x 上进行评估，我的 f 与您的 get_python 的速度大致相同。
In [115]: MM = get_sympy(30)
In [116]: timeit MM(3)
109 µs ± 112 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [117]: timeit get_python(3,30)
241 µs ± 2.06 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [118]: timeit f(3,np.arange(30)[:,None], np.arange(30)).astype(int)
254 µs ± 1.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
我非常喜欢接受的答案。我只是想分享一个我自己的解决方案。
我使用的可调用矩阵示例是玩具示例。实际的可调用矩阵是已将 Logit 变换应用于其行的马尔可夫链。我还将针对某些参数采用这些矩阵的导数。无论如何，关键是我通过 Sympy 而不是手动完成所有这些。因此，对我的结果使用 lambdify 函数是有意义的。额外的好处是可以得到一个非常快速的可调用矩阵。
我只想提一下计算和 lambdify 过程是计算密集型的。但是我们有 cloudpickle 作为我们的朋友。
所以我所做的是为 range(2,500) 中 N 的每个维度计算我的可调用矩阵，然后将它们存储在一个序列化的大字典中。注意：cloudpickle 非常健壮，是 pickle、cpickle 或 dill 中唯一无需任何额外设置要求即可处理此问题的方法。
这是一个小例子：
from cloudpickle import dump,load
callable_arrays = dict()
for N in range(2,500):
callable_arrays[N] = get_sympy(N)
# Serialize the dictionary
with open('callable_array_file','wb') as file:
dump(callable_arrays,file)
# We can write a re-usable function to access callable arrays
def get_callable_array(N):
output = None
with open('callable_array_file','rb') as file:
output = load(file)[N]
return output
这可能可以进一步完善，但我对这个想法很满意。当天的惊喜确实在于 Sympy 可以生成可调用数组。 @hpaulji 给出的公认答案详细说明了为什么会这样。
我想制作一个矩阵，其所有条目都作为某个变量的函数 x。因此 B(x) 将以快速的方式理想地给出 N x N 输出。如果您愿意输入以函数作为条目的矩阵，这实际上是一项简单的任务。例如：
f1 = lambda x: 2*x
f2 = lambda x: x**2
B = lambda x : np.array([[f1(x),f2(x)],
[f2(x),f1(x)]])
这很天真，因为它无法扩展到数组很大且具有多种功能的场景。一个问题打出来就需要很长时间。人们通常会创建一个空数组并使用两个 Python for 循环来计算给定条目的特定函数，然后将输出存放在数组中。然后返回数组。
上述方法的问题在于每次调用函数时，它都是运行那些for循环。如果您想 运行 在 x 值的数据集上运行函数，这会使速度变慢。我尝试使用 Sympy 的 lambdfiy 函数创建静态可调用数组。对于 x 的评估，它似乎比纯 Python 中的 for 循环解决方案更快。但是，安装成本远远超过了这一点。有关详细信息，请参阅下面的代码。
有没有办法使用 Numpy 中的 vectorize 函数来加快速度？您能否找到比 for 循环版本更快的解决方案？
我也玩弄这个想法（或称之为梦想），其中可以评估 X 的整个数据集而不是单独评估每个 x。就像在 Numpy 中广播一样。例如
# Naive
result1 = [np.sin(x) for x in X]
# vs
results2 = np.sin(X)
总之，这太牵强了。
这是我写的代码。请尝试使用 N 的大小，看看速度下降有多迷人。澄清一下，我已经从整体上评估了我的程序，可调用数组的这个问题是瓶颈所在。
import numpy as np
from sympy import symbols,lambdify,zeros
from time import time
def get_sympy(N):
'''
Creates a callable array using Sympys lambdfiy capabilites.
This is truly a callable array.
'''
x = symbols('x')
output = zeros(N,N)
for i in range(N):
for j in range(N):
if i == j:
output[i,j] = x**2
elif i == 1:
output[i,j] = x**3
elif j == 0:
output[i,j] = x**4
else:
output[i,j] = x
return lambdify(x,output,'numpy')
def get_python(x,N):
'''
Uses Python loops to setup an array that mimic that of a callable array.
It is not truly a callable array as the loops run on each call of
this function.
'''
output = np.zeros((N,N))
f1 = lambda x: x**2
f2 = lambda x: x**3
f3 = lambda x: x**4
for i in range(N):
for j in range(N):
if i == j:
output[i,j] = f1(x)
elif i == 1:
output[i,j] = f2(x)
elif j == 0:
output[i,j] = f3(x)
else:
output[i,j] = x
return output
if __name__ == '__main__':
N = 30
X = np.random.uniform()
callable_sympy_array = get_sympy(N)
callable_python_array = lambda x: get_python(x,N)
t1 = time()
sympy_result = callable_sympy_array(X)
t2 = time()
python_result = callable_python_array(X)
t3 = time()
sympy_func = get_sympy(N)
t4 = time()
sympy_time = t2-t1
python_time = t3-t2
sympy_setup_time = t4-t3
print('\nSingle Evaluation Results:\n')
print('Sympy: ',round(sympy_time,5))
print('Python: ',round(python_time,5))
print('Sympy + Setup',round(sympy_setup_time,5))
evals = 100
print('\nResults for {0} evaluations of a {1} by {1} array:\n'.format(evals,N))
print('Sympy: ',sympy_setup_time + evals*sympy_time)
print('Python: ',python_time*evals)
快速 numpy 评估需要将内置编译的 operators/functions 应用于整个数组。任何类型的 python 级迭代都会减慢您的速度，就像在标量上评估（一般）Python 函数一样。快速的东西主要限于运算符（如 **）和 ufunc（np.sin 等）。
您的 sympy 生成的函数说明了这一点：
在 isympy 会话中：
In [65]: M = get_sympy(3)
使用 ipython 代码内省：
In [66]: M??
Signature: M(x)
Docstring:
Created with lambdify. Signature:
func(x)
Expression:
Matrix([[x**2, x, x], [x**3, x**2, x**3], [x**4, x, x**2]])
Source code:
def _lambdifygenerated(x):
return (array([[x**2, x, x], [x**3, x**2, x**3], [x**4, x, x**2]]))
Imported modules:
Source:
def _lambdifygenerated(x):
return (array([[x**2, x, x], [x**3, x**2, x**3], [x**4, x, x**2]]))
File: /<lambdifygenerated-8>
Type: function
所以这是 x 中的一个函数，使用 numpy 运算、** 运算符和数组创建。就像您输入它一样。sympy 通过在其符号代码中进行词法替换来创建它，因此您可以说它确实 'type it in'.
它可以对标量进行运算
In [67]: M(3)
Out[67]:
array([[ 9, 3, 3],
[27, 9, 27],
[81, 3, 9]])
在数组上，此处生成 (3,3,3) 结果：
In [68]: M(np.arange(1,4))
Out[68]:
array([[[ 1, 4, 9],
[ 1, 2, 3],
[ 1, 2, 3]],
[[ 1, 8, 27],
[ 1, 4, 9],
[ 1, 8, 27]],
[[ 1, 16, 81],
[ 1, 2, 3],
[ 1, 4, 9]]])
我希望很容易编写一个 sympy 表达式，该表达式在翻译时不能采用数组参数。 if 测试很难以数组兼容的形式编写，因为 Python if 表达式仅适用于标量布尔值。
您的 get_python 不会采用数组 x，主要是因为
output = np.zeros((N,N))
有固定尺寸；使用 np.zeros((N,N)+x.shape), x.dtype) 可能会解决这个问题。
无论如何，由于每次调用的 python 级迭代，它会很慢。
===
如果您尝试分配元素组会更快。例如在这种情况下：
In [76]: output = np.zeros((3,3),int)
In [77]: output[:] = 3
In [78]: output[:,0]=3**4
In [79]: output[1,:]=3**3
In [80]: output[np.arange(3),np.arange(3)]=3**2
In [81]: output
Out[81]:
array([[ 9, 3, 3],
[27, 9, 27],
[81, 3, 9]])
===
frompyfunc 是处理此类情况的便捷工具。在某些情况下，它提供比直接迭代快 2 倍的速度提升。但即使没有它也可以使代码更简洁。
例如，快速写下您的示例：
In [82]: def foo(x,i,j):
...: if i==j: return x**2
...: if i==1: return x**3
...: if j==0: return x**4
...: return x
In [83]: f = np.frompyfunc(foo, 3, 1)
In [84]: f(3,np.arange(3)[:,None], np.arange(3))
Out[84]:
array([[9, 3, 3],
[27, 9, 27],
[81, 3, 9]], dtype=object)
对于 Out[68] 示例：
In [98]: f(np.arange(1,4),np.arange(3)[:,None,None], np.arange(3)[:,None]).shape
Out[98]: (3, 3, 3)
In [99]: timeit f(np.arange(1,4),np.arange(3)[:,None,None], np.arange(3)[:,None]).shape
23 µs ± 471 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [100]: timeit M(np.arange(1,4))
21.7 µs ± 440 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
在标量 x 上进行评估，我的 f 与您的 get_python 的速度大致相同。
In [115]: MM = get_sympy(30)
In [116]: timeit MM(3)
109 µs ± 112 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [117]: timeit get_python(3,30)
241 µs ± 2.06 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [118]: timeit f(3,np.arange(30)[:,None], np.arange(30)).astype(int)
254 µs ± 1.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
我非常喜欢接受的答案。我只是想分享一个我自己的解决方案。
我使用的可调用矩阵示例是玩具示例。实际的可调用矩阵是已将 Logit 变换应用于其行的马尔可夫链。我还将针对某些参数采用这些矩阵的导数。无论如何，关键是我通过 Sympy 而不是手动完成所有这些。因此，对我的结果使用 lambdify 函数是有意义的。额外的好处是可以得到一个非常快速的可调用矩阵。
我只想提一下计算和 lambdify 过程是计算密集型的。但是我们有 cloudpickle 作为我们的朋友。
所以我所做的是为 range(2,500) 中 N 的每个维度计算我的可调用矩阵，然后将它们存储在一个序列化的大字典中。注意：cloudpickle 非常健壮，是 pickle、cpickle 或 dill 中唯一无需任何额外设置要求即可处理此问题的方法。
这是一个小例子：
from cloudpickle import dump,load
callable_arrays = dict()
for N in range(2,500):
callable_arrays[N] = get_sympy(N)
# Serialize the dictionary
with open('callable_array_file','wb') as file:
dump(callable_arrays,file)
# We can write a re-usable function to access callable arrays
def get_callable_array(N):
output = None
with open('callable_array_file','rb') as file:
output = load(file)[N]
return output
这可能可以进一步完善，但我对这个想法很满意。当天的惊喜确实在于 Sympy 可以生成可调用数组。 @hpaulji 给出的公认答案详细说明了为什么会这样。