numpy如何求标准差和相关系数（基础计算线性代数）

导读：本文介绍一下在Python科学计算中非常重要的一个库——Numpy。

作者：王天庆

如需转载请联系华章科技

Numpy是Numerical Python extensions 的缩写，字面意思是Python数值计算扩展。Numpy是Python中众多机器学习库的依赖，这些库通过Numpy实现基本的矩阵计算，Python的OpenCV库自然也不例外。

Numpy支持高阶、大量计算的矩阵、向量计算，与此同时提供了较为丰富的函数。Numpy采用友好的BSD许可协议开放源代码。它是一个跨平台的科学计算库，提供了与Matlab相似的功能和操作方法。

虽然科学计算领域一直是Matlab的天下，但是Numpy基于更加现代化的编程语言——Python。而且Python凭借着开源、免费、灵活性、简单易学、工程特性好等特点风靡技术圈，已经成为机器学习、数据分析等领域的主流编程语言。

虽然Matlab提供的包非常多，但是Python因其简单灵活、扩展性强等特点，也诞生了一系列优秀的库。例如Scipy具有大多数Matlab所具备的功能，Matplotlib能够简便地进行数据可视化。虽然当前Matlab的地位仍然难以撼动，但是，随着时间的推移，Python在科学计算上的生态系统也会越来越丰富。

安装Numpy的方法也很简单，使用Python的包管理工具pip或者anaconda便可以实现。例如在shell中输入下列命令行便可以通过pip安装Numpy：

pip install numpy

另外，Numpy是OpenCV的一个依赖库，所以，我们使用pip工具安装好OpenCV库之后，Numpy一般也都已经安装好了。

01 array类型

Numpy的array类型是该库的一个基本数据类型，这个数据类型从字面上看是数组的意思，也就意味着它最关键的属性是元素与维度，我们可以这个数据类型来实现多维数组。

因此，通过这个数据类型，我们可以使用一维数组用来表示向量，二维数组来表示矩阵，以此类推用以表示更高维度的张量。

我们通过下面的例子来简单体会一下在Numpy中array类型的使用。

1. Numpy中array类型的基本使用

import numpy as np

# 引入numpy库，并将其别名为np

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

# 通过np.array()方法创建一个名为array的array类型，参数是一个list

array

# 在shell中输入，返回的结果为：

# array([1, 2, 3, 4])

array.max()

# 获取该array类型中最大的元素值，结果为:

# 4

array.mean()

# 求该array中元素的平均值，结果为:

# 2.5

array.min()

# 获取该array中元素的最小值：

# 1

array * 2

# 直接将该array乘以2，Python将该运算符重载，将每一个元素都乘以了2，其输出结果为：

# array([2, 4, 6, 8])

array 1

# 将每一个元素都加上1，输出结果为：

# array([2, 3, 4, 5])

array / 2

# 将每一个元素都除以2，得到浮点数表示的结果为：

# array([ 0.5, 1. , 1.5, 2. ])

array % 2

# Numpy库除了可以对array实现除法运算，还可以实现取模运算，结果为：

# array([1, 0, 1, 0], dtype=int32)

array.argmax()

# 获取该组数据中元素值最大的一个数据的索引，下标从0开始，其结果为：

# 3

通过上面的代码片段，我们可以了解Numpy中array类型的基本使用方法。我们可以看到，array其实是一个类，通过传入一个list参数来实例化为一个对象，也就实现了对数据的封装。这个对象中包含对各个元素进行计算的基本方法，例如求平均值、求最大值等。除此之外，我们再看一下关于更高维度数据的处理。

2. Numpy对更高维度数据的处理

import numpy as np

array = np.array([[1,2],[3,4],[5,6]])

# 创建一个二维数组，用以表示一个3行2列的矩阵，名为array

array

#在交互式编程界面中输入array，返回结果为：

# array([[1, 2],

# [3, 4],

# [5, 6]])

array.shape

# 查看数据的维度属性，下面的输出结果元组，代表的是3行2列

# (3, 2)

array.size

# 查看array中的元素数量，输出结果为：

# 6

array.argmax()

# 查看元素值最大的元素索引，结果为：

