第五讲数组 DONE

1. 数组

1 数组

数组由 Python 的 NumPy 模块定义，含义为“numerical Python”，即 Python 数值计算工具。NumPy 起源于使用 Python 调用 Fortran 进行线性代数计算。历史上 Fortran 是最早的计算机高级语言，是编写数值计算程序的首选。几十年间，Fortran 积累了大量优质的数值计算工具库。Python 被用于科学计算时，最重要的是能无缝调用已有的 Fortran 程序，在巨人的肩膀上前进。 NumPy 从 Fortran 的调用接口开始发展，逐步研发更高级的功能，成为了 Python 科学计算的基础，是 Python 上数值计算的“最佳工具”。所有的 Python 科学计算工具库都沿用了 NumPy 的数据结构定义。 NumPy 还不是 Python 的标准库，需要使用

apt install python3-numpy

来安装。

1.1 创建数组

import numpy as np

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

[1 2 3 4 3 2 1]

np 是约定的 numpy 缩写。第一次用 np.array ，可以用 help 来查看文档学习用法。给 np.array 放进列表，返一个创建数组赋予 nv 。数组的内容确认与列表一样。

1.2 列表与数组的区别

自然的疑问是：既然数组与列表的内容相同，为何还要增加新的数据类型？数组要求其元素的数据类型一致，如果给了不一致的元素，元素会退化成 object 类型，失去大多数运算功能。

print(nv.dtype)

print(np.array([1, "a", None]).dtype)

int64
object

数据类型一致的限制，换取的是数组所占空间可预测的好处。这使得它可以在计算机内存中连续存储，不仅有更高的读写效率，还可以表达更高维的结构。列表由动态链表实现，灵活但是损失了效率。

1.3 数组的索引

数组的索引与列表有同样的基础语法，例如：

nv[2], nv[5:], nv[-1], nv[::2], nv[::-3]

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

::2 代表每两个元素取一次， ::-3 代表从后向前每三个元素取一次。

np.arange 可以直接生连续数字的数组。

a = np.arange(10)
print(a)
print(a[::-1])

[0 1 2 3 4 5 6 7 8 9]
[9 8 7 6 5 4 3 2 1 0]

1.4 二维数组与矩阵

二维数组，先在第一个维度（比如行）排列元素，再在第二个维度排列，与一维数组无本质区别。这样的方法能直接推广到 N 维数组，按维度先后排列元素。二维数组可与矩阵等同，我们来构造一个单位矩阵。

identity = np.array([[1, 0], [0, 1]])
print(identity)

[[1 0]
 [0 1]]

嵌套的列表，先写第一行，再写第二行，可以排成需要的数组。访问二维数组的元素时，下标先写第一个维度，再写第二个。或者先写第一个下标，获得它所指定的一维数组后，再取下标。

print(identity[0, 1])
print(identity[0, 1] is identity[0][1], identity[0, 1] == identity[0][1])

0
False True

前一种方法是从二维数组中取元素，后一种是先生成一个复制了的一维数组，再取其中的元素，因此 is 判断给出 False ，但它们相等。

矩阵的类型是 numpy.ndarray

type(identity)

<class 'numpy.ndarray'>

任意数组都是这个类型。

1.5 数组的生成

用列表把数组元素一个一个写出来较繁琐，多维数组更甚。 NumPy 提供一些生成数组的函数，最简单的生成全1和全0数组。

print(np.ones((3, 3)))
print(np.zeros((3, 4)))

