python scipy库总结

本节导图：https://www.processon.com/view/link/5fd9f5dc5653bb06f344d655

大纲

1. 数值计算库SciPy

SciPy在Numpy的基础上，增加了众多的数学计算、科学计算、以及工程计算中的常用模块，例如：线性代数、方程求解、信号处理、统计学、图像处理、稀疏矩阵等。

接下来我们主要介绍其中的一些重要的子模块，它们在我们后续的NLP生涯中很重要。

2. 拟合与优化

scipy.optimize模块提供了许多数值优化算法，我们可以使用其进行最小二乘拟合和计算函数最小值

最小二乘拟合

最小二乘法是一种数学优化技术，它通过最小化平方误差，来寻求与数据最佳匹配的函数。

如果用p来表示函数中需要确定的参数，目标是找到一组p使得函数S的值最小，于是最小二乘拟合被表示为：

机器学习里线性回归模型的参数估计使用的就是最小二乘法。

线性回归模型，即使用线性模型解决实值预测问题的模型。

这里，通过一个例子，使用optimize.leastsq()函数对二维数据点进行最小二乘拟合，并计算出直线的参数（斜率和截距）。

import numpy as np
from scipy import optimize

# 7个二维数据点，x,y坐标如下
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  # 参数p由斜率k和截距b组成
    return Y - (k * X + b)

# 使用leastsq，输入误差函数和初始参数
r = optimize.leastsq(residuals, [1, 0])
k, b = r[0]
print('斜率k={}, 截距b={}'.format(k, b))

斜率k=0.6134953491930442, 截距b=1.794092543259387

计算函数极值

optimize库还提供了许多求函数最小值的算法，比如：Nelder-Mead、Powell、CG、BFGS、Newton-CG、L-BFGS-B等。

这里使用BFGS算法和CG算法求解Rosenbrock函数的最小值**。

Rosenbrock函数经常被用来测试最小化算法的收敛速度，

• 函数公式如下：

  

• 函数图像如下：

  

此函数的特点是：有一个平坦的山谷区域，收敛到山谷区域较为容易，但在山谷区域搜索到最小点则比较困难。

结合公式和图像，我们易知，最小值发生在(1,1)处，最小值是0

一些优化算法需要指定fprime参数，即给出目标函数的偏导函数。f(x,y)对x, y的偏导函数为：

import numpy as np
from scipy import optimize


def target_func(p):
    """目标函数，rosenbrock函数"""
    x, y = p
    z = (1 - x) ** 2 + 100 * (y - x ** 2) ** 2
    return z


def prime_func(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 = np.array([0, 2])  # 初始点选在(0,2)

# 使用bfgs算法
result1 = optimize.fmin_bfgs(target_func, init_point, prime_func)  # 设置目标函数、初始值和导函数
print(result1)

# 使用cg算法
result2 = optimize.fmin_cg(target_func, init_point, prime_func)
print(result2)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 23
         Function evaluations: 30
         Gradient evaluations: 30
[0.99999989 0.99999977]
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 43
         Gradient evaluations: 43
[1.00000005 1.00000009]

3. 线性代数和矩阵分解

Numpy和SciPy都提供了线性代数库 linalg，但SciPy的线性代数库比Numpy更加全面。

这里，我们主要演示计算矩阵的特征值和特征向量以及矩阵分解

计算矩阵的特征值和特征向量

n*n的矩阵A可以看做一个n维空间的线性变换。如果x为n维空间的一个向量，那Ax（A与x的矩阵乘积）就是对x进行了线性变换，结果是一个新的向量。

如果x为变换矩阵A的特征向量，那么经过Ax的线性变换后，新向量仍与x方向相同，但长度可能发生变换，长度的缩放值由特征值λ决定。

所以，特征向量和特征值有如下公式：

以二维平面上的线性变换矩阵为例。

import numpy as np
from scipy import linalg

# 定义变换矩阵A
A = np.array([[1, -0.3], [-0.1, 0.9]])

# 求解其特征值和特征向量
evalues, evectors = linalg.eig(A)
print(evalues)
print(evectors)

[1.13027756+0.j 0.76972244+0.j]
[[ 0.91724574  0.79325185]
 [-0.3983218   0.60889368]]

奇异值分解

奇异值分解是线性代数中一种重要的矩阵分解技术，在机器学习领域有着重要的应用。

假设X是一个m*n的矩阵，存在一个分解，使得：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Ei96paTy-1637250312971)(C:\Users\winter\AppData\Roaming\Typora\typora-user-images\image-20211118233420758.png)]

