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

Question

假设我有一个函数，它的输出是一个数字（标量），但输入是一个向量（多个数字）。比如说：

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

那么，有什么好的方法可以找到这个函数的一个根（也就是让函数输出为零的输入）呢？scipy.optimize.fsolve 和 scipy.optimize.root 这两个工具希望你提供的函数返回的是一个向量，而不是单个数字；而 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])

这样它就能同时改变两个参数了。不过代码看起来还是有点复杂。

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]

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。