numpy教程06---ndarray的进阶操作
欢迎关注公众号【Python开发实战】, 获取更多内容!
工具-numpy
numpy是使用Python进行数据科学的基础库。numpy以一个强大的N维数组对象为中心，它还包含有用的线性代数，傅里叶变换和随机数函数。
线性代数
numpy中二维的ndarray可以在Python中高效地表示矩阵，下面将介绍一些主要的矩阵运算。
导入numpy
import numpy as np
矩阵转置
当秩大于等于2时，T属性相当于调用transpose()函数。
m1 = np.arange(10).reshape(2, 5)
m1
输出：
array([[0, 1, 2, 3, 4],
[5, 6, 7, 8, 9]])
m1.T
输出：
array([[0, 5],
[1, 6],
[2, 7],
[3, 8],
[4, 9]])
m1.transpose()
输出：
array([[0, 5],
[1, 6],
[2, 7],
[3, 8],
[4, 9]])
T属性对秩为0（空）或秩为1的ndarray没有影响。
m2 = np.arange(5)
m2
输出：
array([0, 1, 2, 3, 4])
m2.T
输出：
array([0, 1, 2, 3, 4])
可以将一维的ndarray重塑为单行的二维矩阵，进而得到转置。
m2r = m2.reshape(1, 5)
m2r
输出：
array([[0, 1, 2, 3, 4]])
m2r.T
输出：
array([[0],
[1],
[2],
[3],
[4]])
矩阵乘法
dot()方法可以计算两个矩阵的乘法，矩阵乘法需要满足左边矩阵的列数必须等于右边矩阵的行数。
n1 = np.arange(10).reshape(2, 5)
n1
输出：
array([[0, 1, 2, 3, 4],
[5, 6, 7, 8, 9]])
n2 = np.arange(15).reshape(5, 3)
n2
输出：
array([[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8],
[ 9, 10, 11],
[12, 13, 14]])
n1.dot(n2)
输出：
array([[ 90, 100, 110],
[240, 275, 310]])
注意：n1*n2不是矩阵乘法，而是按元素乘积。
矩阵的逆与伪逆
许多线性代数函数在numpy中都可用。linalg模块中的inv函数可以计算一个方阵的逆。
import numpy.linalg as linalg
m3 = np.array([[1, 2, 3], [5, 7, 11], [21, 29, 31]])
m3
输出：
array([[ 1, 2, 3],
[ 5, 7, 11],
[21, 29, 31]])
linalg.inv(m3)
输出：
array([[-2.31818182, 0.56818182, 0.02272727],
[ 1.72727273, -0.72727273, 0.09090909],
[-0.04545455, 0.29545455, -0.06818182]])
也可以通过pinv函数来计算伪逆。
linalg.pinv(m3)
输出：
array([[-2.31818182, 0.56818182, 0.02272727],
[ 1.72727273, -0.72727273, 0.09090909],
[-0.04545455, 0.29545455, -0.06818182]])
单位矩阵
矩阵与其逆矩阵相乘返回一个单位矩阵，下面的例子会有很小的浮点误差。
m3.dot(linalg.inv(m3))
输出：
array([[ 1.00000000e+00, -5.55111512e-17, 0.00000000e+00],
[-2.98372438e-16, 1.00000000e+00, -5.55111512e-17],
[ 5.78009862e-15, 1.27675648e-15, 1.00000000e+00]])
可以通过eye函数来创建一个N×N大小的单位矩阵。
np.eye(3)
输出：
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
QR分解
linalg模块中的qr函数可以计算一个矩阵的QR分解（正交三角分解）。
q, r = linalg.qr(m3)
q
输出：
array([[-0.04627448, 0.98786672, 0.14824986],
[-0.23137241, 0.13377362, -0.96362411],
[-0.97176411, -0.07889213, 0.22237479]])
r
输出：
array([[-21.61018278, -29.89331494, -32.80860727],
[ 0. , 0.62427688, 1.9894538 ],
[ 0. , 0. , -3.26149699]])
q.dot(r)
输出：
array([[ 1., 2., 3.],
[ 5., 7., 11.],
[21., 29., 31.]])
矩阵的行列式
linalg模块中的det函数可以计算矩阵的行列式。
linalg.det(m3)
输出：
43.99999999999999
特征值和特征向量
linalg模块中的eig函数可以计算一个方阵的特征值和特征向量。
eigenvalues, eigenvetors = linalg.eig(m3)
eigenvalues # λ
输出：
array([42.26600592, -0.35798416, -2.90802176])
eigenvetors # v
输出：
array([[-0.08381182, -0.76283526, -0.18913107],
[-0.3075286 , 0.64133975, -0.6853186 ],
[-0.94784057, -0.08225377, 0.70325518]])
m3.dot(eigenvetors) - eigenvalues * eigenvetors # m3.v -λ*v= 0
输出：
array([[ 9.76996262e-15, 2.22044605e-16, -3.10862447e-15],
[ 7.10542736e-15, 2.02615702e-15, -1.11022302e-15],
[ 2.84217094e-14, 5.11049536e-15, -4.88498131e-15]])
奇异值分解
linalg模块中的svd函数可以计算矩阵的奇异值分解。
m4 = np.array([[1, 0, 0, 0, 2], [0, 0, 3, 0, 0], [0, 0, 0, 0, 0], [0, 2, 0, 0, 0]])
m4
输出：
array([[1, 0, 0, 0, 2],
[0, 0, 3, 0, 0],
[0, 0, 0, 0, 0],
[0, 2, 0, 0, 0]])
U, S_diag, V = linalg.svd(m4)
U
输出：
array([[ 0., 1., 0., 0.],
[ 1., 0., 0., 0.],
[ 0., 0., 0., -1.],
[ 0., 0., 1., 0.]])
S_diag # Σ对角线上的值
输出：
array([3. , 2.23606798, 2. , 0. ])
V
输出：
array([[-0. , 0. , 1. , 0. , 0. ],
[ 0.4472136 , 0. , 0. , 0. , 0.89442719],
[-0. , 1. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 1. , 0. ],
[-0.89442719, 0. , 0. , 0. , 0.4472136 ]])
svd函数只返回Σ对角线上的值，完整的矩阵可以这样创建：
S = np.zeros((4, 5))
S[np.diag_indices(4)] = S_diag # np.diag_indices函数返回索引以访问（4，4）数组的主对角线。
S # Σ
输出：
array([[3. , 0. , 0. , 0. , 0. ],
[0. , 2.23606798, 0. , 0. , 0. ],
[0. , 0. , 2. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ]])
U.dot(S).dot(V) # U.Σ.V = m4
输出：
array([[1., 0., 0., 0., 2.],
[0., 0., 3., 0., 0.],
[0., 0., 0., 0., 0.],
[0., 2., 0., 0., 0.]])
对角线和轨迹
m3
输出：
array([[ 1, 2, 3],
[ 5, 7, 11],
[21, 29, 31]])
np.diag(m3) # 返回m3对角线上的元素值
输出：
array([ 1, 7, 31])
np.trace(m3) # 相当于 np.diag(m3).sum()
输出：
39
np.diag(m3).sum()
输出：
39
求解线性标量方程组
linalg模块中的solve函数可以求解线性标量方程组, 例如如下方程组
- 2x + 6y = 6
- 5x + 3y = -9
coeffs = np.array([[2, 6], [5, 3]])
depvars = np.array([6, -9])
solution = linalg.solve(coeffs, depvars)
solution
输出：
array([-3., 2.])
可以验证一下求解：
coeffs.dot(solution), depvars
输出：
(array([ 6., -9.]), array([ 6, -9]))
还可以通过另一个方式验证求解。
np.allclose(coeffs.dot(solution), depvars) # np.allclose比较两个ndarray的每一个元素是否都相等
输出：
True
矢量化
如果坚持进行ndarray操作，而不是一次一个地对单个的元素极性操作，那么代码的效率要高的多， 这就叫做矢量化。这样，可以从numpy的许多优化中受益。
例如，要根据公式sin(xy/40.5)生成一个768×1024的ndarray，一个不好的选择就是使用嵌套循环进行计算。
import math
data = np.empty((768, 1024))
for y in range(768):
for x in range(1024):
data[y, x] = math.sin(x*y/40.5)
data
输出：
array([[0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[0. , 0.02468885, 0.04936265, ..., 0.07705885, 0.1016508 ,
0.12618078],
[0. , 0.04936265, 0.09860494, ..., 0.15365943, 0.20224852,
0.25034449],
...,
[0. , 0.03932283, 0.07858482, ..., 0.6301488 , 0.59912825,
0.56718092],
[0. , 0.06398059, 0.12769901, ..., 0.56844086, 0.51463783,
0.45872596],
[0. , 0.08859936, 0.17650185, ..., 0.50335246, 0.42481591,
0.34293805]])
上面的方法虽然可行，但是效率非常低，因为循环是在纯Python中进行的。现在我们把这个方法矢量化，首先，使用numpy的meshgrid函数，该函数可以根据坐标向量生成坐标矩阵。
x_coords = np.arange(1024)
y_coords = np.arange(768)
X, Y = np.meshgrid(x_coords, y_coords)
X
输出：
array([[ 0, 1, 2, ..., 1021, 1022, 1023],
[ 0, 1, 2, ..., 1021, 1022, 1023],
[ 0, 1, 2, ..., 1021, 1022, 1023],
...,
[ 0, 1, 2, ..., 1021, 1022, 1023],
[ 0, 1, 2, ..., 1021, 1022, 1023],
[ 0, 1, 2, ..., 1021, 1022, 1023]])
Y
输出：
array([[ 0, 0, 0, ..., 0, 0, 0],
[ 1, 1, 1, ..., 1, 1, 1],
[ 2, 2, 2, ..., 2, 2, 2],
...,
[765, 765, 765, ..., 765, 765, 765],
[766, 766, 766, ..., 766, 766, 766],
[767, 767, 767, ..., 767, 767, 767]])
X.shape, Y.shape
输出：
((768, 1024), (768, 1024))
X和Y都是768×1024的ndarray，X中所有的值对应水平轴的坐标， Y中所有的值对应垂直轴的坐标。现在可以简单地使用ndarray运算计算结果。
data = np.sin(X*Y/40.5)
data
输出：
array([[0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[0. , 0.02468885, 0.04936265, ..., 0.07705885, 0.1016508 ,
0.12618078],
[0. , 0.04936265, 0.09860494, ..., 0.15365943, 0.20224852,
0.25034449],
...,
[0. , 0.03932283, 0.07858482, ..., 0.6301488 , 0.59912825,
0.56718092],
[0. , 0.06398059, 0.12769901, ..., 0.56844086, 0.51463783,
0.45872596],
[0. , 0.08859936, 0.17650185, ..., 0.50335246, 0.42481591,
0.34293805]])
保存和加载
numpy可以方便地以二进制或文本格式来保存和加载ndarray。
二进制格式 .npy
下面创建一个随机的ndarray，并保存。
a = np.random.rand(2, 3)
a
输出：
array([[0.76953407, 0.79648012, 0.8868019 ],
[0.88131806, 0.57297333, 0.70059907]])
np.save('my_array', a)
由于文件名不包含扩展名，因此numpy会自动添加扩展名为.npy，下面查看一下文件内容：
with open('my_array.npy', 'rb') as f:
content = f.read()
content
输出：
b"\x93NUMPY\x01\x00v\x00{'descr': '<f8', 'fortran_order': False, 'shape': (2, 3), } \n0L\xbb\xe8\x05\xa0\xe8?\x13s\x00\xdf\xc3|\xe9?\xdd\x05\xf4`\xae`\xec?S\x0e\xad\xee\xc13\xec?\x03\xa6E)\xccU\xe2?\x97\xc6\xdc\xbcNk\xe6?"
想要把该文件加载到numpy数组中只需要调用load函数。
a_loaded = np.load('my_array.npy')
a_loaded
输出：
array([[0.76953407, 0.79648012, 0.8868019 ],
[0.88131806, 0.57297333, 0.70059907]])
文本格式
现在以文本格式来保存ndarray。
np.savetxt('my_array.csv', a)
现在查看一下文件内容：
with open('my_array.csv', 'rt') as f:
print(f.read())
输出：
7.695340676822350900e-01 7.964801173689884939e-01 8.868019002549619723e-01
8.813180600777265061e-01 5.729733282179839682e-01 7.005990685194828371e-01
这是以制表符作为分隔符的csv文件，还可以设置不同的分隔符。
np.savetxt('my_array.csv', a, delimiter=',')
再查看一下文件内容
with open('my_array.csv', 'rt') as f:
print(f.read())
输出：
7.695340676822350900e-01,7.964801173689884939e-01,8.868019002549619723e-01
8.813180600777265061e-01,5.729733282179839682e-01,7.005990685194828371e-01
加载该文件，只需要调用loadtxt函数。
a_loaded = np.loadtxt('my_array.csv', delimiter=',')
a_loaded
输出：
array([[0.76953407, 0.79648012, 0.8868019 ],
[0.88131806, 0.57297333, 0.70059907]])
压缩格式 .npz
还可以在一个压缩文件中保存多个ndarray。
b = np.arange(24).reshape(2, 3, 4)
b
输出：
array([[[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]],
[[12, 13, 14, 15],
[16, 17, 18, 19],
[20, 21, 22, 23]]])
np.savez('my_arrays', my_a=a, my_b=b)
现在看一下文件内容，文件的扩展名.npz已经自动添加。
with open('my_arrays.npz', 'rb') as f:
content = f.read()
repr(content)[:180] + '[...]'
输出：
'b"PK\\x03\\x04\\x14\\x00\\x00\\x00\\x00\\x00\\x00\\x00!\\x00&\\xa9j\\xe4\\xb0\\x00\\x00\\x00\\xb0\\x00\\x00\\x00\\x08\\x00\\x00\\x00my_a.npy\\x93NUMPY\\x01\\x00v\\x00{\'descr\': \'<f8\', \'fortran_order\': False, \'s[...]'
然后可以这样加载这个文件。
my_arrays = np.load('my_arrays.npz')
my_arrays
输出：
<numpy.lib.npyio.NpzFile at 0x93b8080>
这是一个类似字典的对象，它会懒加载ndarray。
my_arrays.keys()
输出：
['my_a', 'my_b']
my_arrays['my_a']
输出：
array([[0.76953407, 0.79648012, 0.8868019 ],
[0.88131806, 0.57297333, 0.70059907]])
my_arrays['my_b']
输出：
array([[[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]],
[[12, 13, 14, 15],
[16, 17, 18, 19],
[20, 21, 22, 23]]])