其中，矩阵U为m*m维，Σ为m*n维的对角矩阵，矩阵V为n*n维，这样的分解就是奇异值分解，Σ对角线上的元素为矩阵X的奇异值，奇异值按照从大到小排列。

我们来对一个矩阵X进行奇异值分解：

import numpy as np
from scipy import linalg

# 定义矩阵X (3*5)
X = np.arange(15).reshape(3,5)
print(X, X.shape,'\n')

# 返回奇异值
svd_values = linalg.svdvals(X)
print(svd_values, '\n')

# 奇异值分解，返回三个矩阵
U, sigma, V = linalg.svd(X)
print(U, U.shape, '\n')
print(sigma, sigma.shape, '\n')  # sigma是矩阵Σ对角线元素组成的一维矩阵
print(V, V.shape, '\n')

[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]] (3, 5) 

[3.17420265e+01 2.72832424e+00 9.19513106e-16] 

[[-0.15425367 -0.89974393  0.40824829]
 [-0.50248417 -0.28432901 -0.81649658]
 [-0.85071468  0.3310859   0.40824829]] (3, 3) 

[3.17420265e+01 2.72832424e+00 9.19513106e-16] (3,) 

[[-0.34716018 -0.39465093 -0.44214167 -0.48963242 -0.53712316]
 [ 0.69244481  0.37980343  0.06716206 -0.24547932 -0.55812069]
 [ 0.45650723 -0.61153491  0.23629065 -0.46400549  0.38274252]
 [-0.2754905  -0.21468338  0.83080946  0.08439319 -0.42502878]
 [-0.34015605  0.52908988  0.23221189 -0.69106924  0.26992351]] (5, 5) 

​

奇异值的大小衡量一些潜在概念的重要性，我们往往需要筛选出最重要的一些奇异值，来降维重构U Σ V三个矩阵。

keep_dims = 2  # 保留2个维度
sigma2 = np.diag(sigma[:keep_dims])  # 保留2个奇异值，重构sigma矩阵
U2 = U[:, :keep_dims]  # 降维重构矩阵U
V2 = V[:keep_dims, :]  # 降维重构矩阵V
print(U2.shape, sigma2.shape, V2.shape)

(3, 2) (2, 2) (2, 5)

在NLP中，潜在语义分析（LSA）就是使用SVD来进行的。此时，我们的输入矩阵A为词-文档矩阵，每一行代表一个词，每一列代表一个文档，矩阵中的值A(i,j)为词i在文档j中的频数。

如下图示：

我们将 词-文档矩阵A 通过SVD分解，并截取最重要的100个奇异值后，得到三个矩阵：X B Y

• X是降维后的词矩阵，每一行代表一个词，每一列标识了一个潜在的语义概念；
• Y是降维后的文档矩阵，每一列代表一个文档，每一行标识了一个潜在的主题概念；
• B标识每个语义类和每个主题类之间的相关性；

分解后，你如果需要计算 词和词之间、文档和文档 之间的相似度，就可以只计算100维，否则你需要计算的是 50万维 或 100万维。

4. 统计学分布与检验

SciPy中的stats模块包含了多种概率分布的随机变量，这里的随机变量包括：连续随机变量 和 离散随机变量。

from scipy import stats
import numpy as np

正态分布

我们这里以正态分布为例。定义一个符合正态分布的随机变量X

X = stats.norm(loc=1.0, scale=2.0)  # 均值为1，标准差为2
X.stats()  # 打印随机变量

(array(1.), array(4.))#平均值和方差

