使用SciPy找到多维标量函数的根的最佳方法

Question

假设我有一个函数，其值域为标量，但其定义域为向量。例如：

def func(x):
  return x[0] + 1 + x[1]**2

如何找到这个函数的根？scipy.optimize.fsolve和scipy.optimize.root期望func返回一个向量（而不是标量），并且scipy.optimize.newton只接受标量参数。我可以重新定义func：

def func(x):
  return [x[0] + 1 + x[1]**2, 0]

然后root和fsolve可以找到根，但是雅可比矩阵中的零意味着它不总是能够很好地完成任务。例如：

fsolve(func, array([0,2]))
=> array([-5,  2])

这将仅改变第一个参数而不是第二个参数，意味着它经常找到一个很远的零。

---

编辑：看起来以下对func的重新定义效果更好：

def func(x):
  fx = x[0] + 1 + x[1]**2
  return [fx, fx]

fsolve(func, array([0,5]))
=>array([-16.27342781,   3.90812331])

所以现在它愿意更改这两个参数。不过代码还是有点丑。

- lnmaurer

---

任何形如(-1 - y*2, y)的点都是根，因此询问“the* root”没有意义。在一般情况下，你应该期望f(x,y)=0的解集在(x,y)平面上是一条曲线。如果你想要唯一的解，你需要第二个函数或一个约束条件。 - Warren Weckesser

是的，如果你想要严谨一些，我正在寻找“根”——最好是相对于初始猜测比较接近的。也就是说，我正在寻找像 scipy.optimize.fsolve 或 scipy.optimize.root 这样的东西。 - lnmaurer

Answer 1

您尝试过使用fmin来最小化绝对值函数吗？例如：

>>> import scipy.optimize as op
>>> import numpy as np

>>> def func(x):
>>>     return x[0] + 1 + x[1]**2
>>> func1 = lambda x: np.abs(func(x))

>>> tmp = op.fmin(func1, [10000., 10000.])
>>> func(tmp)
0.0
>>> print tmp
[-8346.12025122    91.35162971]

- Taha

---

我个人不会使用这种解决方案，除非是某些类型的函数，因为起始点的选择很重要，而Jacobi方法并不连续。 如果您知道函数的形状并且它不是非常复杂的话，这可能是一个解决方案。 - Taha

我考虑过尝试最小化绝对值（或平方），但我认为寻找根可能更准确（并且可能更快），而不是最小化函数。 对于我要解决的实际问题，我有一个相当好的初始猜测，并且函数非常平滑，因此我实际上只需要使用牛顿（或割线）方法。 - lnmaurer

Answer 2

对于我的问题，由于有一个很好的初始猜测和一个非疯狂的函数，牛顿法很有效。对于标量、多维函数，牛顿法变成：

下面是一个简单的代码示例：

def func(x): #the function to find a root of
  return x[0] + 1 + x[1]**2

def dfunc(x): #the gradient of that function
  return array([1, 2*x[1]])

def newtRoot(x0, func, dfunc):
  x = array(x0)
  for n in xrange(100): # do at most 100 iterations
    f  = func(x)
    df = dfunc(x)

    if abs(f) < 1e-6: # exit function if we're close enough
      break

    x = x - df*f/norm(df)**2 # update guess
  return x

正在使用:

nsolve([0,2],func,dfunc)
=> array([-1.0052546 ,  0.07248865])

func([-1.0052546 ,  0.07248865])
=> 4.3788225025098715e-09

不错！当然，这个函数很粗糙，但你有了想法。它也不适用于“棘手”的函数或者你没有一个好的起始猜测。我认为我会使用类似这样的方法，但是如果牛顿法不收敛，我会退而求其次使用fsolve或root。

- lnmaurer