我有 4 个非线性方程，其中包含三个未知数X, Y, Z，我想求解。方程的形式为：
F(m) = X^2 + a(m)Y^2 + b(m)XYcosZ + c(m)XYsinZ
...其中a和是取决于四个方程中每个值b的常数。cF
解决这个问题的最佳方法是什么？
问问题
38173 次
2 回答
52
有两种方法可以做到这一点。
- 使用非线性求解器
- 线性化问题并在最小二乘意义上解决它
设置
因此，据我了解您的问题，您知道 F、a、b 和 c 在 4 个不同的点，并且您想要反转模型参数 X、Y 和 Z。我们有 3 个未知数和 4 个观察到的数据点，所以问题是过度确定的。因此，我们将以最小二乘的方式求解。
在这种情况下使用相反的术语更为常见，所以让我们翻转你的等式。代替：
F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ
让我们写：
F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)
我们知道F, X, Y, 和Z4 个不同的点（例如F_0, F_1, ... F_i）。
我们只是更改变量的名称，而不是方程本身。（这比其他任何事情都更便于我思考。）
线性解
实际上可以线性化这个方程。您可以轻松求解a^2、b^2、a b cos(c)和a b sin(c)。为了使这更容易一点，让我们再次重新标记事物：
d = a^2
e = b^2
f = a b cos(c)
g = a b sin(c)
现在等式要简单得多：F_i = d + e X_i + f Y_i + g Z_i. d对、e、f和进行最小二乘线性反演很容易g。然后我们可以得到a, b, 和cfrom：
a = sqrt(d)
b = sqrt(e)
c = arctan(g/f)
好的，让我们把它写成矩阵形式。我们将翻译 4 个观察结果（我们将编写的代码将接受任意数量的观察结果，但现在让我们保持具体）：
F_i = d + e X_i + f Y_i + g Z_i
进入：
|F_0| |1, X_0, Y_0, Z_0| |d|
|F_1| = |1, X_1, Y_1, Z_1| * |e|
|F_2| |1, X_2, Y_2, Z_2| |f|
|F_3| |1, X_3, Y_3, Z_3| |g|
或者：（F = G * m我是一名地球物理学家，所以我们使用G“格林函数”和m“模型参数”。通常我们也会使用d“数据”而不是F，。）
在 python 中，这将转换为：
def invert(f, x, y, z):
G = np.vstack([np.ones_like(x), x, y, z]).T
m, _, _, _ = np.linalg.lstsq(G, f)
d, e, f, g = m
a = np.sqrt(d)
b = np.sqrt(e)
c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
return a, b, c
非线性解决方案
scipy.optimize正如@Joe 建议的那样，您也可以使用 解决这个问题。最容易访问的函数scipy.optimize是scipy.optimize.curve_fit默认使用 Levenberg-Marquardt 方法。
Levenberg-Marquardt 是一种“爬山”算法（好吧，在这种情况下它会走下坡路，但无论如何都会使用该术语）。从某种意义上说，您对模型参数（默认情况下scipy.optimize全部observed - predicted在
警告：选择正确的非线性反演方法、初始猜测和调整方法的参数是非常“黑魔法”。你只能通过做来学习它，并且有很多情况下事情不会正常工作。如果您的参数空间相当平滑（这个应该是），Levenberg-Marquardt 是一种很好的通用方法。还有很多其他方法（包括遗传算法、神经网络等，除了模拟退火等更常见的方法）在其他情况下更好。我不打算在这里深入研究那部分。
有一个常见的问题是一些优化工具包试图纠正它scipy.optimize并没有尝试处理。如果您的模型参数具有不同的大小（例如a=1, b=1000, c=1e-8），您需要重新调整大小以使它们的大小相似。否则scipy.optimize，“爬山”算法（如 LM）将无法准确计算局部梯度的估计值，并且会给出非常不准确的结果。目前，我假设a、b和c具有相对相似的量级。另外，请注意，基本上所有非线性方法都需要您进行初始猜测，并且对该猜测很敏感。我将它留在下面（只需将其作为p0kwarg 传递给curve_fit），因为默认值a, b, c = 1, 1, 1是对a, b, c = 3, 2, 1.
排除警告后，curve_fit期望传递一个函数、一组进行观察的点（作为单个ndim x npoints数组）以及观察到的值。
所以，如果我们这样写函数：
def func(x, y, z, a, b, c):
f = (a**2
+ x * b**2
+ y * a * b * np.cos(c)
+ z * a * b * np.sin(c))
return f
在将它传递给curve_fit.
简而言之：
def nonlinear_invert(f, x, y, z):
def wrapped_func(observation_points, a, b, c):
x, y, z = observation_points
return func(x, y, z, a, b, c)
xdata = np.vstack([x, y, z])
model, cov = opt.curve_fit(wrapped_func, xdata, f)
return model
两种方法的独立示例：
为了给你一个完整的实现，这里有一个例子
- 生成随机分布的点来评估函数，
- 在这些点上评估函数（使用设置的模型参数），
- 给结果增加噪音，
- 然后使用上述线性和非线性方法对模型参数进行反演。
import numpy as np
import scipy.optimize as opt
def main():
nobservations = 4
a, b, c = 3.0, 2.0, 1.0
f, x, y, z = generate_data(nobservations, a, b, c)
print 'Linear results (should be {}, {}, {}):'.format(a, b, c)
print linear_invert(f, x, y, z)
print 'Non-linear results (should be {}, {}, {}):'.format(a, b, c)
print nonlinear_invert(f, x, y, z)
def generate_data(nobservations, a, b, c, noise_level=0.01):
x, y, z = np.random.random((3, nobservations))
noise = noise_level * np.random.normal(0, noise_level, nobservations)
f = func(x, y, z, a, b, c) + noise
return f, x, y, z
def func(x, y, z, a, b, c):
f = (a**2
+ x * b**2
+ y * a * b * np.cos(c)
+ z * a * b * np.sin(c))
return f
def linear_invert(f, x, y, z):
G = np.vstack([np.ones_like(x), x, y, z]).T
m, _, _, _ = np.linalg.lstsq(G, f)
d, e, f, g = m
a = np.sqrt(d)
b = np.sqrt(e)
c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
return a, b, c
def nonlinear_invert(f, x, y, z):
# "curve_fit" expects the function to take a slightly different form...
def wrapped_func(observation_points, a, b, c):
x, y, z = observation_points
return func(x, y, z, a, b, c)
xdata = np.vstack([x, y, z])
model, cov = opt.curve_fit(wrapped_func, xdata, f)
return model
main()
于 2013-10-23T14:35:21.663 回答
2
您可能想使用 scipy 的非线性求解器，它们非常简单：http ://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html
于 2013-10-23T14:15:31.617 回答