Python库之SciPy

工具：Pycharm2019.01，Python3.5

关于ScipyScipy包中有许多工具用来解决科学计算的一些共同问题。不同的子模块负责不同的应用。比如插值、拟合、最优化、图像处理、数值计算以及特殊函数等等。

Scipy中许多模块都是基于numpy的，scipy命名空间中许多主要的函数实际上就是numpy函数比如scipy.cos就是np.cos，故通常一起导入numpy库和scipy库

File input/output: scipy.io

import numpy as np

from scipy import io as siosio.savemat('file_path',{'key',value}) : 将一个用命名和数组的字典保存到matlab样式的.mat文件中。

sio.loadmat('file_path')：导入.mat格式的文件1

2

3

4

5

6

7

8

9

10

11import numpy as np

import matplotlib.pyplot as plt

from scipy import io as sio

a = np.ones((3, 3))

sio.savemat('file.mat',{'a':a}) # 将命名为‘a’的字典保存为.mat格式的文件，默认保存到当前目录下，(.mat)后缀可以不加

data = sio.loadmat('file.mat')

print(data['a'])

- - - - - - - - -

[[1. 1. 1.]

[1. 1. 1.]

[1. 1. 1.]]

Linear algebra operations: scipy.linalg

线形代数运算

import numpy as np

from scipy import linalg1

2

3

4

5

6

7

8

9a=np.array([[1,2],[3,4]])

print(linalg.det(a))  # 求方阵的行列式，要求必须是方阵

print(linalg.inv(a)) # 求逆，必须可逆，否则报错提示奇异(singular不可逆)矩阵

print(np.allclose(np.dot(a,linalg.inv(a)),np.eye(2))) # 判断是否为单位矩阵

- - - - - - - -

-2.0

[[-2. 1. ]

[ 1.5 -0.5]]

True

奇异值分解：

关于奇异值分解：SVD简介1

2

3

4

5import numpy as np

from scipy import linalg

arr = np.arange(9).reshape((3, 3)) + np.diag([1, 0, 1])

print(arr)

uarr, spec, vharr = linalg.svd(arr)  # uarr应该是伪逆

Optimization and fit: scipy.optimize

scipy.optimize提供一种算法来最小化函数值(规模或多维度)，曲线拟合和根查找1

2

3import numpy as np

import matplotlib.pyplot as plt

import scipy.optimize as opt

Curve fitting

曲线拟合，这里以sin()曲线为例1

2

3

4

5

6

7

8

9

10def test_func(x, a, b):  # 定义函数

return a * np.sin(b * x)

x_data = np.linspace(-5, 5, num=50)

y_data = 2.9 * np.sin(1.5 * x_data) + np.random.normal(size=50) # sin函数存在偏差

params, params_covariance =opt.curve_fit(test_func, x_data, y_data, p0=[2, 2]) # 求上面函数参数a，b

plt.scatter(x_data,y_data,label='Data')

plt.legend(loc='upper right')

plt.plot(x_data,test_func(x_data,params[0],params[1]),label='fitted function') # 拟合

plt.legend(loc='upper right')

plt.show() # 显示效果如下

Finding the minimum of a scalar function

寻找标量函数最小值。

例如：

下面函数全局最优(小)值为-1.3，而局部最优值大概3.8。

用opt.minimize()函数，需要定义函数f()，以及传入起点坐标x0，返回找到的最小值信息(包括函数值，坐标等)。

为什么说是局部最优呢，当x0=5时，返回的函数值是8.3，位于x=3.83处。1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18def f(x):

return x**2 + 10*np.sin(x)

x = np.arange(-10, 10, 0.1)

plt.plot(x, f(x))

plt.show()

result=opt.minimize(f,x0=0) # 寻找局部最小值的x坐标

print(result) # 打印返回结果信息，

- - - - - - - - - - -

fun: -7.945823375615215  # 局部最小值

hess_inv: array([[0.08589237]])

jac: array([-1.1920929e-06])

message: 'Optimization terminated successfully.'

nfev: 18

nit: 5

njev: 6

status: 0

success: True

x: array([-1.30644012])  # 局部最小值的坐标

如果我们不知道全局最优点附近信息也就无法合理选择起点X0了，而opt.basinhopping()函数可以找到全局最优解。将上面opt.minimize(f,x0=0)替换成opt.basinhopping(f,x0)即可，x0还是需要赋初值，经测试当x0=7时找到的是全局最优(-1.3,-7.9)，而当x0=5时找到的是局部最优值(3.8,8.3)，故为了找到全局最优保险的方法还是需要多试几个参数。

试试opt.minimize_scalar()方法：能找到全局最优解

opt.minimize_scalar()是只有一个变量时的最小函数值，而opt.minimize()是变量作为一个参数向量传进去的。1

2

3

4

5

6

7def f(x):

return x**2 + 10*np.sin(x)

x = np.arange(-10, 10, 0.1)

result=opt.minimize_scalar(f)

print(result.x, result.fun)

- - - - - - - - - - - - - - - -  - - - - -

-1.3064400120612139 -7.945823375615284

Finding the roots of a scalar function

需求：想找到x=？时，f(x)=0。

方法：opt.root(f,x0)，x0需要赋初值1

2

3

4

5

6

7

8

9

10def f(x):

return x**2 + 10*np.sin(x)

x = np.arange(-10, 10, 0.1)

root=opt.root(f,x0=-3) # 设置起点，只能找到一个，要找到另外一个需要调参使得x0=1

print(root.x,root.fun)

root=opt.root(f,x0=1)

print(root.x,root.fun)

- - - - - - - - - - - - - - - - -

[-2.47948183] [-1.77635684e-15]

[0.] [0.]

将以上三个函数绘制到一张图中：1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43def f(x): # 定义函数

return x**2 + 10*np.sin(x)

x = np.arange(-10, 10, 0.1)

# Global optimization

grid = (-10, 10, 0.1)

xmin_global = opt.brute(f, (grid, )) # 暴力找出参数列表中所得解的最优值

print("Global minima found %s" % xmin_global)

# Constrain optimization

xmin_local = opt.fminbound(f, 0, 10) # 在取区间(0,10)找局部最优解

print("Local minimum found %s" % xmin_local)

root = opt.root(f, 1) # our initial guess is 1

print("First root found %s" % root.x)

root2 = opt.root(f, -2.5)

print("Second root found %s" % root2.x)

fig = plt.figure(figsize=(6, 4)) # 新建画布对象指定size

ax = fig.add_subplot(111) # 将画布分成一行一列指定第一块区域作为对象赋给ax

# Plot the function

ax.plot(x, f(x), 'b-', label="f(x)")

# Plot the minima

xmins = np.array([xmin_global[0], xmin_local])

ax.plot(xmins, f(xmins), 'go', label="Minima")

# Plot the roots

roots = np.array([root.x, root2.x])

ax.plot(roots, f(roots), 'kv', label="Roots")

# Decorate the figure

ax.legend(loc='best') # 自适应选择位置

ax.set_xlabel('x')

ax.set_ylabel('f(x)')

ax.axhline(0, color='gray') # 设置y=0画横线

plt.show()

- - - - - - - - - - - - - - - - - -

Global minima found [-1.30641113]

Local minimum found 3.8374671194983834

First root found [0.]

Second root found [-2.47948183]

未完待续…..