接下来调用随机变量X的rvs()方法，得到一个多次随机采样值的数组samples

samples = X.rvs(size=1000)  # 采样1000次，得到数组samples
samples.shape

(1000,)

接下来通过np.mean()方法和np.var()方法验证下

np.mean(samples), np.var(samples)

(1.038664111644651, 3.831619958926069)

我们也可以使用fit()方法对随机取样序列samples进行拟合，它返回正态分布的参数（均值和标准差）

stats.norm.fit(samples)

(1.038664111644651, 1.9574524154947086)

卡方分布与卡方检验

卡方分布

若n个相互独立的随机变量ξ₁，ξ₂，…,ξn ，均服从标准正态分布，则这n个服从标准正态分布的随机变量的平方和构成一新的随机变量，其分布规律称为卡方分布。

卡方检验

可以通过卡方分布来进行独立性检验，即：判断观测值与理论值的差异是否只是因为随机误差造成的。独立性检验，通俗点说就是判断是否无关。

我们假设有2个袋子，有5个球，每个袋子中球的分布不同。现在要判断两个袋子中的球是否是平均分布的；

def choose_balls(probs, size):
    """设置当前袋子里球的分布为probs，然后取size次，返回每个球取到的次数"""
    n_ball = len(probs)  # 球的个数
    rv = stats.rv_discrete(values=(range(n_ball), probs))  # 以概率分布probs定义离散随机变量，可取值为球的id号（0 ~ n_ball）
    samples = rv.rvs(size=size)  # 使用随机变量rv进行size次采样
    counts = np.bincount(samples)  # 统计每个球出现的次数
    return counts


counts1 = choose_balls([0.18, 0.24, 0.25, 0.16, 0.17], 500)  # 第1个袋子，5个球以不等概率分布，取500次
counts2 = choose_balls([0.2] * 5, 500)  # 第2个袋子，5个球均匀分布，取500次

np.random.seed(42)
# 将每个球的频数输入chisquare函数
chi1, p1 = stats.chisquare(counts1)  # 卡方检验第1个袋子
print("counts =", counts1, "chi1 =", chi1, "p1 =", p1)

chi2, p2 = stats.chisquare(counts2)  # 卡方检验第2个袋子
print("counts =", counts2, "chi2 =", chi2, "p2 =", p2)

counts = [ 95 106 137  80  82] chi1 = 21.54 p1 = 0.00024741435438580667
counts = [ 99  91 101 111  98] chi2 = 2.08 p2 = 0.7210475511959114

卡方检验的零假设为样本符合目标概率，由上面的检验结果可知：

• 第一个袋子对应的p值只有0.017，也就是说如果第一个袋子中的球真的符合平均分布，那么得到的观测结果[103 109 122 78 88]的概率只有0.017, 因此可以推翻零假设，即袋子中的球不太可能是平均分布的。
• 反之，第二个袋子均匀分布下，得到观测[111 103 106 87 93]的概率为0.428，因此不能推翻零假设，即不能否定袋子中的球是均匀分布的；

我们可以将其概率密度函数画出来：

chil和chi2对应的位置分别为χ^2_1和χ^2_2；p1和p2分别对应χ^2_1右侧的面积和χ^2_2右侧的面积

在机器学习中，我们常常要做特征选择，即删除一些不重要的、与类别无关的特征，我们可以使用卡方检验来选择特征。

我们使用卡方统计度量 词特征和类别 独立性的缺乏程度，卡方越大，独立性越小，相关性越大，特征越重要。以下是邮件二分类（是否垃圾）任务下，特征词**“赌博”**出现的文档数情况：

table = [  # 输入上述表格
    [40, 10],
    [60, 90]
]
chi2, p, _, _ = stats.chi2_contingency(table)  # 卡方检验
print(chi2, p)

22.42666666666667 2.1832165337148537e-06

5. 稀疏矩阵

为什么要用稀疏矩阵

在机器学习建模时，经常会出现许多大型的矩阵，这些矩阵中大部分元素都为0，这就是稀疏矩阵。用数组直接存储稀疏矩阵非常浪费，由于矩阵的稀疏性，我们可以只保存其非零元素，从而节约内存。

