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

Posted

技术标签:

【中文标题】使用 numpy 线性代数计算二次响应【英文标题】:Computing a quadratic response using numpy linear algebra 【发布时间】:2019-10-17 04:41:37 【问题描述】:

我在 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.])

【问题讨论】:

【参考方案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]

【讨论】:

谢谢!这里使用@ 符号的符号的名称是什么?我以前没见过。 是矩阵乘法运算符。对应的dunder是.__matmul__

以上是关于使用 numpy 线性代数计算二次响应的主要内容,如果未能解决你的问题,请参考以下文章

numpy linalg模块

科学计算基础软件包NumPy入门讲座:线性代数子模块

科学计算基础软件包NumPy入门讲座:线性代数子模块

python---Numpy模块中线性代数运算,统计和数学函数

numpy中的线性代数

numpy/scipy中非线性函数的数值梯度