为多个矩阵 numpy 应用每个矩阵单元格的函数

如何解决为多个矩阵 numpy 应用每个矩阵单元格的函数

我会尽量简化我的实际问题。假设我有 3 个形状相似的矩阵：

m1
array([[0.1,0.12,0.32],[0.22,0.15,0.16],[0.42,0.52,0.62],[0.44,0.35,0.36]])

m2
array([[0.12,0.144,0.384],[0.264,0.18,0.192],[0.504,0.624,0.744],[0.528,0.42,0.432]])

m3
array([[0.106,0.1272,0.3392],[0.2332,0.159,0.1696],[0.4452,0.5512,0.6572],[0.4664,0.371,0.3816]])

每个矩阵都有一个对应的日期，我将其转换为整数并将其称为增量，因此：

delta_m1  = 2
delta_m2 = 5
delta_m3 = 10

我想做的是对每个基于单元格的时间序列应用一个函数，这意味着如果我考虑位置 (0,0) 处的所有单元格，则 m1_cell 为 0.1，m2_cell 为 0.12，而 m3_cell是 0.106。所以现在我有一个基于单元格的时间序列，其中时间（x 轴）是增量 (2,5,10)，y 轴是单元格 (0.1,0.106) 的值。我想应用一个平滑函数（它需要两个一维数组并输出一维数组），它采用这个基于单元格的时间序列，这意味着采用 X=(2,10) 和 y=(0.1,0.106) 并输出一个新值每个日期/增量/矩阵，因此，例如，输出可以是 (0.09,0.1,0.8)，其中每个值对应于每个矩阵 (0,0) 处的单元格。我描述的整个过程只涉及一个单元格。我想将它应用于所有单元格。然后我可以根据结果创建新矩阵。我想我可以简单地遍历每个矩阵的每个单元格，并应用该函数，但这似乎效率不高。我怎样才能以 Pythonic 高效的方式实现这一目标？

编辑

我想在每个基于单元格的时间序列上应用的函数：

def local_regression(x0,X,Y,k):
    # add bias term
    x0 = np.r_[1,x0]
    X = np.c_[np.ones(len(X)),X]
    
    # fit model: normal equations with kernel
    xw = X.T * radial_kernel(x0,k)
    beta = np.linalg.pinv(xw @ X) @ xw @ Y
    
    # predict value
    return x0 @ beta

def radial_kernel(x0,k):
    return np.exp(np.sum((X - x0) ** 2,axis=1) / (-2 * k * k))

def smooth(x,y):
    return [local_regression(x[p],x,y,k) for p in range(len(x))]

例如：

y = np.array([2,10])
x = np.array([0.1,0.106])
smooth(x,y)
>>> array([5.14556944,6.34810085,5.50632924])

其中 x 是一个整数数组（尽管它可以更改为浮点数）而 y 是一个浮点数数组。 k 是函数的参数，这里 k=20

解决方法

广播的问题主要是因为您在最后一个函数中使用了列表理解。因此发生的情况是您的函数仅在最后一个维度（轴 = -1）上运行。这意味着即使使用 np.vectorize，您也被迫转置您的输入，然后将它们转回。

此外，因此您还必须在 np.vectorize 中明确提及签名。

检查下面的代码。我已经标记了我进行更改的区域。

#Stack all the matrices (because similar shape)
M = np.stack([m1,m2,m3])                   #(3,4,3)
D = np.stack([delta_m1,delta_m2,delta_m3]) #(3,)

##################### FUNCTION ######################
k=20

def local_regression(x0,X,Y,k):
    # add bias term
    x0 = np.r_[1,x0]
    X = np.c_[np.ones(len(X)),X]
    
    # fit model: normal equations with kernel
    xw = X.T * radial_kernel(x0,k)
    beta = np.linalg.pinv(xw @ X) @ xw @ Y
    
    # predict value
    return x0 @ beta

def radial_kernel(x0,k):
    return np.exp(np.sum((X - x0) ** 2,axis=1) / (-2 * k * k))

def smooth(x,y):
    return np.array([local_regression(x[p],x,y,k) for p in range(len(x))]) #<-----

smooth_v = np.vectorize(smooth,signature='(i),(i)->(i)') #<----
##################################################

smooth_v(D[None,None,:],M.transpose(2,1,0)).transpose(2,0)  #<----- 
#First transpose brings axis=0 to axis=-1,and the second transpose reverses it.

array([[[0.10713757,0.12856508,0.34284022],[0.23570265,0.16070635,0.17142011],[0.44997779,0.55711536,0.66425293],[0.4714053,0.37498149,0.38569525]],[[0.10847502,0.13017003,0.34712007],[0.23864505,0.16271253,0.17356004],[0.4555951,0.56407012,0.67254514],[0.4772901,0.37966258,0.39051008]],[[0.11017182,0.13220618,0.35254982],[0.242378,0.16525773,0.17627491],[0.46272164,0.57289345,0.68306527],[0.484756,0.38560136,0.39661854]]])