Scipy教程 - 优化和拟合库scipy.optimize

野性酷女 2022-07-26 08:51 527阅读 0赞

http://[blog.csdn.net/pipisorry/article/details/51106570][blog.csdn.net_pipisorry_article_details_51106570]

**最优化函数库Optimization**

优化是找到最小值或等式的数值解的问题。scipy.optimization子模块提供了函数最小值(标量或多维)、曲线拟合和寻找等式的根的有用算法。

from scipy import optimize

[皮皮blog][blog]

# 最小二乘拟合 #

假设有一组实验数据(xi，yi ), 事先知道它们之间应该满足某函数关系yi=f(xi)，通过这些已知信息，需要确定函数f的一些参数。例如，如果函数f是线性函数f(x)=kx+b,那么参数 k和b就是需要确定的值。  
如果用p表示函数中需要确定的参数，那么目标就是找到一组p,使得下面的函数S的值最小：   
![Center][]  
这种算法被称为最小二乘拟合(Least-square fitting)。

在optimize模块中可以使用leastsq()对数据进行最小二乘拟合计算。leastsq() 只需要将计算误差的函数和待确定参数 的初始值传递给它即可。

leastsq()函数传入误差计算函数和初始值，该初始值将作为误差计算函数的第一个参数传入；

计算的结果r是一个包含两个元素的元组，第一个元素是一个数组，表示拟合后的参数k、b；第二个元素如果等于1、2、3、4中的其中一个整数，则拟合成功，否则将会返回mesg。

### 示例 ###

使用最小二乘法对带噪声的正玄波数据进行拟合：

"""  
使用leastsq()对带噪声的正弦波数据进行拟合。拟合所得到的参数虽然和实际的参数有可能完全不同，但是由于正弦函数具有周期性，实际上拟合的结果和实际的函数是一致的。  
"""  
import numpy as np  
from scipy.optimize import leastsq  
  
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(-2\*np.pi, 0, 100)  
A, k, theta = 10, 0.34, np.pi/6 \# 真实数据的函数参数  
y0 = func(x, \[A, k, theta\]) \# 真实数据  
\# 加入噪声之后的实验数据  
y1 = y0 + 2 \* np.random.randn(len(x))   
  
p0 = \[7, 0.2, 0\] \# 第一次猜测的函数拟合参数

\# 调用leastsq进行数据拟合, residuals为计算误差的函数  
\# p0为拟合参数的初始值,\# args为需要拟合的实验数据  
plsq = leastsq(residuals, p0, args=(y1, x))   
\# 除了初始值之外，还调用了args参数，用于指定residuals中使用到的其他参数（直线拟合时直接使用了X,Y的全局变量）,同样也返回一个元组，第一个元素为拟合后的参数数组；

这里将 (y1, x)传递给args参数。Leastsq()会将这两个额外的参数传递给residuals()。因此residuals()有三个参数，p是正弦函数的参数，y和x是表示实验数据的数组。

print u"真实参数:", \[A, k, theta\]   
print u“拟合参数”, plsq\[0\] \# 实验数据拟合后的参数  
  
import pylab as pl  
pl.plot(x, y0, label=u"真实数据")  
pl.plot(x, y1, label=u"带噪声的实验数据")  
pl.plot(x, func(x, plsq\[0\]), label=u"拟合数据")  
pl.legend()  
pl.show()

>>> 真实参数: \[10, 0.34000000000000002, 0.52359877559829882\]  
>>> 拟合参数\[-9.84152775 0.33829767 -2.68899335\]

![Center 1][]

[皮皮blog][blog]

# 标量函数的最小值 #

## Local Optimization ##

最新的最小化函数有这些：

**minimize**(fun, x0\[, args, method, jac, hess, ...\])Minimization of scalar function of one or more variables.

Note: 参数bounds=((0,None), (0,None))代表x的每个维度对应的界限，如果只有一维，也要写成bounds= ((0, None), )，或者使用下面的minimize\_scalar函数。

bounds= ((0, None), (0, None))对应x的初始值x0就应该是二维的，如(1, 0).

example见\[ [scipy.optimize.minimize][]\]

