无法直接应用 meshgrid 创建 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,)

很清楚为什么上面的代码不起作用，但我不知道如何 fix/redefine 函数 S 以便我可以应用 meshgrid 函数来获得所需的曲面图的 Z 值。有人知道如何解决这个问题吗？

Answer 1

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

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

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

[..., None] 中的省略号 (...) 表示 "take all the existing dimensions"。在这种情况下，它是 shorthand 用于执行 [:, :] ，但对于特定数组需要多次。 None 等同于 np.newaxis，意思是 "add a unit dimension"。总的来说，索引表达式是一种向形状附加维度的惯用方法。另一种方法是做类似 a1.reshape(a1.shape + (1,)).

的事情

现在你有以下操作：

形状 (1000, 1000, 1)（扩展 a2）乘以形状 (6,)（x）-> 结果形状：(1000, 1000, 6)。 a2 的最后一个维度是广播。 x 被赋予两个大小为 1 的隐式前导维度，它们也被广播。

这是您的错误发生的地方。您不能以有意义的方式逐元素地乘以 (1000, 1000) 和 (6,) 形状的数组，但可以使用新的单位维度。
形状的指数 (1000, 1000, 6) 数组保持形状。
将结果乘以形状 (1000, 1000, 1)（扩展 a1）。最后一个维度是 broadcasted.

现在 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)

无法直接应用 meshgrid 创建 Surfaceplot

Can't apply meshgrid directly to create Surfaceplot

python

numpy

matplotlib

array-broadcasting