SciPy的optimize模块提供了许多数值优化算法，本节对其中的非线性方程组求解、数据拟合、函数最小值等进行简单介绍。

一、非线性方程组求解

fsolve( )可以对非线性方程组进行求解，它的基本调用形式为fsolve(func, x0)。其中func是计算方程组误差的函数，它的参数x是一个数组，其值为方程组的一组可能的解。func返回将x代入方程组之后得到的每个方程的误差，x0为未知数的一组初始值。假设要对下面的方程组进行求解：

f1(u1,u2,u3)=0, f2(u1,u2,u3)=0, f3(u1,u2,u3)=0

那么func函数可以定义如下：

def func(x):
u1,u2,u3 = x
return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]

下面我们看一个对下列方程组求解的例子：

from math import sin, cos
from scipy import optimize

def f(x): ?
x0, x1, x2 = x.tolist( ) ?
return [
5*x1+3,
4*x0*x0 - 2*sin(x1*x2),
x1*x2 - 1.5
]

# f计算方程组的误差，[1,1,1]是未知数的初始值
result = optimize.fsolve(f, [1,1,1]) ?
print result
print f(result)
[-0.70622057 -0.6 -2.5 ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]

?f( )是计算方程组的误差的函数，x参数是一组可能的解。fsolve( )在调用f( )时，传递给f( )的参数是一个数组。?先调用数组的tolist( )方法，将其转换为Python的标准浮点数列表，然后调用math模块中的函数进行运算。因为在进行单个数值的运算时，标准浮点类型比NumPy的浮点类型要快许多，所以把数值都转换成标准浮点数类型，能缩短一些计算时间。?调用fsolve( )时，传递计算误差的函数f( )以及未知数的初始值。

在对方程组进行求解时，fsolve( )会自动计算方程组在某点对各个未知数变量的偏导数，这些偏导数组成一个二维数组，数学上称之为雅可比矩阵。如果方程组中的未知数很多，而与每个方程有关联的未知数较少，即雅可比矩阵比较稀疏时，将计算雅可比矩阵的函数作为参数传递给fsolve( )，这能大幅度提高运算速度。笔者在一个模拟计算的程序中需要求解有50个未知数的非线性方程组。每个方程平均与6个未知数相关，通过传递计算雅可比矩阵的函数使fsolve( )的计算速度提高了4倍。

雅可比矩阵

雅可比矩阵是一阶偏导数以一定方式排列的矩阵，它给出了可微分方程与给定点的最优线性逼近，因此类似于多元函数的导数。例如前面的函数f1、f2、f3和未知数u1、u2、u3的雅可比矩阵如下：

下面使用雅可比矩阵对方程组进行求解。?计算雅可比矩阵的函数j( )和f( )一样，其x参数是未知数的一组值，它计算非线性方程组在x处的雅可比矩阵。?通过fprime参数将j( )传递给fsolve( )。由于本例中的未知数很少，因此计算雅可比矩阵并不能显著地提高计算速度。

def j(x): ?
x0, x1, x2 = x.tolist( )
return [
[0, 5, 0],
[8*x0, -2*x2*cos(x1*x2), -2*x1*cos(x1*x2)],
[0, x2, x1]
]

result = optimize.fsolve(f, [1,1,1], fprime=j) ?
print result
print f(result)
[-0.70622057 -0.6 -2.5 ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]

二、最小二乘拟合

假设有一组实验数据(xi,yi)，我们事先知道它们之间应该满足某函数关系：yi=f(xi)。通过这些已知信息，需要确定函数f( )的一些参数。例如，如果函数f( )是线性函数f(x)=kx+b，那么参数k和b就是需要确定的值。

如果用p表示函数中需要确定的参数，则目标是找到一组p使得函数S的值最小：

这种算法被称为最小二乘拟合(least-square fitting)。在optimize模块中，可以使用leastsq( )对数据进行最小二乘拟合计算。leastsq( )的用法很简单，只需要将计算误差的函数和待确定参数的初始值传递给它即可。下面是用leastsq( )对线性函数进行拟合的程序：

import numpy as np
from scipy import optimize

X = np.array([ 8.19, 2.72, 6.39, 8.71, 4.7 , 2.66, 3.78])
Y = np.array([ 7.01, 2.78, 6.47, 6.71, 4.1 , 4.23, 4.05])

def residuals(p): ?
"计算以p为参数的直线和原始数据之间的误差"
k, b = p
return Y - (k*X + b)

# leastsq使得residuals( )的输出数组的平方和最小，参数的初始值为[1,0]
r = optimize.leastsq(residuals, [1, 0]) ?
k, b = r[0]
print "k =",k, "b =",b
k = 0.613495349193 b = 1.79409254326