**minimize\_scalar**(fun\[, bracket, bounds, ...\])Minimization of scalar function of one variable.

Note: minimize\_scalar等这些函数返回的是一个OptimizeResult对象，print()结果如：

fun: -2.3627208257312922  
    nfev: 9  
     nit: 8  
 success: True  
       x: 1.333333333395553

可以使用同字典取key来得到最小值x,如OptimizeResult\['x'\]

小示例：

def test():
        import math
        from scipy import optimize
        func = lambda x: x * math.log(x, 2)
        a = optimize.minimize_scalar(func, bounds=[0.1, 1], method='bounded')
        print(a)

fun: -0.5307378454230427  
 message: 'Solution found.'  
    nfev: 9  
  status: 0  
 success: True  
       x: 0.3678794534141907

trunc = [0, 4]max_x = optimize.minimize_scalar(lambda x: -norm.pdf(x, p[0], p[1]) / norm.pdf(x, q[0], q[1]),bounds=trunc)['x']

\[ [scipy.optimize.minimize\_scalar][scipy.optimize.minimize_scalar]\]

OptimizeResultRepresents the optimization result.  
OptimizeWarning

**The minimize function supports the following methods:**  
minimize(method=’Nelder-Mead’)  
minimize(method=’Powell’)  
minimize(method=’CG’)  
minimize(method=’BFGS’)  
minimize(method=’Newton-CG’)  
minimize(method=’L-BFGS-B’)  
minimize(method=’TNC’)  
minimize(method=’COBYLA’)  
minimize(method=’SLSQP’)  
minimize(method=’dogleg’)  
minimize(method=’trust-ncg’)  
**The minimize\_scalar function supports the following methods:**  
minimize\_scalar(method=’brent’)  
minimize\_scalar(method=’bounded’)  
minimize\_scalar(method=’golden’)

\[[module-scipy.optimize][]\]

官网推荐使用上面的这些函数，而不是下面的函数。\{The specific optimization method interfaces below in this subsection are not recommended for use in new scripts; all of these methods are accessible via a newer, more consistent interface provided by the functions above.\}不过lz将使用下面的函数进行说明使用。

optimize模块提供了许多求函数最小值的算法：fmin、fmin\_powell、 fmin\_cg、fmin\_bfgs 等。不同的“fmin\*()”参数有所不同，为了提高运算速度和精度，有些“fmin()” 带有一个fprime参数，它是计算目标函数f对各个自变量的偏导数的函数。

### 求函数最大/最小值的示例1 ###

首先定义一个在区间(p\[0\], p\[1\])上要求最大值点x的函数norm.pdf(x, p\[0\], p\[1\])

Note: 要求最大值前面加一个-号就可以了，好像只有求最小值点的函数

加上求值区间

optimize.fminbound(lambda x: -norm.pdf(x, p[0], p[1]) / norm.pdf(x, q[0], q[1]), trunc[0], trunc[1])

这样就可以计算出最小值点x在哪

代入原来的函数中就可以求出最大值y了

\[[Mathematical optimization: finding minima of functions][Mathematical optimization_ finding minima of functions]\]

### 示例2 ###

找到函数f(x,y)的最小值的，并且绘制出f所表示的曲面和寻找最小值时的搜索路径。  
f(x, y) = (1-x)^2+100(y-x^2)^2  
f(x,y)对变量x和y的偏导函数为：  
![Center 2][]  
这个函数叫做Rosenbrock(罗森布罗克)函数，它经常用来测试最小化算法的收敛速度。它有一个十分平坦的山谷区域，收敛到此山谷区域比较容易，但是在山谷区域搜索到最小点则比较困难。根据函数的计算公式不难看出此函数的最小值是0,在(1,1)处。

