在python中的多维数组中进行插值，并将矩阵作为输出

Question

在我的问题中，我有一个函数在一个或多个参数上变化，并输出一个矩阵。现在对该函数进行采样。假设函数F依赖于x和y，那么对于某个x-y网格，我对每个网格点都有一个矩阵。

在MATLAB中，我使用4D矩阵来解决这个问题，其中前两个维度是矩阵本身，其余的是x-y网格的坐标。参见以下示例：

n = 2;
% x-y grid 
xg = [-1, 1, 7];
sz_xg = size(xg,2);
yg = [3, 4, 9];
sz_yg = size(yg,2);
% some function:
fcn = @(x,y) [3*x+sin(y), cos(x); exp(-x), 5+y]; % some n x n matrix function
% Define grid
A = zeros(n,n,sz_xg,sz_yg);
for i = 1:sz_xg
for j = 1:sz_yg
A(1:n, 1:n, i, j) = fcn(xg(i), yg(j));
end
end
xpoint = 5;
ypoint = 6;
Ainterp = interpn(1:n, 1:n, xg, yg, A, 1:n, 1:n, xpoint, ypoint)

我也尝试在Python中修复这个问题，但没有成功。当我试图实现完全相同的例子时，我得到了一个2x1的向量。。。

我的尝试：

import numpy as np
from scipy.interpolate import interpn
n = 2
xg = np.array([-1, 1, 7])
yg = np.array([3, 4, 9])
# some function
def fcn(x,y):
return np.array([[3*x+np.sin(y), np.cos(x)],[np.exp(-x), 5+y]])
A = np.zeros((n,n,xg.size,yg.size))
for i in range(xg.size):
for j in range(yg.size):
A[:,:,i,j]=fcn(xg[i],yg[j])
xpoint = 5
ypoint = 6
At = interpn((np.arange(n), np.arange(n), xg, yg),
A,(np.arange(0,n),np.arange(0,n),xpoint,ypoint))

你能帮我一下吗？

Answer 1

interpn的实现在Matlab和scipy之间略有不同。

Matlab返回输入向量的排列，因此当您为第一维度和第二维度指定[1，2]时，它会在处返回值

[(1,1) (1,2);
(2,1) (2,2)]

然而，scipy interpn的情况并非如此，scipy使用并行阵列来确定插值的位置，这意味着第一输出是由两个阵列中的第一输入构建的，第二输出是由这两个阵列的第二输入构建的等等

为了得到你想要的行为，你必须通过上面列出的整个矩阵，使用下面描述的更有效的方法：

grid_points = np.meshgrid(np.arange(0,n),np.arange(0,n),sparse=True)
At = interpn((np.arange(n), np.arange(n), xg, yg),
A,(grid_points[0],grid_points[1],xpoint,ypoint))

基本上，meshgrid(稀疏=True(将以一种高效且内存压缩的方式为您生成我上面写的矩阵。

在我看来，python方法比matlab更容易预测，因为它更容易确定输出数组的维度，而不像matlab，输出的维度是基于有多少输入不是一个单一的值，这使得很难预测输出中维度的顺序。

OP注意：当插值的输出是矢量时，您可能会出现错误。例如，At应该是3x1，并且它是用5的网格插值的3个变量。因此，A将为3x1x5x5x5。那么在误差计算中可能会出现误差。我通过在interpn函数之前添加来避免这种情况

np.seterr(divide='ignore', invalid='ignore')

并且在CCD_ 3功能之后

np.seterr(divide='warn', invalid='warn')