优化函数，包括大矩阵操作

Question

我可能解释得不好，所以我将使用一个与我的问题有点相似的例子，但这里是:

我需要计算一个复杂的操作(重复)，它是一个只有几个标量的函数，比如x1, x2和x3。然后计算一个矩阵，其元素由x1, x2, x3(以及行和列位置)的函数决定，称为x。

中间步骤A、B和c中涉及几个非常大的矩阵积，这些积中涉及的矩阵的值是常数。

假设我们有大小为

的矩阵:

• A: k by n
• B: n by n
• C: n by k
• X: 3 by k

，其中n非常非常大，k非常大。

函数做(伪代码):

func = function(x1, x2, x3)
  X = make_X(x1, x2, x3)
return X * A * B * C * X^t

，其中*为常规点积。输出将只是一个3 × 3的矩阵!

我(可能天真地)认为一定有某种自动的方法可以有效地将其编译成本质上的9个x1, x2, x3的函数——每个输出矩阵的元素一个。

我经常使用numpy/scipy，但对sympy或theano没有经验，尽管它们似乎在大致范围内。对于如何解决这个问题有什么建议吗?

注:，任何解决这个问题的代码包最好是用python编写的，但它们不必是这样，只要它们可以从python调用即可。

Answer 1

用einsum术语表达问题可能会有所帮助:

np.einsum('ij,jk,kl,lm,nm->in', X, A, B, C, X)

可以分成两步:

ABC = np.einsum('jk,kl,lm->jm', A, B, C) # k by k
np.einsum('ij,jm,nm->in', X, ABC, X)

所以结果的i,j元素是:

R[i,j] = np.einsum('j,jm,m->', X[i,:], X[j,:])

和新的@操作符(在这种情况下只是点的操作符版本)

R = X@A@B@C@(X.T)

对于k,n=10,20，最后一个是最快的。

如果你是为x1,x2,x3的不同组合做这个，但是一组A,B,C，那么做ABC=A@B@C应该是一个节省时间的方法。但我怀疑将X@ABC@(X.T)分解成9个步骤的价值。ABC是kxk，所以您已经完成了涉及B的较大计算。

Answer 2

既然您提到了SymPy，那么您可以使用

一次轻松地计算整个内容。

x1, x2, x3 = symbols('x1, x2, x3')
X = Matrix(...)
A = Matrix(...)
B = Matrix(...)
C = Matrix(...)

(将...替换为实际的矩阵项)。然后计算乘积

result = X*A*B*C*X.T

然后可以在这个矩阵上使用lambdify将其转换为一个可以与numpy或scipy一起使用的数字函数

func = lambdify([x1, x2, x3], result)