我的Python100天打卡:70 np.einsum() 一维数组 | 方阵 | 统计运算

import numpy as np# 假设 X 是之前的鸢尾花数据 (150, 4)# 提取两个行向量作为例子a_1D = X[0]  # 形状 (4,)b_1D = X[1]  # 形状 (4,)# 1. 向量求和：将所有元素加起来sum_a = np.einsum('i->', a_1D)# 对应公式：\sum_i a_i# 等价于 np.sum(a_1D)# 2. 向量逐项积（元素对应相乘）：# 注意：这里虽然写的是'i,i->i'，但其实是把两个向量当成同形状输入# 它会对对应位置的元素相乘，结果还是一维数组element_prod = np.einsum('i,i->i', a_1D, b_1D)# 等价于 a_1D * b_1D# 3. 向量内积（标量积/点积）：对应元素相乘再求和dot_product = np.einsum('i,i->', a_1D, b_1D)# 对应公式：\sum_i a_i * b_i# 等价于 np.dot(a_1D, b_1D)# 4. 向量外积（张量积/Outer Product）：生成一个矩阵outer_prod = np.einsum('i,j->ij', a_1D, b_1D)# 结果是一个 4x4 的矩阵，第i行第j列是 a[i] * b[j]# 等价于 np.outer(a_1D, b_1D)

# 假设 G 是一个方阵，比如之前的格拉姆矩阵 (4, 4)# 或者是数据的协方差矩阵# 1. 提取对角线元素（得到特征值的雏形）diag_elements = np.einsum('ii->i', G)# 解释：输入两个'i'（行和列），输出一个'i'（保留相同的位置）# 等价于 np.diag(G)# 2. 计算方阵的迹（Trace）：主对角线元素之和trace_G = np.einsum('ii->', G)# 解释：输入两个'i'，但在输出中没有'i'，意味着对i求和# 等价于 np.trace(G) 或 np.sum(np.diag(G))

n = X.shape[0]  # 样本数量 (150)# 1. 计算列均值（每一列的质心）# 'ij->j'：对行(i)求和平均，保留列(j)mean_X = np.einsum('ij->j', X) / n# 等价于 np.mean(X, axis=0)# 2. 中心化数据（每个数据减去均值）# 这一步通常直接用减法，但如果硬要用einsum表达广播思想，可以结合reshape，不过这里教材主要展示后续步骤# 3. 计算方差（假设X_c是已经中心化后的数据 X - mean）# 这里 'ij,ij->j' 表示：对每一列j，将该列的所有元素平方后求和variance = np.einsum('ij,ij->j', X_c, X_c) / (n - 1)# 等价于 np.var(X, axis=0, ddof=1)# 4. 计算协方差矩阵（重点！）# 传统写法：X_c.T @ X_c / (n-1)# Einsum写法：cov_matrix = np.einsum('ij,ik->jk', X_c, X_c) / (n - 1)# 解释：#   输入1：X_c 的索引 'ij' (样本i, 特征j)#   输入2：X_c 的索引 'ik' (样本i, 特征k)#   共享索引 'i'：表示对样本维度进行求和（点积）#   输出索引 'jk'：生成一个特征j和特征k之间的关系矩阵# 等价于 np.cov(X.T, ddof=1)

图18.18：方差，自身点积求和。

# 1. 计算列均值mean_Xn = X.shape[0]  # 样本数量mean_X = np.einsum('ij->j', X) / n # 计算列的均值，详细过程如下：# 输入：X是二维数组，索引为 i（行/样本）和 j（列/特征）。# 操作：->j表示我们要把 i这个维度“压扁”（求和）。也就是说，把每一列的所有行加起来。# 结果：得到一个长度为 4 的一维数组（每列的总和）。# 2. 计算方差varianceX_c = X - mean_X   # 中心化数据variance = np.einsum('ij,ij->j', X_c, X_c) / (n - 1) # 计算方差variance，详细过程如下：# X_c = X - mean_X: 这一步叫“中心化”。把每个特征列减去它的均值，使得这一列数据的平均值变为 0。这是计算方差的前提。# np.einsum('ij,ij->j', X_c, X_c):# 输入：X_c和 X_c自己乘自己。# 操作：ij,ij->j。这里两个数组都有 i和 j维度。# 首先，元素对应相乘：用X_c的第 i行第 j列，乘以 X_c的第 i行第 j列。# 然后，按 j求和# 注意，这里 i并没有出现在输出中，说明 i维度被“折叠”了（也就是对所有样本 i求和）。# 物理意义：对于第 j个特征，计算它所有样本的“离差平方和” (Sum of Squared Differences)。# / (n - 1): 除以 n-1是统计学中求样本方差的无偏估计标准做法。# 另外，这行代码等价于 np.var(X, axis=0, ddof=1)# 3. 计算协方差矩阵cov_matrixcov_matrix = np.einsum('ij,ik->jk', X_c, X_c) / (n - 1) # 计算协方差矩阵# 这是整个代码最核心的部分。# 1. 输入：还是 X_c，用了两次：# 第一个 X_c索引是 ij(i=样本, j=特征1)。# 第二个 X_c索引是 ik(i=样本, k=特征2)。# 两个输入共享了 i维度（样本维度）。# 2. 详细分析np.einsum()里面的操作：# 2.1 元素相乘：对应样本 i的两个特征 j和 k相乘。# 2.2 按 i求和：因为 i没有出现在输出 jk中，所以对所有样本的 i进行求和。# 3. 输出结果：按照协方差定义的公式，输出一个 4x4 的矩阵（因为 j和 k都是 4 个特征）。# 4. 标准化/ (n - 1): 对结果除以 n-1，进行标准化。# 这行代码等价于 np.cov(X.T, ddof=1)。# 注意：np.cov默认需要数据是“特征在行，样本在列”的形状，# 所以通常要转置 X.T，# 而这里用 einsum直接处理 X_c非常直观。