NumPy基础：数组创建与操作、矩阵运算、广播机制、随机数生成、性能优化技巧

各位做材料计算的朋友，大家好。今天我们来聊聊NumPy。说实话，我刚开始接触材料计算那会儿，觉得NumPy就是个数学库，没啥了不起。直到有一次处理一个300万原子的分子动力学轨迹文件，用纯Python列表跑了整整两天没跑完，换成NumPy数组后，半小时出结果。嗯，从那以后，我再也不敢小看它了。

NumPy是Python科学计算的基石。做材料计算，不管是处理晶格常数、计算能带结构，还是分析径向分布函数，都离不开它。说白了，它就是给Python装上了「矩阵引擎」，让我们的计算飞起来。

一、数组创建：从零开始搭积木

创建数组是第一步。我个人习惯用 np.array() 从列表直接转，但更常用的是那些「批量生成」的函数。

import numpy as np

# 从列表创建
a = np.array([1, 2, 3, 4, 5])

# 全零数组，常用于初始化
zeros = np.zeros((3, 4))  # 3行4列

# 全一数组
ones = np.ones((2, 3))

# 等差数列，做能带计算时经常用
x = np.linspace(0, 10, 100)  # 从0到10均匀取100个点

# 单位矩阵，计算弹性常数时必备
I = np.eye(4)
我的小经验： 用 np.linspace 代替 np.arange 做等间距采样，能避免浮点数精度带来的边界问题。我在处理布里渊区路径时吃过这个亏，后来就只用 linspace 了。

二、数组操作：切片、变形、拼接

数组操作是日常最频繁的工作。你想想看，一个三维的电荷密度数据，你要提取某个晶面上的值，或者把多个体系的能带数据拼在一起，都离不开这些操作。

# 切片操作，和列表类似
data = np.random.rand(10, 5)
first_row = data[0, :]      # 取第一行
first_col = data[:, 0]      # 取第一列
block = data[2:5, 1:4]      # 取子块

# 变形，reshape 是神器
a = np.arange(12)
b = a.reshape(3, 4)         # 变成3行4列

# 拼接，做数据合并时常用
x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
z = np.concatenate([x, y])  # 水平拼接

# 转置
matrix = np.random.rand(4, 4)
T = matrix.T
注意： reshape 返回的是视图还是副本？这取决于内存布局。我曾经在修改 reshape 后的数组时，意外改变了原数组的数据，排查了半天才发现。建议不确定时用 .copy() 显式复制。

三、矩阵运算：材料计算的核心引擎

矩阵运算在材料计算中无处不在。计算晶格矢量的点积、求解本征值、做线性变换，都是矩阵运算的活。

# 矩阵乘法
A = np.random.rand(3, 3)
B = np.random.rand(3, 3)
C = np.dot(A, B)        # 或者 A @ B

# 求逆矩阵，计算倒格矢时必用
inv_A = np.linalg.inv(A)

# 特征值和特征向量，能带计算的核心
eigenvalues, eigenvectors = np.linalg.eig(A)

# 解线性方程组，做结构优化时常用
b = np.array([1, 2, 3])
x = np.linalg.solve(A, b)

我记得有一次做声子谱计算，需要反复求解大型线性方程组。用 np.linalg.solve 比手动求逆再乘要快一个数量级。为什么？因为 solve 用的是 LU 分解，避免了显式求逆的额外开销。

四、广播机制：隐形的性能加速器

广播机制是NumPy最强大的特性之一，也是初学者最容易踩坑的地方。说白了，就是让不同形状的数组能「自动对齐」进行运算。

# 经典例子：给每行加一个偏置
data = np.random.rand(5, 3)   # 5个样本，3个特征
bias = np.array([0.1, 0.2, 0.3])
result = data + bias          # 广播自动发生