[[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]
[[0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]

单位矩阵用 np.eye 生成， help(np.eye) 可查更多的参数，例如

print(np.eye(4))
print(np.eye(4, k=2))

[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]
[[0. 0. 1. 0.]
 [0. 0. 0. 1.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]

1.6 索引综合练习

熟练使用索引，是发挥数组强大功能的基础。生成一个 (10, 10) 的矩阵。

square = np.arange(100)
square.shape = (10, 10)
print(square)

[[ 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 43 44 45 46 47 48 49]
 [50 51 52 53 54 55 56 57 58 59]
 [60 61 62 63 64 65 66 67 68 69]
 [70 71 72 73 74 75 76 77 78 79]
 [80 81 82 83 84 85 86 87 88 89]
 [90 91 92 93 94 95 96 97 98 99]]

我首先生成了一个长度为 100 的一维数组，随后在保持数据不变的前提下，把它的形状改成了 (10, 10) ，即把它解读成二维方阵。这个操作也可以调用 reshape 函数实现。

np.arange(100).reshape((10, 10))

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, 24, 25, 26, 27, 28, 29],
       [30, 31, 32, 33, 34, 35, 36, 37, 38, 39],
       [40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
       [50, 51, 52, 53, 54, 55, 56, 57, 58, 59],
       [60, 61, 62, 63, 64, 65, 66, 67, 68, 69],
       [70, 71, 72, 73, 74, 75, 76, 77, 78, 79],
       [80, 81, 82, 83, 84, 85, 86, 87, 88, 89],
       [90, 91, 92, 93, 94, 95, 96, 97, 98, 99]])

这个操作有一处反直觉的地方：当我们在数学上定义向量、矩阵、张量时，它的维度已经固定了。但是计算机不论存储多少维的数组，本质都是线性一维的，8GB 的内存就是有从0到80亿的位置线性排列。高维数组不过是把多个标号整齐地一一映射到一维空间而已，是一维数组的另一种形式的表征。例如 (10, 10) 的二维数组，第一个指标增加1时对应内存中的地址加10，第二个指标加1时对应的地址加1。这体现了一种重要原理性的构造思想：简单的数据结构，配合不同的描述，衍生出丰富的形式。这种方法可以构造多阶的张量。

取数组的第0列的所有行，

square[:, 0]

array([ 0, 10, 20, 30, 40, 50, 60, 70, 80, 90])

取标号为3的倍数的行和5的倍数的列，

square[::3, ::5]

array([[ 0,  5],
       [30, 35],
       [60, 65],
       [90, 95]])

倒过来取列

square[::3, ::-1]

array([[ 9,  8,  7,  6,  5,  4,  3,  2,  1,  0],
       [39, 38, 37, 36, 35, 34, 33, 32, 31, 30],
       [69, 68, 67, 66, 65, 64, 63, 62, 61, 60],
       [99, 98, 97, 96, 95, 94, 93, 92, 91, 90]])

从第2行开始，每3行取一次，同时取第3到5列（左闭右开），

square[2::3, 3:5]

array([[23, 24],
       [53, 54],
       [83, 84]])

1.7 数组运算

数组可进行各类运算，如取平方，

square ** 2

array([[   0,    1,    4,    9,   16,   25,   36,   49,   64,   81],
       [ 100,  121,  144,  169,  196,  225,  256,  289,  324,  361],
       [ 400,  441,  484,  529,  576,  625,  676,  729,  784,  841],
       [ 900,  961, 1024, 1089, 1156, 1225, 1296, 1369, 1444, 1521],
       [1600, 1681, 1764, 1849, 1936, 2025, 2116, 2209, 2304, 2401],
       [2500, 2601, 2704, 2809, 2916, 3025, 3136, 3249, 3364, 3481],
       [3600, 3721, 3844, 3969, 4096, 4225, 4356, 4489, 4624, 4761],
       [4900, 5041, 5184, 5329, 5476, 5625, 5776, 5929, 6084, 6241],
       [6400, 6561, 6724, 6889, 7056, 7225, 7396, 7569, 7744, 7921],
       [8100, 8281, 8464, 8649, 8836, 9025, 9216, 9409, 9604, 9801]])

是把数组的每个元素都平方了。对比二重循环的实现，

for r in square:
    for c in r:
	print(c**2, end=" ")
    print()

0 1 4 9 16 25 36 49 64 81 
100 121 144 169 196 225 256 289 324 361 
400 441 484 529 576 625 676 729 784 841 
900 961 1024 1089 1156 1225 1296 1369 1444 1521 
1600 1681 1764 1849 1936 2025 2116 2209 2304 2401 
2500 2601 2704 2809 2916 3025 3136 3249 3364 3481 
3600 3721 3844 3969 4096 4225 4356 4489 4624 4761 
4900 5041 5184 5329 5476 5625 5776 5929 6084 6241 
6400 6561 6724 6889 7056 7225 7396 7569 7744 7921 
8100 8281 8464 8649 8836 9025 9216 9409 9604 9801

数组的表达非常直观，而且更本质。

对数组的元素做加法

square[::-1, ::-1] + square

array([[99, 99, 99, 99, 99, 99, 99, 99, 99, 99],
       [99, 99, 99, 99, 99, 99, 99, 99, 99, 99],
       [99, 99, 99, 99, 99, 99, 99, 99, 99, 99],
       [99, 99, 99, 99, 99, 99, 99, 99, 99, 99],
       [99, 99, 99, 99, 99, 99, 99, 99, 99, 99],
       [99, 99, 99, 99, 99, 99, 99, 99, 99, 99],
       [99, 99, 99, 99, 99, 99, 99, 99, 99, 99],
       [99, 99, 99, 99, 99, 99, 99, 99, 99, 99],
       [99, 99, 99, 99, 99, 99, 99, 99, 99, 99],
       [99, 99, 99, 99, 99, 99, 99, 99, 99, 99]])

结合索引的倒取，非常直观地构造出了较复杂的运算，对多个数字操作，形式上与一个数一样。NumPy 的便利性，使用得当，容易表达出简洁优美的运算逻辑。初学者值得适合在线文档多多练习。

相比于按元素运算，总结性运算是把数组化成低维的，极端情形是化成0维，即数字。例如，取平均、中位数与总和，

np.mean(square), np.median(square), np.sum(square)

(49.5, 49.5, 4950)

可选只针对一个维度，使用 axis 参数。

np.sum(square, axis=0), np.mean(square, axis=1)

(array([450, 460, 470, 480, 490, 500, 510, 520, 530, 540]),
 array([ 4.5, 14.5, 24.5, 34.5, 44.5, 54.5, 64.5, 74.5, 84.5, 94.5]))

同样方法推广到3维数组，3阶张量，

cube = np.arange(64).reshape((4, 4, 4))
np.sum(cube, axis=0)

array([[ 96, 100, 104, 108],
       [112, 116, 120, 124],
       [128, 132, 136, 140],
       [144, 148, 152, 156]])

对第0维求和，就还剩两个维度，组成矩阵。也可以对两个维度操作

np.sum(cube, axis=(0, 1))

array([480, 496, 512, 528])

1.8 Pauli 矩阵运算练习

定义三个 Pauli 矩阵，并放到列表里。

pauli = []
pauli.append(np.array([0,1,1,0]).reshape(2,2))
pauli.append(np.array([0,-1j,1j,0]).reshape(2,2))
pauli.append(np.array([1,0,0,-1]).reshape(2,2))
for m in pauli:
    print(m)

[[0 1]
 [1 0]]
[[ 0.+0.j -0.-1.j]
 [ 0.+1.j  0.+0.j]]
[[ 1  0]
 [ 0 -1]]

1j 是由 Python 定义的虚数单位。

Pauli 矩阵的平方都是单位阵。

for m in pauli:
    print(np.dot(m, m))

[[1 0]
 [0 1]]
[[1.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]
[[1 0]
 [0 1]]

它们的对易关系是 Pauli 矩阵的核心性质。

def commute(a, b):
    '''
    给出 a, b 的对易子 [a, b]:=ab - ba
    '''
    return a@b - b@a
for i in range(3):
    l = (i+1) % 3
    m = (i+2) % 3
    if np.all(commute(pauli[i], pauli[l]) == 2j * pauli[m]):
	print(f"[ pauli_{i} , pauli_{l} ] == 2i pauli_{m}")

[ pauli_0 , pauli_1 ] == 2i pauli_2
[ pauli_1 , pauli_2 ] == 2i pauli_0
[ pauli_2 , pauli_0 ] == 2i pauli_1

其中 np.all 仅当数组所有元素都为 True 时返回 True 。

1.8.1 张量运算

\( \sigma_i \) 看起来是三个元素的向量，但算上 pauli 的方阵，本质上是 3, 2, 2 的张量。如果我们全盘使用张量，有可能把所有的 for 循环去掉。\( \epsilon_{ijk} \) 是 3, 3, 3 的全反称张量，对所有的指标奇排列值为 -1 ，所有的指标偶排列值为 1。

# 把 pauli 矩阵构造成 3, 2, 2 张量
p_tensor = np.array(pauli)

# 全反称张量尝试

# epsilon_(0,1,2) = 120 = 201 = 1
# epsilon_(2,1,0) = 102 = 021 = -1
# 其它 = 0

eps = np.zeros((3,3,3))
eps[0,1,2] = eps[1,2,0] = eps[2,0,1] = 1
eps[0,2,1] = eps[1,0,2] = eps[2,1,0] = -1

# 传进来的是 (3, 1, 2, 2) 和 (1, 3, 2, 2) 的张量，得到 (3,3,2,2)
commuted = commute(p_tensor[:, None, :, :] , p_tensor[None, :, :, :])
if np.all(commuted == 2j * np.tensordot(eps, p_tensor, axes=1)):
    print("Pauli 矩阵对易关系验证成功")

Pauli 矩阵对易关系验证成功

1.8.2 特征值和迹

Pauli 阵的两个特征值分别是 \(\pm 1\)，要使用 NumPy 的 linalg （意为 linear algebra）子模块的特征值函数。

for m in pauli:
    print(f"本征值是{np.linalg.eigvals(m)}，迹是{np.trace(m)}")

本征值是[ 1. -1.]，迹是0
本征值是[ 1.+0.j -1.+0.j]，迹是0j
本征值是[ 1. -1.]，迹是0

使用高阶张量，省去迹中的 for 循环。

np.trace(p_tensor, axis1 = 1, axis2 = 2)

array([0.+0.j, 0.+0.j, 0.+0.j])

1.9 Python 工具的学习策略

NumPy 的工具非常丰富，建议大家边学边用，带着问题实现程序。这样学到的东西都会马上应用。用得越多，印象越深。记得住的，都是有用的，不会学到无用的知识。因为工具实在浩如烟海，比如在Python软件库中，不下10万种工具，不可能都了解。只有学会在线调用文档，查阅和甄别学习资料是不变的。

第五讲 数组 DONE

Table of Contents