多元数据矩阵与向量化运算

小象实战讲义 · 多元统计分析

在多元统计分析中，我们面对的核心研究对象不再是单个变量，而是由多个变量构成的随机向量及其观测数据——多元数据矩阵。本节作为现代统计计算与数据科学的基础，将系统讲解如何从数学上描述多元数据，并介绍其核心的向量化运算逻辑。学完本节，你将能够理解多元数据的基本结构，掌握数据标准化、协方差矩阵与相关矩阵的计算原理，并熟练使用Python进行向量化操作，为后续的聚类、判别、降维等高级分析奠定坚实的计算基础。

💡 核心导读

• 多元数据矩阵：理解 $n \times p$ 数据矩阵 $X$ 的数学定义，其中 $n$ 为样本量，$p$ 为变量维数。
• 数据中心化与标准化：掌握将原始数据转化为零均值（中心化）和单位方差（Z-score标准化）的数学过程及其几何意义。
• 散度矩阵与协方差矩阵：推导样本协方差矩阵 $S$ 的两种等价计算公式，理解其作为变量间线性关系度量的核心作用。
• 相关矩阵：从标准化数据出发，定义并计算皮尔逊相关系数矩阵 $R$，理解其与协方差矩阵的内在联系。
• Python向量化实战：使用NumPy实现上述所有计算，体验向量化运算相比循环的巨大效率优势。

多元数据矩阵：结构与表示

在多元统计分析中，我们收集到的原始数据通常以表格形式呈现。设有 $p$ 个变量（特征），对 $n$ 个样本（观测对象）进行测量，则所有数据构成一个 $n \times p$ 的矩阵 $X$：

$$X = \begin{pmatrix} x_{11} & x_{12} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{np} \end{pmatrix} \in \mathbb{R}^{n \times p}$$

为了进行严格的向量运算，我们通常将每个样本视为一个 $p$ 维列向量。令第 $i$ 个样本的观测向量为 $\mathbf{x}^{(i)} \in \mathbb{R}^p$，它是数据矩阵 $X$ 第 $i$ 行的转置：

$$\mathbf{x}^{(i)} = (x_{i1}, x_{i2}, \ldots, x_{ip})' = \begin{pmatrix} x_{i1} \\ x_{i2} \\ \vdots \\ x_{ip} \end{pmatrix}$$

这样，$X'$ 可看作由 $n$ 个样本列向量 $\mathbf{x}^{(i)}$ 水平堆叠而成（即 $X' = [\mathbf{x}^{(1)},\ldots,\mathbf{x}^{(n)}]$，为 $p \times n$ 矩阵），而 $X$ 本身则是对应行向量 $(\mathbf{x}^{(i)})'$ 的纵向堆叠（尽管存储时是行优先）。第 $j$ 个变量的 $n$ 次观测则构成列向量 $\mathbf{x}_j = (x_{1j}, x_{2j}, \ldots, x_{nj})' \in \mathbb{R}^n$。

样本均值向量 $\bar{\mathbf{x}} \in \mathbb{R}^p$ 概括了所有变量的中心位置，通过对所有样本向量求平均得到：

$$\bar{\mathbf{x}} = \frac{1}{n} \sum_{i=1}^{n} \mathbf{x}^{(i)} = \begin{pmatrix} \bar{x}_1 \\ \bar{x}_2 \\ \vdots \\ \bar{x}_p \end{pmatrix}, \quad \text{其中 } \bar{x}_j = \frac{1}{n} \sum_{i=1}^{n} x_{ij}$$

数据中心化、标准化与散度矩阵

数据中心化

数据中心化是许多多元分析的第一步，其目的是使数据围绕原点分布、消除各变量均值水平（绝对数值）的差异；注意中心化只平移数据位置，不改变方差与量纲，量纲的消除由随后的Z-score标准化完成。定义中心化矩阵 $H$：

$$H = I_n - \frac{1}{n} \mathbf{1}_n \mathbf{1}_n'$$

其中 $I_n$ 是 $n$ 阶单位矩阵，$\mathbf{1}_n = (1,1,\ldots,1)' \in \mathbb{R}^n$。$H$ 是一个幂等对称矩阵（$H^2 = H$），其作用是将任意向量投影到与 $\mathbf{1}_n$ 正交的子空间。

对原始数据矩阵 $X$ 进行中心化，得到中心化数据矩阵 $X_c$：

$$X_c = H X = X - \mathbf{1}_n \bar{\mathbf{x}}'$$

$X_c$ 的每一列（即每个变量的观测向量）均值均为0。

数据标准化（Z-score）

为进一步消除各变量量纲和方差差异的影响，常进行Z-score标准化。令 $s_j$ 为第 $j$ 个变量的样本标准差：