# 5

array.flatten()

# 将shape为(3,2)的array转换为一行表示，输出结果为：

# array([1, 2, 3, 4, 5, 6])

# 我们可以看到，flatten()方法是将多维数据“压平”为一维数组的过程

array.reshape(2,3)

# 将array数据从shape为(3,2)的形式转换为(2,3)的形式：

# array([[1, 2, 3],

# [4, 5, 6]])

除此之外，Numpy还包含了创建特殊类别的array类型的方法，例如。

3. Numpy创建特殊类别的array类型

import numpy as np

array_zeros = np.zeros((2,3,3))

#生成结果为：

# array([[[ 0., 0., 0.],

# [ 0., 0., 0.],

# [ 0., 0., 0.]],

#

# [[ 0., 0., 0.],

# [ 0., 0., 0.],

# [ 0., 0., 0.]]])

array_ones = np.ones((2,3,3))

# 生成所有元素都为1的array，其shape是(2,3,3)

array_ones.shape

# (2, 3, 3)

array_arange = np.arange(10)

# 生成一个array，从0递增到10，步长为1，结果为：

# array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

array_linespace = np.linspace(0,10,5)

# 生成一个array从0到10递增，步长为5，结果为：

# array([ 0. , 2.5, 5. , 7.5, 10. ])

Numpy作为Python的一款著名数值计算库，其在基础计算上的功能也是非常完备的，代码如下。

4. Numpy基础计算演示

import numpy as np

np.abs([1,-2,-3,4])

# 取绝对值，结果为：array([1, 2, 3, 4])

np.sin(np.pi / 2)

# 求余弦值，结果为：1.0

np.arctan(1)

# 求反正切值，结果为：0.78539816339744828

np.exp(2)

# 求自然常数e的2次方，结果为：7.3890560989306504

np.power(2,3)

# 求2的3次方，结果为：8

np.dot([1,2],[3,4])

# 将向量[1,2]与[3,4]求点积，结果为：11

np.sqrt(4)

# 将4开平方，结果为：2.0

np.sum([1,2,3,4])

# 求和，结果为：10

np.mean([1,2,3,4])

# 求平均值，结果为：2.5

np.std([1,2,3,4])

# 求标准差，结果为：1.1180339887498949

除此之外，Numpy所包含的基本计算功能还有很多，例如将array切分、拼接、倒序等。

02 线性代数相关

我们在前面介绍了array类型及其基本操作方法，了解到使用array类型可以表示向量、矩阵和多维张量。线性代数计算在科学计算领域非常重要，在机器学习和数据挖掘领域，线性代数相关函数的使用也是非常频繁的。下面，我们介绍一下Numpy为我们提供的线性代数操作。

5. Numpy提供的线性代数操作

import numpy as np

vector_a = np.array([1,2,3])

vector_b = np.array([2,3,4])

# 定义两个向量vector_a与vector_b

np.dot(vector_a,vector_b)

# 将两个向量相乘，在这里也就是点乘，结果为20

vector_a.dot(vector_b)

# 将vector_a与vector_b相乘，结果为20

np.dot(vector_a,vector_b.T)

'''

将一个行向量与一个列向量叉乘的结果相当于将两个行向量求点积，在这里我们测试了dot()方法。其中array类型的T()方法表示转置。

测试结果表明：

dot()方法对于两个向量默认求其点积。对于符合叉乘格式的矩阵，自动进行叉乘。

我们看一下下面这个例子：

'''

matrix_a = np.array([[1,2],

[3,4]])

# 定义一个2行2列的方阵

matrix_b = np.dot(matrix_a,matrix_a.T)

# 这里将该方阵与其转置叉乘，将结果赋予matrix_b变量

matrix_b

'''

结果为：

array([[ 5, 11],

[11, 25]])

'''

np.linalg.norm([1,2])

# 求一个向量的范数的值，结果为2.2360679774997898

# 如果norm()方法没有指定第二个参数，则默认为求2范数。

np.linalg.norm([1,-2],1)

# 指定第二个参数值为1，即求1范数，我们在前面介绍过，1范数的结果即为向量中各元素绝对值之和，结果为3.0

