python - 从计算运算符/矩阵向量积的代码构造（稀疏）NumPy 数组？

Question

在科学计算中，我们经常需要构建计算微分算子的矩阵。编写应用运算符的代码通常比显式构造矩阵更容易。是否有一个库可以接收代码（假设它只使用线性运算）并输出矩阵，理想情况下是稀疏形式？

例子：

# computes finite differences
def myop(a):
  return a[1:]-a[:-1]

a = np.array(5)
myop(a)

计算“a”向量的有限差分. 我现在想要这样的东西

a = some_library.array(5)
op = myop(a)
print(op.as_matrix())

这应该给我矩阵表示：

[[-1, 1, 0, 0, 0],
 [0, -1, 1, 0, 0],
 [0, 0, -1, 1, 0],
 [0, 0, 0, -1, 1]]

拥有矩阵非常有用，例如用于计算转置算子或分析稀疏模式。从技术上讲，应该可以使用自动微分工具并提取 op() 的雅可比行列式，但我没有找到任何有效处理整个雅可比行列式的 AD 库，特别是如果它是稀疏的。他们似乎要么每行一次，要么每列一次雅可比行列式，即使只有几百个变量，这也非常慢。

Answer 1

一种简单的方法：

In [568]: arr = np.zeros((4,5),int)
In [569]: arr[np.arange(4),np.arange(4)]=-1
In [570]: arr[np.arange(4),np.arange(1,5)]=1
In [571]: arr
Out[571]: 
array([[-1,  1,  0,  0,  0],
       [ 0, -1,  1,  0,  0],
       [ 0,  0, -1,  1,  0],
       [ 0,  0,  0, -1,  1]])

有各种np.diag...功能，但由于您不想要方阵，我认为这会更容易。

该scipy.sparse库通常用于有限差分问题，并具有各种dia构造函数。但这需要更多的阅读。