"""  
使用fmin()计算函数最小值，并用matplotlib绘制搜索最小值的路径。  
"“”  
import scipy.optimize as opt   
import numpy as np   
import sys  
  
points = \[\]  
def f(p):   
    x, y = p  
    z = (1-x)\*\*2 + 100\*(y-x\*\*2)\*\*2  
    points.append((x,y,z))  
    return z  
def fprime(p):   
    x, y = p  
    dx = -2 + 2\*x - 400\*x\*(y - x\*\*2)  
    dy = 200\*y - 200\*x\*\*2  
    return np.array(\[dx, dy\])  
  
init\_point =(-2,-2)   
  
try:  
    method = sys.argv\[1\]  
except:  
    method = "fmin\_bfgs”  
  
fmin\_func = opt.\_\_dict\_\_\[method\]   
if method in \["fmin", "fmin\_powell"\]:  
    result = fmin\_func(f, init\_point) \#参数为目标函数和初始值  
elif method in \["fmin\_cg", "fmin\_bfgs", "fmin\_l\_bfgs\_b", "fmin\_tnc"\]:  
    result = fmin\_func(f, init\_point, fprime) \#参数为目标函数、初始值和导函数（导函数可选？）  
elif method in \["fmin\_cobyla"\]:  
    result = fmin\_func(f, init\_point, \[\])  
else:  
    print "fmin function not found"  
    sys.exit(0)

f()计算f(x,y)的函数值，为了记录下最小化过程中的计算轨迹，在f()中将每个计算过的 点都添加进全局列表points中。fprime()计算f(x,y)对两个自变量在p处的偏导函数的值。最小化的初值设置为(-2,2)，此程序从optimize模块的\_\_dict\_字典中获得由命令行 参数指定的最小值函数。

### 示例3 ###

通过求解卷积的逆运算演示fmin的功能,比较fmin, fmin\_powell, fmin\_cg, fmin\_bfgs。  
    对于一个离散的线性时不变系统h, 如果它的输入是x，那么其输出y可以用x和h的卷积表示：y = x (\*) h。已知系统的输入x和输出y，计算传递函数h；或已知系统的传递函数h和输出y，计算系统的输入x。这种运算称为反卷积运算。fmin计算反卷积，这种方法只能用在很小规模的数列之上，不过用来评价fmin函数的性能还是不错的。  
import scipy.optimize as opt  
import numpy as np  
  
def test\_fmin\_convolve(fminfunc, x, h, y, yn, x0):  
\#x (\*) h = y, (\*)表示卷积,yn为在y的基础上添加一些干扰噪声的结果,x0为求解x的初始值  
def convolve\_func(h):  
\#计算yn - x (\*) h 的power,fmin将通过计算使得此power最小  
return np.sum((yn - **np.convolve(x, h)**)\*\*2)  
  
\# 调用fmin函数，以x0为初始值  
h0 = fminfunc(convolve\_func, x0)  
  
print fminfunc.\_\_name\_\_  
print “---------------------“  
\# 输出x (\*) h0 和y 之间的相对误差  
print "error of y:", np.sum((np.convolve(x, h0)-y)\*\*2)/np.sum(y\*\*2)  
\# 输出h0 和h 之间的相对误差  
print "error of h:", np.sum((h0-h)\*\*2)/np.sum(h\*\*2)  
print  
  
def test\_n(m, n, nscale):  
"""  
随机产生x, h, y, yn, x0等数列，调用各种fmin函数求解b  
m为x的长度, n为h的长度, nscale为干扰的强度  
""“  
x = np.random.rand(m)  
h = np.random.rand(n)  
y = np.convolve(x, h)  
yn = y + np.random.rand(len(y)) \* nscale  
x0 = np.random.rand(n)  
  
test\_fmin\_convolve(opt.fmin, x, h, y, yn, x0)  
test\_fmin\_convolve(opt.fmin\_powell, x, h, y, yn, x0)  
test\_fmin\_convolve(opt.fmin\_cg, x, h, y, yn, x0)  
test\_fmin\_convolve(opt.fmin\_bfgs, x, h, y, yn, x0)  
  
if \_\_name\_\_ == "\_\_main\_\_":  
test\_n(200, 20, 0.1)  
  
程序的输出：  
fmin  
\---------------------  
error of y: 0.00223334655127  
error of h: 0.0883631777326  
  
fmin\_powell  
\---------------------  
error of y: 0.000121416075668  
error of h: 0.000191138253652   
  