np.linalg.norm([1,2,3,4],np.inf)

# 求向量的无穷范数，其中np.inf表示正无穷，也就是向量中元素值最大的那个，其结果为4.0

np.linalg.norm([1,2,3,4],-np.inf)

# 同理，求负无穷范数的结果为1，也就是向量中元素的最小值

np.linalg.norm(matrix_b)

# 除了向量可以求范数，矩阵也可以有类似的运算，即为F范数，结果为29.866369046136157

np.linalg.det(matrix_a)

# 求矩阵matrix_a的行列式，结果为-2.0000000000000004

np.trace(matrix_a)

# 求矩阵matrix_a的迹，结果为5

np.linalg.matrix_rank(matrix_a)

# 求矩阵的秩，结果为2

vector_a * vector_b

# 使用*符号将两个向量相乘，是将两个向量中的元素分别相乘，也就是前面我们所讲到的哈达马乘积，结果为array([ 2, 6, 12])

vector_a ** vector_b

# 使用二元运算符**对两个向量进行操作，结果为array([ 1, 8, 81], dtype=int32)

# 表示将向量vector_a中元素对应vector_b中的元素值求幂运算。例如最终结果[1,8,81]可以表示为：

# [1*1,2*2*2,3*3*3*3]

np.linalg.pinv(matrix_a)

'''

求矩阵的逆矩阵，方法pinv()求的是伪逆矩阵，结果为：

array([[-2. , 1. ],

[ 1.5, -0.5]])

不使用伪逆矩阵的算法，直接使用逆矩阵的方法是inv()，即

np.linalg.inv(matrix_a)

结果相同，也为：

array([[-2. , 1. ],

[ 1.5, -0.5]])

'''

03 矩阵的高级函数

我们在前面学习了Numpy的基本数据类型array，同时了解了一些基本的数学运算方法。其实除了前面我们所提到的对矩阵求逆、求秩、求转置等基本运算之外，Numpy还为我们提供了矩阵的分解等更高级的函数。

• 矩阵分解

矩阵分解(decomposition, factorization)是将矩阵拆解为若干个矩阵的相乘的过程。在数值分析，常常被用来实现一些矩阵运算的快速算法，在机器学习领域有非常重要的作用。

例如我们在前面介绍过线性降维的PCA算法，其中就涉及矩阵分解的步骤。今日头条、亚马逊网上商城这类互联网产品，总会根据我们的个人喜好给我们推送一些它认为我们会感兴趣的资讯或商品，这类用于推送消息的系统称为推荐系统(Recommendation System)。

在推荐系统的实现过程中，就用到了矩阵分解算法。例如主流的开源大数据计算引擎Spark在ml机器学习库中通过ALS算法实现了推荐系统，也有的推荐系统采用SVD算法来实现整套系统中的矩阵分解过程。

在Numpy中，为我们提供了基于SVD算法的矩阵分解，SVD算法即为奇异值分解法，相对于矩阵的特征值分解法，它可以对非方阵形式的矩阵进行分解，将一个矩阵A分解为如下形式：

A = U∑V^T

式中，A代表需要被分解的矩阵，设其维度是m×n。U矩阵是被分解为的三个矩阵之一，它是一个m×m的方阵，构成这个矩阵的向量是正交的，被称为左奇异向量；∑是一个m×n的向量，它的特点是除了对角线中的元素外，其余元素都为0。V是一个n×n的方阵，它的转置也是一个方阵，与U矩阵类似，构成这个矩阵的向量也是正交的，被称为右奇异向量。整个奇异值分解算法矩阵的形式如图4-1所示，具体算法实现在此不再赘述。

▲图4-1 SVD算法的矩阵形式

我们使用Numpy演示一下SVD算法的使用。

6. SVD算法演示

import numpy as np

matrix = np.array([

[1,2],

[3,4]])

another_matrix = np.dot(matrix,matrix.T)

# 生成一个矩阵 another_matrix

print(another_matrix)

'''

该矩阵为：

array([[ 5, 11],

[11, 25]])

'''

U,s,V = np.linalg.svd(another_matrix,2)

