使用 numpy 线性代数计算二次响应

Question

我在 numpy.xml 中处理矩阵时遇到问题。我有一个包含 N 个输入的向量 theta，我想使用二次模型 theta^T A theta + theta^T b + c 计算标量响应，其中 ^T 表示转置，A 是 NxN 方阵，b 是一个 N 维向量，c 是一个标量。当 theta 是一个 (NxM) 矩阵时，这意味着我有 M 个 theta 值要传播，我必须计算 theta^T A theta 以产生一个 M 维矩阵。在索引符号中，计算是 theta_{mj} A_{ji} theta_{im}，其中 m 不求和。

如果我只有一组 theta 值，numpy 的线性代数会按预期工作（这里，N=10 和 M=1）：

In [1]: import numpy as np                                            

In [2]: theta = np.ones(10)                                           

In [3]: A = np.ones((10, 10))                                         

In [4]: b = np.ones(10)                                               

In [5]: theta.T.dot(A).dot(theta) + theta.T.dot(b) + 1                
Out[5]: 111.0

我认为操作 theta^T A theta 可能与这里的点积不同。当我将 theta 设为 NxM 矩阵时，我不明白有什么根本不同。我认为额外的维度会自然地贯穿这个代码，就像它对 b 和 c 术语一样。

如何让 numpy 从操作 theta^T A theta 返回一个 M 维数组？

我只能让它返回方阵。不幸的是，点积函数将此操作视为矩阵乘法（这里，N=10 和 M=5）：

In [6]: theta = np.ones((10, 5))
#       theta.T.dot(A).dot(theta) is equivalient to:
#             
#                 (M x N)           (N x N)  (N x M)
In [7]: np.matmul(theta.T, np.matmul(A,     theta))                      
Out[7]: 
array([[100., 100., 100., 100., 100.],
       [100., 100., 100., 100., 100.],
       [100., 100., 100., 100., 100.],
       [100., 100., 100., 100., 100.],
       [100., 100., 100., 100., 100.]]) 

相比之下，b 和 c 项自然带有额外的 theta 项并提供我想要的 M 维输出：

In [8]: theta.T.dot(b) + 1                                           
Out[8]: array([11., 11., 11., 11., 11.])

Answer 1

两种可能：

N,M = 10,5
A = np.random.randint(0,10,(N,N))
theta = np.random.randint(0,10,(N,M))
b = np.random.randint(0,10,N)

1) 利用 matmul 将 >2D 操作数视为堆栈这一事实：

(theta.T[:,None]@A@theta.T[...,None])[...,0,0] + b@theta + 1
# array([ 8188, 14837,  7697,  9719,  7262])

或

2) 使用einsum

np.einsum("ik,ij,jk->k",theta,A,theta) + b@theta + 1
# array([ 8188, 14837,  7697,  9719,  7262])

比较验证的一对一评估：

[t@A@t + b@t + 1 for t in theta.T]
# [8188, 14837, 7697, 9719, 7262]