$$s_j = \sqrt{\frac{1}{n-1} \sum_{i=1}^{n} (x_{ij} - \bar{x}_j)^2}$$

则标准化后的元素为：

$$z_{ij} = \frac{x_{ij} - \bar{x}_j}{s_j}$$

标准化后数据矩阵 $Z$ 的每一列均值为0，标准差为1。

散度矩阵与协方差矩阵

样本散度矩阵 $T$（也称为总平方和与叉积矩阵）定义为所有中心化样本向量的外积之和：

$$T = \sum_{i=1}^{n} (\mathbf{x}^{(i)} - \bar{\mathbf{x}})(\mathbf{x}^{(i)} - \bar{\mathbf{x}})' \in \mathbb{R}^{p \times p}$$

利用矩阵形式，可等价写为 $T = X_c' X_c$。

样本协方差矩阵 $S$ 是散度矩阵的无偏估计：

$$S = \frac{1}{n-1} T = \frac{1}{n-1} \sum_{i=1}^{n} (\mathbf{x}^{(i)} - \bar{\mathbf{x}})(\mathbf{x}^{(i)} - \bar{\mathbf{x}})' \in \mathbb{R}^{p \times p}$$

其元素 $s_{jk}$ 为变量 $j$ 与 $k$ 的样本协方差：

$$s_{jk} = \frac{1}{n-1} \sum_{i=1}^{n} (x_{ij} - \bar{x}_j)(x_{ik} - \bar{x}_k)$$

相关矩阵

从标准化数据 $Z$ 出发，可计算样本相关矩阵 $R$：

$$R = \frac{1}{n-1} Z' Z \in \mathbb{R}^{p \times p}$$

$R$ 的元素 $r_{jk}$ 即变量 $j$ 与 $k$ 的皮尔逊相关系数：

$$r_{jk} = \frac{s_{jk}}{s_j s_k} = \frac{\sum_{i=1}^{n} (x_{ij} - \bar{x}_j)(x_{ik} - \bar{x}_k)}{\sqrt{\sum_{i=1}^{n} (x_{ij} - \bar{x}_j)^2 \sum_{i=1}^{n} (x_{ik} - \bar{x}_k)^2}}$$

Python实战：多元数据矩阵运算

下面我们使用Python的NumPy库，完整演示从生成多元数据到计算协方差矩阵和相关矩阵的全过程。代码严格遵循上述数学公式，并验证与NumPy内置函数结果的一致性。

"""
实战演示 2.1：多元数据矩阵运算
演示中心化、标准化、协方差矩阵与相关矩阵的手动计算与验证
"""
import numpy as np

# 固定随机种子，确保结果可复现
np.random.seed(42)

# 生成模拟多元数据：n=50个样本，p=4个变量
# 假设数据来自均值为10，标准差为3的正态分布
n, p = 50, 4
X = np.random.randn(n, p) * 3 + 10  # 形状 (50, 4)

print("数据矩阵形状:", X.shape)
print("前5行数据预览:\n", np.round(X[:5], 2))

# 1. 计算样本均值向量 (每个变量的均值)
mean_vector = np.mean(X, axis=0)  # 形状 (4,)
print("\n样本均值向量:", np.round(mean_vector, 4))

# 2. 数据中心化：使用中心化矩阵 H
# H = I_n - (1/n) * 1_n * 1_n'
I_n = np.eye(n)
ones_n = np.ones((n, 1))
H = I_n - np.ones((n, n)) / n  # 等价于 H = I_n - (1/n)*ones_n@ones_n.T

X_centered = H @ X  # 中心化数据矩阵
print("\n中心化后数据均值验证:", np.round(np.mean(X_centered, axis=0), 8))

# 3. 计算样本协方差矩阵 S (两种方法)
# 方法1：基于中心化数据手动计算 S = (X_c' X_c) / (n-1)
S_manual = (X_centered.T @ X_centered) / (n - 1)

# 方法2：使用NumPy内置函数
S_numpy = np.cov(X, rowvar=False)  # rowvar=False 表示每列是一个变量

# 验证两种方法结果一致性
diff_norm = np.linalg.norm(S_manual - S_numpy)
print("\n手动计算与 numpy.cov 结果差异范数:", round(float(diff_norm), 8))

# 4. 数据标准化 (Z-score)
std_devs = np.std(X, axis=0, ddof=1)  # 样本标准差，ddof=1 对应 n-1 分母
Z = (X - np.mean(X, axis=0)) / std_devs

print("\n标准化数据 Z 的均值:", np.round(np.mean(Z, axis=0), 4).tolist())
print("标准化数据 Z 的标准差:", np.round(np.std(Z, axis=0, ddof=1), 4).tolist())

# 5. 计算相关矩阵 R (两种方法)
# 方法1：基于标准化数据 R = (Z' Z) / (n-1)
R_from_Z = (Z.T @ Z) / (n - 1)