# 材料计算中的实际应用：计算所有原子到某点的距离
positions = np.random.rand(100, 3)  # 100个原子的坐标
center = np.array([0.5, 0.5, 0.5])
distances = np.sqrt(np.sum((positions - center)**2, axis=1))
广播规则就三条：
  1. 如果维度不同，在左边补1
  2. 如果某个维度大小为1，就沿着该维度复制
  3. 如果维度大小既不相等也不为1，报错

我曾经在计算径向分布函数时，因为广播机制理解不到位，写了个三重循环，跑了10分钟。后来改成向量化+广播，0.5秒搞定。你想想看，这差距有多大。

五、随机数生成：蒙特卡洛模拟的基石

材料计算中，随机数用于蒙特卡洛模拟、随机初始结构生成、数据增强等场景。NumPy的随机模块功能很全。

# 设置种子，保证可重复性
np.random.seed(42)

# 均匀分布 [0,1)
rand_vals = np.random.rand(100)

# 正态分布，模拟热噪声
normal_vals = np.random.randn(100)

# 整数随机，用于随机采样
indices = np.random.randint(0, 100, size=10)

# 从给定分布中采样
# 比如生成符合玻尔兹曼分布的随机能量
beta = 1.0 / 0.025  # 室温下
energies = np.random.exponential(scale=1/beta, size=1000)
避坑指南： 我曾经在并行计算中用了默认的随机数生成器，结果每个进程都产生了相同的随机序列。后来改用 np.random.SeedSequence 为每个进程生成独立的种子，才解决了这个问题。

六、性能优化技巧：让代码飞起来

做材料计算，数据量动不动就是几百万个点。性能优化不是锦上添花，而是刚需。下面是我总结的几个实用技巧。

技巧 说明 加速比（经验值）
向量化运算 用数组运算代替循环 10~100倍
使用内置函数 用 np.sum 代替 Python sum 5~20倍
指定数据类型 用 float32 代替 float64 2倍（内存减半）
避免复制 用视图代替副本 视情况而定
使用 numba JIT编译加速循环 10~100倍 
# 坏写法：用循环
def slow_sum(data):
    total = 0
    for i in range(len(data)):
        total += data[i] ** 2
    return total

# 好写法：向量化
def fast_sum(data):
    return np.sum(data ** 2)

# 实测：100万个点，慢写法0.3秒，快写法0.003秒

我个人最常用的优化手段是「先向量化，再考虑数据类型」。如果向量化后还慢，再考虑用 np.float32 降低精度。做能带计算时，把电子密度从 float64 降到 float32，内存占用直接减半，计算速度也快了近一倍。

注意： 降低精度要谨慎。在需要高精度的场景（如求逆矩阵、解特征值），float32 可能导致数值不稳定。我建议先做小规模测试，确认误差在可接受范围内再大规模使用。

知识体系总览

下面这张图是我自己整理的NumPy知识体系，涵盖了本章的核心内容。你可以把它当作一个「导航图」，遇到问题时先看看属于哪个模块。

NumPy 知识体系总览 NumPy 核心 数组创建 np.array / np.zeros np.ones / np.linspace np.arange / np.eye 数组操作 切片 / 索引 reshape / transpose concatenate / split 矩阵运算 np.dot / @ 乘法 np.linalg.inv 求逆 eig / solve 广播机制 维度自动对齐 标量与数组运算 向量化计算 随机数 & 优化 np.random 模块 向量化 / 数据类型 numba JIT 加速 核心思想：用向量化思维替代循环思维 材料计算中，90% 的循环都可以用 NumPy 向量化操作替代

好了，以上就是NumPy在材料计算中的核心用法。说实话，这些内容我每次做新项目都会重新翻一遍。NumPy这东西，用多了自然就熟了。关键是要养成「向量化思维」——遇到循环先想想能不能用数组运算替代。这个习惯一旦养成，你的代码速度和可读性都会上一个台阶。

公众号:蓝海资料掘金营，微信deep3321

← 上一章 📖 返回目录 下一章 →