图1（左）直观地显示了原始数据、拟合直线以及它们之间的误差。?residuals( )的参数p是拟合直线的参数，函数返回的是原始数据和拟合直线之间的误差。图中用数据点到拟合直线在Y轴上的距离表示误差。?leastsq( )使得这些误差的平方和最小，即图中所有正方形的面积之和最小。

由前面的函数S的公式可知，对于直线拟合来说，误差的平方和是直线参数k和b的二次多项式函数，因此可以用如图1（右）所示的曲面直观地显示误差平方和与两个参数之间的关系。图中用红色圆点表示曲面的最小点，它的X-Y轴的坐标就是leastsq( )的拟合结果。

图1　最小化正方形面积之和（左），误差曲面（右）

接下来，让我们再看一个对正弦波数据进行拟合的例子：

def func(x, p): ?
"""
数据拟合所用的函数: A*sin(2*pi*k*x + theta)
"""
A, k, theta = p
return A*np.sin(2*np.pi*k*x+theta)

def residuals(p, y, x): ?
"""
实验数据x, y和拟合函数之间的差，p为拟合需要找到的系数
"""
return y - func(x, p)

x = np.linspace(0, 2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6 # 真实数据的函数参数
y0 = func(x, [A, k, theta]) # 真实数据
# 加入噪声之后的实验数据
np.random.seed(0)
y1 = y0 + 2 * np.random.randn(len(x)) ?

p0 = [7, 0.40, 0] # 第一次猜测的函数拟合参数

# 调用leastsq进行数据拟合
# residuals为计算误差的函数
# p0为拟合参数的初始值
# args为需要拟合的实验数据
plsq = optimize.leastsq(residuals, p0, args=(y1, x)) ?

print u"真实参数:", [A, k, theta]
print u"拟合参数", plsq[0] # 实验数据拟合后的参数

pl.plot(x, y1, "o", label=u"带噪声的实验数据")
pl.plot(x, y0, label=u"真实数据")
pl.plot(x, func(x, plsq[0]), label=u"拟合数据")
pl.legend(loc="best")
真实参数: [10, 0.34, 0.5235987755982988]
拟合参数 [ 10.25218748 0.3423992 0.50817424]

图2显示了带噪声的正弦波拟合。

图2　带噪声的正弦波拟合

程序中，?要拟合的目标函数func( )是一个正弦函数，它的参数p是一个数组，包含决定正弦波的三个参数A、k、theta，分别对应正弦函数的振幅、频率和相角。?待拟合的实验数据是一组包含噪声的数据(x, y1)，其中数组y1为标准正弦波数据y0加上随机噪声。

?用leastsq( )对带噪声的实验数据(x, y1)进行数据拟合，它可以找到数组x和y0之间的正弦关系，即确定A、k、theta等参数。和前面的直线拟合程序不同的是，这里我们将(y1, x)传递给args参数。leastsq( )会将这两个额外的参数传递给residuals( )。?因此residuals( )有三个参数，p是正弦函数的参数，y和x是表示实验数据的数组。

对于这种一维曲线拟合，optimize库还提供了一个curve_fit( )函数，下面使用此函数对正弦波数据进行拟合。它的目标函数与leastsq( )稍有不同，各个待优化参数直接作为函数的参数传入。

def func2(x, A, k, theta):
return A*np.sin(2*np.pi*k*x+theta)

popt, _ = optimize.curve_fit(func2, x, y1, p0=p0)
print popt
[ 10.25218748 0.3423992 0.50817424]

如果频率的初值和真实值的差别较大，拟合结果中的频率参数可能无法收敛于实际的频率。在下面的例子中，由于频率初值的选择不当，导致curve_fit( )未能收敛到真实的参数。这时可以通过其他方法先估算一个频率的近似值，或者使用全局优化算法。在后面的例子中，我们会使用全局优化算法重新对正弦波数据进行拟合。

popt, _ = optimize.curve_fit(func2, x, y1, p0=[10, 1, 0])
print u"真实参数:", [A, k, theta]
print u"拟合参数", popt
真实参数: [10, 0.34, 0.5235987755982988]
拟合参数 [ 0.71093473 1.02074599 -0.1277666 ]

三、计算函数局域最小值

optimize库还提供了许多求函数最小值的算法：Nelder-Mead、Powell、CG、BFGS、Newton-CG、L-BFGS-B等。下面我们用一个实例观察这些优化函数是如何找到函数的最小值的。在本例中，要计算最小值的函数f(x,y)为：

f(x,y)=(1-x)2+100(y-x2)2

这个函数叫作Rosenbrock函数，它经常用来测试最小化算法的收敛速度。它有一个十分平坦的山谷区域，收敛到此山谷区域比较容易，但是在山谷区域搜索到最小点则比较困难。根据函数的计算公式不难看出此函数的最小值是0，在(1,1)处。

为了提高运算速度和精度，有些算法带有一个fprime参数，它是计算目标函数f( )对各个自变量的偏导数的函数。f(x,y)对变量x和y的偏导函数为：

而Newton-CG算法则需要计算海森矩阵，它是一个由自变量为向量的实值函数的二阶偏导数构成的方块矩阵，对于函数f(x1, x2, …, xn)，其海森矩阵的定义如下：

对于本例来说，海森矩阵为一个二阶矩阵：

下面使用各种最小值优化算法计算f(x,y)的最小值，根据其输出可知有些算法需要较长的收敛时间，而有些算法则利用导数信息更快地找到最小点。

def target_function(x, y):
return (1-x)**2 + 100*(y-x**2)**2

class TargetFunction(object):

def __init__(self):
self.f_points = [ ]
self.fprime_points = [ ]
self.fhess_points = [ ]

def f(self, p):
x, y = p.tolist( )
z = target_function(x, y)
self.f_points.append((x, y))
return z

def fprime(self, p):
x, y = p.tolist( )
self.fprime_points.append((x, y))
dx = -2 + 2*x - 400*x*(y - x**2)
dy = 200*y - 200*x**2
return np.array([dx, dy])

def fhess(self, p):
x, y = p.tolist( )
self.fhess_points.append((x, y))
return np.array([[2*(600*x**2 - 200*y + 1), -400*x],
[-400*x, 200]])

def fmin_demo(method):
target = TargetFunction( )
init_point =(-1, -1)
res = optimize.minimize(target.f, init_point,
method=method,
jac=target.fprime,
hess=target.fhess)
return res, [np.array(points) for points in
(target.f_points, target.fprime_points, target.fhess_points)]

methods = ("Nelder-Mead", "Powell", "CG", "BFGS", "Newton-CG", "L-BFGS-B")
for method in methods:
res, (f_points, fprime_points, fhess_points) = fmin_demo(method)
print "{:12s}: min={:12g}, f count={:3d}, fprime count={:3d}, "\
"fhess count={:3d}".format(
method, float(res["fun"]), len(f_points),
len(fprime_points), len(fhess_points))
Nelder-Mead : min= 5.30934e-10, f count=125, fprime count= 0, fhess count= 0
Powell : min= 0, f count= 52, fprime count= 0, fhess count= 0
CG : min= 7.6345e-15, f count= 34, fprime count= 34, fhess count= 0
BFGS : min= 2.31605e-16, f count= 40, fprime count= 40, fhess count= 0
Newton-CG : min= 5.22666e-10, f count= 60, fprime count= 97, fhess count= 38
L-BFGS-B : min= 6.5215e-15, f count= 33, fprime count= 33, fhess count= 0

图3显示了各种优化算法的搜索路径，图中用圆点表示调用f( )时的坐标点，圆点的颜色表示调用顺序；叉点表示调用fprime( )时的坐标点。图中用图像表示二维函数的值，值越大则颜色越浅，值越小则颜色越深。为了更清晰地显示函数的山谷区域，图中显示的实际上是通过对数函数log10( )对f(x,y)进行处理之后的结果。

图3　各种优化算法的搜索路径

四、计算全域最小值

前面介绍的几种最小值优化算法都只能计算局域的最小值，optimize库还提供了几种能进行全局优化的算法，下面以前面的正弦波拟合为例，演示全局优化函数的用法。在使用leastsq( )对正弦波进行拟合时，误差函数residuals( )返回一个数组，表示各个取样点的误差。而函数最小值算法则只能对一个标量值进行最小化，因此最小化的目标函数func_error( )返回所有取样点的误差的平方和。

def func(x, p):
A, k, theta = p
return A*np.sin(2*np.pi*k*x+theta)

def func_error(p, y, x):
return np.sum((y - func(x, p))**2)

x = np.linspace(0, 2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6
y0 = func(x, [A, k, theta])
np.random.seed(0)
y1 = y0 + 2 * np.random.randn(len(x))

我们使用optimize.basinhopping( )全域优化函数找出正弦波的三个参数。它的前两个参数和其他求最小值的函数一样：目标函数和初始值。由于它是全局优化函数，因此初始值的选择并不是太重要。niter参数是全域优化算法的迭代次数，迭代的次数越多，就越有可能找到全域最优解。

在basinhopping( )内部需要调用局域最小值函数，其minimizer_kwargs参数决定了所采用的局域最小值算法以及传递给此函数的参数。下面的程序指定使用L-BFGS-B算法搜索局域最小值，并且将两个对象y1和x传递给该局域最小值求解函数的args参数，而该函数会将这两个参数传递给func_error( )。

result = optimize.basinhopping(func_error, (1, 1, 1),
niter = 10,
minimizer_kwargs={"method":"L-BFGS-B",
"args":(y1, x)})
print result.x
[ 10.25218694 -0.34239909 2.63341582]

虽然频率和相位与原系数不同，但是由于正弦函数的周期性，其拟合曲线是和原始数据重合的，如图4所示。

图4　用basinhopping( )拟合正弦波