# 使用奇异值分解法将该矩阵进行分解，分解得到三个子矩阵U,s,V

# 在s矩阵的基础上，生成S矩阵为：

S = np.array([[s[0],0],

[0,s[1]]])

# 我们在下面看一下生成的几个矩阵的样子

print(U)

'''

[[-0.40455358 -0.9145143 ]

[-0.9145143 0.40455358]]

'''

print(s)

'''

[ 29.86606875 0.13393125]

'''

print(V)

'''

[[-0.40455358 -0.9145143 ]

[-0.9145143 0.40455358]]

'''

# 利用生成的U,S,V三个矩阵，我们可以重建回原来的矩阵another_matrix

np.dot(U,np.dot(S,V))

# 输出结果为:

'''

array([[ 5., 11.],

[ 11., 25.]])

'''

在上面的代码片段中，s向量表示的是分解后的∑矩阵中对角线上的元素，所以我们在这里面引入了一个S矩阵，将s向量中的元素放置在这个矩阵中，用以验证分解后的矩阵重建回原先的矩阵A的过程。

仔细的读者可能会注意到，为什么这里使用SVD算法生成的矩阵U与V^T是相同的。大家可能会注意到在上面的代码片段中，为何多了一个生成矩阵another_matrix的过程。这是因为一个矩阵与其转置相乘之后的矩阵是对称矩阵（矩阵中的元素沿着对角线对称），将对称矩阵进行分解后的结果可以表示为：

A = V∑V^T

通过观察上式，我们不难发现U与V矩阵是相同的，因为这个例子中，U与V矩阵本身也是对称矩阵，不论它的转置与否形式都是一样的。

我们在第2章介绍过用于线性降维的PCA算法，该算法中有一个步骤是将协方差矩阵分解然后重建，下面我们演示一下使用Numpy的SVD算法来实现PCA算法的例子：

7. 基于SVD实现PCA算法

import numpy as np

# 零均值化，即中心化，是数据的预处理方法

def zero_centered(data):

matrix_mean = np.mean(data, axis=0)

return data - matrix_mean

def pca_eig(data, n):

new_data = zero_centered(data)

cov_mat = np.dot(new_data.T, new_data) # 也可以用 np.cov() 方法

eig_values, eig_vectors = np.linalg.eig(np.mat(cov_mat)) # 求特征值和特征向量，特征向量是列向量

value_indices = np.argsort(eig_values) # 将特征值从小到大排序

n_vectors = eig_vectors[:, value_indices[-1:-(n 1):-1]] # 最大的n个特征值对应的特征向量

return new_data * n_vectors # 返回低维特征空间的数据

def pca_svd(data, n):

new_data = zero_centered(data)

cov_mat = np.dot(new_data.T, new_data)

U, s, V = np.linalg.svd(cov_mat) # 将协方差矩阵奇异值分解

pc = np.dot(new_data, U) # 返回矩阵的第一个列向量即是降维后的结果

return pc[:, 0]

def unit_test():

data = np.array(

[[2.5, 2.4], [0.5, 0.7], [2.2, 2.9], [1.9, 2.2], [3.1, 3.0], [2.3, 2.7], [2, 1.6], [1, 1.1], [1.5, 1.6],

[1.1, 0.9]])

result_eig = pca_eig(data, 1) # 使用常规的特征值分解法，将2维数据降到1维

print(result_eig)

result_svd = pca_svd(data, 1) # 使用奇异值分解法将协方差矩阵分解，得到降维结果

print(result_svd)

if __name__ == '__main__':

unit_test()

经过降维的数据为：

[-0.82797019 1.77758033 -0.99219749 -0.27421042 -1.67580142 -0.9129491

0.09910944 1.14457216 0.43804614 1.22382056]

我们可以看到，数据已经从2维的变为1维的了，这两个PCA算法的计算结果是相同的。其中pca_eig() 函数是使用常规的特征值分解方法来求解的，读者可以参照前面讲述的PCA算法过程来理解这段代码。pca_svd() 函数使用奇异值分解法来求解的。这段代码虽然相对精简，但是背后是经过复杂的数学推导的，下面简要阐述一下PCA算法中奇异值分解的步骤。