fmin\_cg  
\---------------------  
error of y: 0.000121010258418  
error of h: 0.000191002626186   
  
fmin\_bfgs  
\---------------------  
error of y: 0.000121010214197  
error of h: 0.000191001784589   
[皮皮blog][blog]

# 非线性方程组求解 #

optimize模块中的fsolve()可以对非线性方程组进行求解，它的基本调用形式：fsolve(func, x0)   
参数：func是计算方程组误差的函数，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)\]

### 示例 ###

使用folve求非线性方程组的解，方程如下：(scipy\_fsolve.py)  
![Center 3][]  
from scipy.optimize import fsolve  
from math import sin  
  
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 = fsolve(f, \[1,1,1\])   
print result  
print f(result)  
  
输出为：  
\[-0.70622057 -0.6 -2.5 \]  
\[0.0, -9.1260332624187868e-14, 5.3290705182007514e-15\]

**Note**: 由于fsolve函数在调用函数f时，传递的参数为数组，因此如果直接使用数组中的元素计算的话，计算速度将会有所降低，因此这里先用float函数将数组中的元素转换为Python中的标准浮点数，然后调用标准math库中的函数进行运算。

### 使用雅可比行列式求非线性方程组的解 ###

**雅可比矩阵jacobian**  
雅可比矩阵是一阶偏导数以一定方式排列的矩阵，它给出了可微分方程与给定点的最优线性逼近，因此类似于多元函数的导数。例如前面的函数f1,f2,f3和未知数u1,u2,u3的雅可比矩阵如下：  
![Center 4][]  
在对 **方程组进行求解**时，fsolve会自动计算方程组的雅可比矩阵，如果方程组中的未知数很多，而与每个方程有关的未知数较少时，即雅可比矩阵比较稀疏时，传递一个计算雅可比矩阵的函数将能大幅度提高运算速度。在一个模拟计算的程序中需要大量求解近有50个未知数的非线性方程组的解。每个方程平均与6个未知数相关，通过传递雅可比矩阵的计算函数使计算速度提高了4倍。  
计算雅可比矩阵的函数j()和f()一样，其x参数是未知数的一组值，它计算非线性方程组在x处的雅可比矩阵。通过fprime参数将j()传递给fsolve()。  
  
使用雅可比行列式求非线性方程组的解:  
from scipy.optimize import fsolve  
from math import sin,cos  
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 = fsolve(f, \[1,1,1\], fprime=j)   
  
print result  
print f(result)  
Note: 由于本例中的未知数很少， 因此计算雅可比矩阵并不能显著地提高计算速度。

[皮皮blog][blog]

## 其它使用实例 ##

### 绘制温度限 ###

定义函数来描述最小和最大温度。提示：这个函数以一年为周期。提示：包括时间偏移。  
对数据使用这个函数scipy.optimize.curve\_fit()  
绘制结果。是否拟合合理？如果不合理，为什么？  
拟合精度的最大最小温度的时间偏移是否一样？  
f(x,y) = (4 - 2.1x^2 + \\frac\{x^4\}\{3\})x^2 + xy + (4y^2 - 4)y^2  
  
变量应该限制在-2 < x < 2 , -1 < y < 1.  
使用numpy.meshgrid()和plt.imshow来可视地搜寻区域。  
使用scipy.optimize.fmin\_bfgs()或其它多维极小化器。

然后绘制它。

![optimization][]

该函数在大约-1.3有个全局最小值,在3.8有个局部最小值。

### 关于不同最（极）小化方法的讨论 ###

找到这个函数最小值一般而有效的方法是从初始点使用梯度下降法。BFGS算法是做这个的好方法：这个方法一个可能的问题在于，如果函数有局部最小值，算法会因初始点不同找到这些局部最小而不是全局最小.  
如果我们不知道全局最小值的邻近值来选定初始点，我们需要借助于耗费资源些的全局优化。为了找到全局最小点，最简单的算法是蛮力算法，该算法求出给定格点的每个函数值。对于大点的格点，scipy.optimize.brute()变得非常慢。  
scipy.optimize.anneal()提供了使用模拟退火的替代函数。对已知的不同类别全局优化问题存在更有效率的算法，但这已经超出scipy的范围。一些有用全局优化软件包是 [OpenOpt][]、 [IPOPT][]、 [PyGMO][]和 [PyEvolve][]。