在传统的向量空间模型中，我们将每个文档表示成一个固定顺序排列的词向量，向量中的每个值 是当前词在当前文档中的词频。如下图：

这个矩阵是稀疏的，因为它中的很多元素都为0，即词频为0。词库中的词往往有数十万，而文档数也可能成千上万，如果不用稀疏矩阵存储，将会造成巨大的内存浪费，甚至，你的机器会因为内存不够而退出。

常用的几种稀疏矩阵

scipy.sparse中提供了许多稀疏矩阵的格式，每种格式都有不同的用途。其中dok_matrix和lil_matrix适合逐渐添加元素。

from scipy import sparse

先看dok_matrix，

• 它从dict继承而来，采用字典保存矩阵中的非零元素；
• 字典的key为保存(行,列)信息的元组；value为对应的元素值；
• 显然dok_matrix很适合 增删改取 元素；

a = sparse.dok_matrix((10, 5))  # 定义shape为(10,5)的dok稀疏矩阵，所有元素初始为0

a[1, :3] = 1.0, 2.0, 3.0  # 设置第1行，前3列元素

print((list(a.keys())))  # 打印稀疏矩阵的key，即非零元素的索引
print((list(a.values())))  # 打印稀疏矩阵的value，即对应非零元素的值
print(a[0,0], a[1,1])  # 打印(0,0)位置的元素、(1,1)位置的元素

[(1, 0), (1, 1), (1, 2)]
[1.0, 2.0, 3.0]
0.0 2.0

我们再来看lil_matrix

• lil_matrix使用两个列表保存非零元素；
• data保存每行中的非零元素，rows保存非零元素所在的列；
• 这种格式也很适合 增删改取 元素，并且能快速获取行相关的数据；

b = sparse.lil_matrix((10, 5))  # 定义shape为(10,5)的lil稀疏矩阵，初始为0

# 设置非0元素
b[0, 1] = 1
b[1, 0] = 2
b[1, 1] = 3

print(b.data, b.data.shape, '\n')  # 打印data属性，10行，每行为非零元素值的list
print(b.rows, b.data.shape)  # 打印rows属性，10行，每行为非零元素列号的list

[list([1.0]) list([2.0, 3.0]) list([]) list([]) list([]) list([]) list([])
 list([]) list([]) list([])] (10,) 

[list([1]) list([0, 1]) list([]) list([]) list([]) list([]) list([])
 list([]) list([]) list([])] (10,)

最后来看coo_matrix

• coo_matrix使用三个数组row col data来保存非零元素的信息；
• 三个数组等长，row保存非零元素的行，col保存非零元素的列，data保存非零元素的值；
• coo_matrix不支持元素的增删改取，一旦创建后，几乎只能转成其它矩阵去处理；
• 它支持重复元素，即同一位置元素可以出现多次，转化为其它矩阵时，会将同一位置的多个值求和；

# 分别定义4个非零元素的row col data
row = [2, 3, 3, 2]
col = [3, 4, 2, 3]
data = [1, 2, 3, 10]

c = sparse.coo_matrix((data, (row, col)), shape=(5, 6))  # 根据row col data，生成coo稀疏矩阵

# 打印row col data三个属性
print(c.col)
print(c.row)
print(c.data, '\n')

print(c.toarray())  # 将其转化为数组，(2,3)位置的两个元素1,10被求和

[3 4 2 3]
[2 3 3 2]
[ 1  2  3 10] 

[[ 0  0  0  0  0  0]
 [ 0  0  0  0  0  0]
 [ 0  0  0 11  0  0]
 [ 0  0  3  0  2  0]
 [ 0  0  0  0  0  0]]

c[0, 0] = 1  # 尝试修改元素，报错

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

TypeError Traceback (most recent call last)

<ipython-input-35-c3fc3415dffd> in <module>()
----> 1 c[0, 0] = 1  # 尝试修改元素，报错


TypeError: 'coo_matrix' object does not support item assignment

c[0,0]#报错，不能查看