# 方法2：使用NumPy内置函数
R_numpy = np.corrcoef(X, rowvar=False)

# 验证两种方法结果一致性
max_diff = np.max(np.abs(R_from_Z - R_numpy))
print("\n两种方法计算的相关矩阵最大绝对差异:", round(float(max_diff), 8))

# 6. 展示协方差矩阵和相关矩阵（保留3位小数）
print("\n样本协方差矩阵 S (前3x3子矩阵):")
print(np.round(S_numpy[:3, :3], 3))

print("\n样本相关矩阵 R (前3x3子矩阵):")
print(np.round(R_numpy[:3, :3], 3))

运行上述代码，你将得到如下输出（与提供的验证输出一致）：

数据矩阵形状: (50, 4)
前5行数据预览:
 [[11.49  9.59 11.94 14.57]
 [ 9.3   9.3  14.74 12.3 ]
 [ 8.59 11.63  8.61  8.6 ]
 [10.73  4.26  4.83  8.31]
 [ 6.96 10.94  7.28  5.76]]

样本均值向量: [ 9.7304  9.9968  9.5762 10.2073]

中心化后数据均值验证: [ 0. -0. -0.  0.]

手动计算与 numpy.cov 结果差异范数: 0.0

标准化数据 Z 的均值: [-0.0, 0.0, -0.0, -0.0]
标准化数据 Z 的标准差: [1.0, 1.0, 1.0, 1.0]

两种方法计算的相关矩阵最大绝对差异: 0.0

样本协方差矩阵 S (前3x3子矩阵):
[[ 4.361  0.43  -0.506]
 [ 0.43   8.152 -0.722]
 [-0.506 -0.722  8.96 ]]

样本相关矩阵 R (前3x3子矩阵):
[[ 1.     0.072 -0.081]
 [ 0.072  1.    -0.084]
 [-0.081 -0.084  1.   ]]

📝 动手练一练

1. 协方差矩阵的性质验证：基于上面生成的协方差矩阵 $S$，请验证： a) $S$ 是否对称？ b) $S$ 的对角线元素 $s_{jj}$ 与各变量的样本方差 $s_j^2$ 是否相等？ c) 计算 $S$ 的行列式 $\det(S)$ 和迹 $\operatorname{tr}(S)$。

2. 从相关矩阵反推协方差矩阵：假设你只有相关矩阵 $R$ 和各变量的标准差向量 $\mathbf{s} = (s_1, s_2, s_3, s_4)$，请写出从 $R$ 和 $\mathbf{s}$ 重构协方差矩阵 $S$ 的公式，并用Python代码验证。

参考答案：

1. a) 是，协方差矩阵总是对称的，因为 $s_{jk} = s_{kj}$。 b) 是，对角线元素 $s_{jj} = \frac{1}{n-1}\sum_i (x_{ij}-\bar{x}_j)^2 = s_j^2$。 c) 使用代码 np.linalg.det(S_numpy) 和 np.trace(S_numpy) 计算。
2. 重构公式：$S = \text{diag}(\mathbf{s}) \cdot R \cdot \text{diag}(\mathbf{s})$，其中 $\text{diag}(\mathbf{s})$ 是以 $\mathbf{s}$ 为对角线元素的对角矩阵。验证代码：

   s = np.std(X, axis=0, ddof=1)
   D = np.diag(s)
   S_reconstructed = D @ R_numpy @ D
   print("重构误差:", np.linalg.norm(S_reconstructed - S_numpy))

本章小结

本节系统建立了多元数据矩阵的数学表示与基本运算框架：

1. 多元数据矩阵 $X \in \mathbb{R}^{n \times p}$ 是多元分析的起点，其行对应样本，列对应变量。每个样本应明确表示为 $p$ 维列向量 $\mathbf{x}^{(i)}$，以确保后续向量运算的维度自洽。
2. 中心化与标准化是数据预处理的基石，分别通过减去均值、除以标准差实现，为后续分析提供统一尺度。
3. 协方差矩阵 $S$ 和相关矩阵 $R$ 是描述变量间线性关系的核心工具，$S$ 包含方差与协方差信息，$R$ 是标准化后的 $S$，专门度量相关性。
4. 向量化运算是高效计算的关键，利用矩阵乘法一次性完成所有样本的聚合计算，比循环快数个数量级。

行动清单：

• 在Python中生成你自己的多元数据集，重复本节所有计算步骤，特别注意样本向量的列向量表示。
• 尝试修改数据生成参数（如均值、方差、样本量），观察协方差矩阵和相关矩阵的变化。
• 查阅NumPy官方文档，深入了解 np.cov、np.corrcoef 等函数的参数与用法。

— 小象教研组