1) PCA算法中得到样本的协方差矩阵是经过零均值化处理的，将其去掉常数部分，则也可表示为：

C = X^TX

其中，X是经过中心化处理后的样本矩阵X. 前面我们介绍过，一个矩阵与其转置矩阵相乘的结果是一个对称矩阵。观察到协方差矩阵C便是一个对称矩阵，那么将其进行奇异值分解后则可以表示为：

C = V∑V^T

2) 将经过中心化的样本矩阵X进行奇异值分解，可以得到：

X = U∑V^T

因此，我们可以得到：

X^TX

= (U∑V^T)^T(U∑V^T)

= V∑^TU^TU∑V^T

= V∑2V^T

奇异矩阵V中的列对应着PCA算法主成分中的主方向，因此可以得到主成分为：

XV = U∑V^TV = U∑

关于更详细的数学推倒过程，读者可参考该网址：

https://stats.stackexchange.com/questions/134282/relationship-between-svd-and-pca-how-to-use-svd-to-perform-pca

• 随机数矩阵

Numpy除了为我们提供常规的数学计算函数和矩阵相关操作之外，还提供了很多功能丰富的模块，随机数模块就是其中一部分。利用随机数模块可以生成随机数矩阵，比Python自带的随机数模块功能要强大，我们看一下下面这个例子。

8. Numpy的随机数功能演示

import numpy as np

# 置随机数种子：

np.random.seed()

# 从[1,3)中生成一个整数的随机数，连续生成10个

np.random.randint(1,3,10)

# 返回：array([2, 2, 1, 1, 1, 1, 2, 1, 2, 2])

# 若要连续产生[1,3)之间的浮点数，可以使用下述方法：

2*np.random.random(10) 1

'''

返回：

array([ 1.25705585, 2.38059578, 1.73232769, 2.12303283, 2.33946996,

2.28020734, 2.15724069, 1.32845829, 2.91361293, 1.78637408])

'''

np.random.uniform(1,3,10)

'''

返回：

array([ 1.37993226, 1.38412227, 1.18063785, 1.75985962, 1.42775752,

1.62100074, 1.71768721, 1.50131522, 2.20297121, 1.08585819])

'''

# 生成一个满足正太分布(高斯分布)的矩阵，其维度是4*4

np.random.normal(size=(4,4))

'''

返回：

array([[-1.81525915, -2.02236963, 0.90969106, 0.25448426],

[-1.04177298, -0.35408201, 1.67850233, -0.70361323],

[-0.30710761, 0.57461312, -0.37867596, -0.74010685],

[-0.94046747, 2.37124816, -0.78503777, -0.33485225]])

'''

# 随机产生10个，n=5,p=0.5的二项分布数据:

np.random.binomial(n=5,p=0.5,size=10)

# 返回：array([2, 0, 1, 3, 3, 1, 3, 3, 4, 2])

data = np.arange(10) # 产生一个0到9的序列

np.random.choice(data,5) # 从data数据中随机采集5个样本，采集过程是有放回的

# 返回：array([0, 0, 1, 6, 2])

np.random.choice(data,5,replace=False)

#从data数据中随机采集5个样本，采集过程是没有放回的

# 返回：array([0, 4, 3, 9, 7])

np.random.permutation(data) # 对data进行乱序，返回乱序结果

# 返回：array([2, 8, 6, 4, 9, 1, 3, 5, 7, 0])

np.random.shuffle(data) # 对data进行乱序，并替换为新的data

print(data)

# 输出：[1 2 8 4 3 6 9 0 5 7]

关于作者：王天庆，长期从事分布式系统、数据科学与工程、人工智能等方面的研究与开发，在人脸识别方面有丰富的实践经验。现就职某世界100强企业的数据实验室，从事数据科学相关技术领域的预研工作。

本文摘编自《Python人脸识别：从入门到工程实践》，经出版方授权发布。

延伸阅读《Python人脸识别》

推荐语：华为资深AI工程师撰写，全面讲解人脸识别各项基础技术、原理和算法，从零实现工程级人脸识别引擎。

,