SciPy 核心库,其基于 NumPy 构建,继承了 NumPy 对多维数据的高效计算能力。
import scipy
scipy.__version__
从某种意义上将,SciPy 可以被看作是 NumPy 的拓展和补充。
常量模块
scipy.constants 模块,该模块下包含了常用的物理和数学常数及单位。
数学中的圆周率和黄金分割常数
from scipy import constants
constants.pi
constants.golden
物理学经常会用到的真空中的光速、普朗克系数等
constants.c, constants.speed_of_light # 两种写法
constants.h, constants.Planck
线性代数
SciPy 中提供的详细而全面的线性代数计算函数。这些函数基本都放置在模块 scipy.linalg 下方。
大致分为:基本求解方法,特征值问题,矩阵分解,矩阵函数,矩阵方程求解,特殊矩阵构造等几个小类。
求给定矩阵的逆,就可以用到 scipy.linalg.inv 函数。传入函数的一般为 NumPy 给出的数组 np.array 或者矩阵 np.matrix 类型。
import numpy as np
from scipy import linalg
linalg.inv(np.matrix([[1, 2], [3, 4]]))
使用 SciPy 提供的 scipy.linalg.svd 函数可以十分方便地完成奇异值分解。例如,我们对一个 5 \times 45×4 的随机矩阵完成奇异值分解。
U, s, Vh = linalg.svd(np.random.randn(5, 4))
U, s, Vh
最终返回酉矩阵 U 和 Vh,以及奇异值 s。
除此之外,scipy.linalg 还包含像最小二乘法求解函数 scipy.linalg.lstsq。现在尝试用其完成一个最小二乘求解过程。
首先,我们给出样本的 xx 和 yy 值。然后假设其符合

分布。
接下来,我们完成

计算,并添加截距项系数 1。
M = x[:, np.newaxis]**[0, 2]
M
然后使用 linalg.lstsq 执行最小二乘法计算,返回的第一组参数即为拟合系数。
p = linalg.lstsq(M, y)[0]
p
我们可以通过绘图来查看最小二乘法得到的参数是否合理,绘制出样本和拟合曲线图。
from matplotlib import pyplot as plt
%matplotlib inline
plt.scatter(x, y)
xx = np.linspace(0, 10, 100)
yy = p[0] + p[1]*xx**2
plt.plot(xx, yy)
插值函数
插值,是数值分析领域中通过已知的、离散的数据点,在范围内推求新数据点的过程或方法。SciPy 提供的 scipy.interpolate 模块下方就包含了大量的数学插值方法。
首先,我们给出一组 x和 y的值。
x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
y = np.array([0, 1, 4, 9, 16, 25, 36, 49, 64, 81])
plt.scatter(x, y)
接下来,我们想要在上方两个点与点之间再插入一个值。怎样才能最好地反映数据的分布趋势呢?这时,就可以用到线性插值的方法。
from scipy import interpolate
xx = np.array([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5]) # 两点之间的点的 x 坐标
f = interpolate.interp1d(x, y) # 使用原样本点建立插值函数
yy = f(xx) # 映射到新样本点
plt.scatter(x, y)
plt.scatter(xx, yy, marker='*')
可以看到,星号点是我们插值的点,而圆点是原样本点,插值的点能准确符合已知离散数据点的趋势。
图像处理
SciPy 集成了大量针对图像处理的函数和方法。一张彩色图片是由 RGB 通道组成,而这实际上就是一个多维数组。所以,SciPy 针对图像的处理的模块 scipy.ndimage,实际上也是针对多维数组的处理过程,你可以完成卷积、滤波,转换等一系列操作。

