python - 无法直接应用网格网格来创建 Surfaceplot

Question

考虑以下示例:

from numpy import array, exp, linspace, meshgrid
import matplotlib.pyplot as plt

x = array([1, 2, 3, 4, 5, 6])
y = array([10, 5.49, 0.89, -0.14, -1.07, 0.84])
a1, a2 = 1, 4

def model(x, a1, a2):
    return a1 * exp(-a2 * x)

def residuen(x, y, a1, a2):
    return modell(x, a1, a2) - y

def S(x, y, a1, a2): 
    return 0.5 * sum(residuen(x, y, a1, a2) ** 2)

输入值x和y是固定的。我想针对 a1、a2 的不同值计算函数 S 的值，并将结果显示为 3D 曲面图。在 matplotlib tutorial ，我找到了一个使用 meshgrid 函数的示例。我尝试将此功能应用于我的问题:

a1 = linspace(-100, 100, 1000)
a2 = linspace(-100, 100, 1000)
A1, A2 = meshgrid(a1, a2)
Z = S(x, y, A1, A2) 

当我运行此代码时，出现错误

ValueError: operands could not be broadcast together with shapes (1000,1000) (6,)

很明显为什么上面的代码不起作用，但我不知道如何修复/重新定义函数 S 以便我可以应用 meshgrid 函数来获取曲面图的所需值 Z。有人知道如何解决这个问题吗？

Answer 1

正如您的错误所示，您的问题实际上与 matplotlib 无关。您正在尝试计算模型与 1e6 个位置中每个位置的 6 个数据点之间的残差。由于残差是一个总和，您可以想象在某个时刻有一个 1000x1000x6 的数组，并将其在最后一个维度上求和。您还可以创建一个 6x1000x1000 数组并在第一个维度上求和，但我将展示前一种情况。

您要使用的技术称为 broadcasting 。这意味着从最后一个尺寸开始排列形状，并填写单元尺寸以匹配。引入广播的最简单的更改是 model :

a1[..., None] * exp(-a2[..., None] * x)

... 中的省略号 ( [..., None] )意思是“获取所有现有的维度”。它是“doing [:, :]”的简写。在本例中，但可以根据特定数组的需要多次。 None相当于 np.newaxis ，表示“添加单位尺寸”。总而言之，索引表达式是将尺寸附加到形状的惯用方法。另一种方法是执行类似 a1.reshape(a1.shape + (1,)) 的操作.

现在您可以进行以下操作:

1. 乘以形状 (1000, 1000, 1) (扩展 a2 )按形状 (6,) ( x ) -> 结果形状: (1000, 1000, 6) 。 a2的最后一个维度被广播。 x给出两个隐式前导维度，大小为 1 ，也广播。

   这是您发生错误的地方。您不能将 (1000, 1000) 形状的数组相乘和(6,)以一种有意义的方式逐元素进行，但您可以使用新的单位维度。

2. 形状指数 (1000, 1000, 6)数组保持形状。
3. 将结果乘以形状 (1000, 1000, 1) (扩展a1)。最后一个维度是广播。

现在residuen中的残差完全没问题:modell(x, a1, a2) - y查找 (1000, 1000, 6) 之间的差异数组和 (6,)大批。最终尺寸完美对齐，剩下的由广播负责。

唯一的区别是你必须告诉 sum 由于数组不再平坦，因此要在哪个轴上进行操作:

0.5 * sum(residuen(x, y, a1, a2)**2, axis=1)

但这并不完全是正确的 RMS 计算。我本来希望看到

sum(residuen(x, y, a1, a2)**2, axis=-1)**0.5

查找数组 RMS 的更简单方法是使用 np.linalg.norm :

np.linalg.norm(residuen(x, y, a1, a2), axis=-1)

TL;DR

您的程序需要使用广播:

x = array([1, 2, 3, 4, 5, 6])
y = array([10, 5.49, 0.89, -0.14, -1.07, 0.84])
a1 = np.linspace(-100, 100, 1000)[None, :]
a2 = np.linspace(-100, 100, 1000)[:, None]

def model(x, a1, a2):
    return a1[..., None] * exp(-a2[..., None] * x)

def residuen(x, y, a1, a2):
    return model(x, a1, a2) - y

def S(x, y, a1, a2): 
    return np.linalg.norm(residuen(x, y, a1, a2), axis=-1)

Z = S(x, y, A1, A2)