为了找到局部最小，我们把变量限制在(0, 10)之间，使用scipy.optimize.fminbound()

**Note**：在高级章节部分[数学优化：找到函数最小值][Link 1]中有关于寻找函数最小值更详细的讨论。

找到标量函数的根

为了寻找根，例如令f(x)=0的点，对以上的用来示例的函数f我们可以使用scipy.optimize.fsolve()

注意仅仅一个根被找到。检查f的图像在-2.5附近有第二个根。我们可以通过调整我们的初始猜测找到这一确切值：

曲线拟合

假设我们有从被噪声污染的f中抽样到的数据：

如果我们知道函数形式(当前情况是x^2 + sin(x))，但是不知道幅度。我们可以通过最小二乘拟合拟合来找到幅度。首先我们定义一个用来拟合的函数：

然后我们可以使用scipy.optimize.curve\_fit()来找到a和b：

现在我们找到了f的最小值和根并且对它使用了曲线拟合。我们将一切放在一个单独的图像中:

![function][]

注意：**Scipy>=0.11中提供所有最小化和根寻找算法的统一接口**scipy.optimize.minimize(),scipy.optimize.minimize\_scalar()和scipy.optimize.root()。它们允许通过method关键字方便地比较不同算法。

你可以在scipy.optimize中找到用来解决多维问题的相同功能的算法。

练习：曲线拟合温度数据

在阿拉斯加每个月的温度上下限，从一月开始，以摄氏单位给出。

练习：2维最小化

[source code][]

\[[非线性最小二乘拟合：在点抽取地形激光雷达数据上的应用][Link 2]\]

\[[Scipy:高端科学计算][Scipy]\]

from: [http://blog.csdn.net/pipisorry/article/details/51106570][blog.csdn.net_pipisorry_article_details_51106570]

ref: Optimization ([scipy.optimize][module-scipy.optimize])

[tutorial/optimize][tutorial_optimize]

[blog.csdn.net_pipisorry_article_details_51106570]: http://blog.csdn.net/pipisorry/article/details/51106570
[blog]: http://blog.csdn.net/pipisorry
[Center]: /images/20220724/88bfc3b3df554356afb489295daa23c0.png
[Center 1]: /images/20220724/0b3139268cbc430cbbceb4f8f4bca96e.png
[scipy.optimize.minimize]: http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize
[scipy.optimize.minimize_scalar]: http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize_scalar.html
[module-scipy.optimize]: http://docs.scipy.org/doc/scipy/reference/optimize.html#module-scipy.optimize
[Mathematical optimization_ finding minima of functions]: http://www.scipy-lectures.org/advanced/mathematical_optimization/index.html#mathematical-optimization
[Center 2]: /images/20220724/e0a631e250f44218be5ef6d837a0a891.png
[Center 3]: /images/20220724/4d50c23fef5e4011beb6d33372e27dc4.png
[Center 4]: /images/20220724/1f904e17fa2b4ed6a08e125e4dcbed9c.png
[optimization]: http://scipy-lectures.github.com/_images/scipy_optimize_example1.png
[OpenOpt]: http://openopt.org/Welcome
[IPOPT]: https://github.com/xuy/pyipopt
[PyGMO]: http://pagmo.sourceforge.net/pygmo/index.html
[PyEvolve]: http://pyevolve.sourceforge.net/
[Link 1]: http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html#mathematical-optimization
[function]: http://scipy-lectures.github.com/_images/scipy_optimize_example2.png
[source code]: http://scipy-lectures.github.com/plot_directive/pyplots/scipy_optimize_sixhump.py
[Link 2]: http://scipy-lectures.github.com/intro/summary-exercises/optimize-fit.html#summary-exercise-optimize
[Scipy]: http://reverland.org/python/2012/10/22/scipy/#scipyoptimize
[tutorial_optimize]: http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html