先使用
scipy.misc 模块中的 face 方法导入一张浣熊的示例图片。scipy.misc 是一个杂项模块,包含了一些无法被准确归类的方法。
from scipy import misc
face = misc.face()
face
face 默认是图片的 RGB 数组,我们可以对其进行可视化还原
plt.imshow(face)
对图片进行高斯模糊处理
from scipy import ndimage
plt.imshow(ndimage.gaussian_filter(face, sigma=5))
针对图像进行旋转变换
plt.imshow(ndimage.rotate(face, 45))
或者对图像执行卷积操作,首先随机定义一个卷积核 k,然后再执行卷积
k = np.random.randn(2, 2, 3)
plt.imshow(ndimage.convolve(face, k))
优化方法
最优化,是应用数学的一个分支,一般我们需要最小化(或者最大化)一个目标函数,而找到的可行解也被称为最优解。最优化的思想在工程应用领域非常常见。举例来讲,机器学习算法都会有一个目标函数,我们也称之为损失函数。而找到损失函数最小值的过程也就是最优化的过程。我们可能会用到最小二乘法,梯度下降法,牛顿法等最优化方法来完成。
SciPy 提供的 scipy.optimize 模块下包含大量写好的最优化方法。例如上面我们用过的 scipy.linalg.lstsq 最小二乘法函数在 scipy.optimize 模块下也有一个很相似的函数 scipy.optimize.least_squares。这个函数可以解决非线性的最小二乘法问题。
scipy.optimize 模块下最常用的函数莫过于 scipy.optimize.minimize,因为只需要指定参数即可以选择大量的最优化方法。例如 scipy.optimize.minimize(method="BFGS") 准牛顿方法。
接下来,我们沿用上面 scipy.linalg.lstsq 演示时同样的数据,使用 scipy.linalg.lstsq 最小二乘法来搞定最小二乘法计算过程。这里会比上面麻烦一些,首先定义拟合函数 func 和残差函数 err_func,实际上我们需要求解残差函数的最小值。
def func(p, x):
w0, w1 = p
f = w0 + w1*x*x
return f
def err_func(p, x, y):
ret = func(p, x) - y
return ret
from scipy.optimize import leastsq
p_init = np.random.randn(2) # 生成 2 个随机数
x = np.array([1, 2.5, 3.5, 4, 5, 7, 8.5])
y = np.array([0.3, 1.1, 1.5, 2.0, 3.2, 6.6, 8.6])
# 使用 Scipy 提供的最小二乘法函数得到最佳拟合参数
parameters = leastsq(err_func, p_init, args=(x, y))
parameters[0]
信号处理
信号处理(英语:Signal processing)是指对信号表示、变换、运算等进行处理的过程,其在计算机科学、药物分析、电子学等学科中应用广泛。几十年来,信号处理在诸如语音与数据通信、生物医学工程、声学、声呐、雷达、地震、石油勘探、仪器仪表、机器人、日用电子产品以及其它很多的这样一些广泛的领域内起着关键的作用。
SciPy 中关于信号处理的相关方法在 scipy.signal 模块中,其又被划分为:卷积,B-样条,滤波,窗口函数,峰值发现,光谱分析等 13 个小类,共计百余种不同的函数和方法。
首先,我们尝试生成一些规则的波形图,例如锯齿波、方形波等。
from scipy import signal
t = np.linspace(-1, 1, 100)
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
axes[0].plot(t, signal.gausspulse(t, fc=5, bw=0.5))
axes[0].set_title("gausspulse")
t *= 5*np.pi
axes[1].plot(t, signal.sawtooth(t))
axes[1].set_title("chirp")
axes[2].plot(t, signal.square(t))
axes[2].set_title("gausspulse")
尝试应用几个滤波函数,首先生成一组包含噪声的波形图
def f(t): return np.sin(np.pi*t) + 0.1*np.cos(7*np.pi*t+0.3) + \
0.2 * np.cos(24*np.pi*t) + 0.3*np.cos(12*np.pi*t+0.5)
t = np.linspace(0, 4, 400)
plt.plot(t, f(t))
然后,我们使用 SciPy 中的中值滤波函数 scipy.signal.medfilt 和维纳滤波函数 scipy.signal.wiener。
plt.plot(t, f(t))
plt.plot(t, signal.medfilt(f(t), kernel_size=55), linewidth=2, label="medfilt")
plt.plot(t, signal.wiener(f(t), mysize=55), linewidth=2, label="wiener")
plt.legend()
统计函数
scipy.stats 模块包含大量概率分布函数,主要有连续分布、离散分布以及多变量分布。除此之外还有摘要统计、频率统计、转换和测试等多个小分类。基本涵盖了统计应用的方方面面。
比较有代表性的 scipy.stats.norm 正态分布连续随机变量函数为代表进行介绍。我们尝试使用 .rvs 方法随机抽取 1000 个正态分布样本,并绘制出条形图。
from scipy.stats import norm
plt.hist(norm.rvs(size=1000))
除了 norm.rvs 返回随机变量,norm.pdf 返回概率密度函数,norm.cdf 返回累计分布函数,norm.sf 返回残存函数,norm.ppf 返回分位点函数,norm.isf 返回逆残存函数,norm.stats 返回均值,方差,(费舍尔)偏态,(费舍尔)峰度,以及 norm.moment 返回分布的非中心矩。对于这些提到的方法,对于大多数 SciPy 连续变量分布函数(柯西分布,泊松分布等)都是通用的。
可以基于概率密度函数绘制出正态分布曲线
x = np.linspace(norm.ppf(0.01), norm.ppf(0.99), 1000)
plt.plot(x, norm.pdf(x))
返回数据的摘要
from scipy.stats import describe
describe(x)
describe 可以返回给定数组最大值、最小值、平均值、方差等信息。scipy.stats 模块所涉及专业性较强,应用门槛比较高,需要你对统计学理论比较熟悉,才能知晓其中函数的释义。
稀疏矩阵
数值分析中,元素大部分为零的矩阵被称为 * 稀疏矩阵*。反之,如果大部分元素都非零,则这个矩阵是稠密的。在科学与工程领域中求解线性模型时经常出现大型的稀疏矩阵。由于其自身的稀疏特性,通过压缩可以大大节省稀疏矩阵的内存代价。更为重要的是,由于过大的尺寸,标准的算法经常无法操作这些稀疏矩阵。
SciPy 中的 scipy.sparse 模块提供了关于稀疏矩阵的储存方法,而 scipy.sparse.linalg 中由包含了专门针对稀疏线性代数的处理方法。此外,scipy.sparse.csgraph 还包含了一些稀疏矩阵的拓扑图理论。
scipy.sparse 中大致有七类稀疏矩阵储存方法,分别为:csc_matrix 压缩稀疏列矩阵,csr_matrix 压缩稀疏行矩阵,bsr_matrix 块稀疏行矩阵,lil_matrix 基于行的链表稀疏矩阵,dok_matrix 基于字典稀疏矩阵,coo_matrix 坐标格式的稀疏矩阵,dia_matrix 对角线稀疏矩阵。
使用一个简单的 NumPy 数组作为示例,将其转换为稀疏矩阵储存。首先是 csr_matrix 压缩稀疏行矩阵。
from scipy.sparse import csr_matrix
array = np.array([[2, 0, 0, 3, 0, 0], [1, 0, 1, 0, 0, 2], [0, 0, 1, 2, 0, 0]])
csr = csr_matrix(array)
print(csr)
稀疏矩阵保存时,实际上是通过坐标的形式表示出矩阵中非零元素的位置。下面尝试使用 csc_matrix 压缩稀疏列矩阵储存同一个数组。
from scipy.sparse import csc_matrix
csc = csc_matrix(array)
print(csc)
实际上,对比就可以看出二者储存的顺序不一样,一个是按行顺序一个是按列顺序。稀疏矩阵表示,可以通过 todense() 变为稠密矩阵表示。
csc.todense()
通过对比稀疏矩阵和稠密矩阵储存时所耗费的内存大小。首先,我们生成一个 1000 \times 10001000×1000 的随机矩阵,然后将其中小于 1 的值替换为 0 人为变成稀疏矩阵。
from scipy.stats import uniform
data = uniform.rvs(size=1000000, loc=0, scale=2).reshape(1000, 1000)
data[data < 1] = 0
data
接下来,我们打印该稠密矩阵在内存中所耗费的存储大小,单位为 MB。
data.nbytes/(1024**2)
接下来,我们以稀疏矩阵存储,并以同样的单位查看大小。
data_csr = csr_matrix(data)
data_csr.data.size/(1024**2)
不出意外的话,以稀疏矩阵格式存储所耗费的内存空间应该远小于以稠密矩阵存储的方式。