典型相关系数阵计算与代码实战

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

在多元统计分析中，当两组指标存在不同量纲或数量级差异显著时，若直接使用原始协方差矩阵求解典型相关，方差较大的变量将会主导分析结果，从而扭曲真实的统计关联。因此，基于数据标准化后的相关系数矩阵进行典型相关分析（Canonical Correlation Analysis, CCA）是工业界与学术界的通用标准。本节系统推导由相关系数矩阵求解典型相关系数的特征值方程，并基于 Python 与 sklearn 库完整实现典型相关系数、典型变量得分及因子载荷的计算与实战解读。

💡 核心导读

• 消除量纲影响：通过 $Z$-score 标准化消除指标量纲与数量级差异，由样本相关阵 $R$ 代替协方差阵 $\Sigma$
• 分块矩阵形式：将相关系数矩阵分块为 $R_{11}, R_{12}, R_{21}, R_{22}$，明确组内与组间关联结构
• 代数特征值方程：典型相关系数平方 $\lambda = \rho^2$ 为矩阵 $M = R_{11}^{-1}R_{12}R_{22}^{-1}R_{21}$ 的非零特征值
• 典型载荷与结构矩阵：通过原始变量与典型变量的相关系数解释提取因子的实际业务含义
• Python 实战落地：基于 sklearn.cross_decomposition.CCA 完成从模型拟合到典型载荷提取的全流程

一、从协方差矩阵到相关系数矩阵

1.1 为何必须进行标准化？

设两组随机向量分别为 $X = (X_1, X_2, \dots, X_p)^\top$ 和 $Y = (Y_1, Y_2, \dots, Y_q)^\top$。在实际业务中，$X$ 组可能是人体的生理指标（如身高以厘米计、体重以公斤计、血压以毫米汞柱计），而 $Y$ 组可能是运动机能指标（如肺活量以毫升计、引体向上以次数计）。

不同变量的方差可能相差数个数量级。若不对数据做无量纲化处理，典型相关分析所求得的线性组合系数将严重受各变量方差大小的牵制。

为此，对各变量进行零均值、单位方差的标准化处理： $$X_i^* = \frac{X_i - \mu_{Xi}}{\sigma_{Xi}}, \quad Y_j^* = \frac{Y_j - \mu_{Yj}}{\sigma_{Yj}}$$

经过标准化后，变量的协方差矩阵即转化为样本相关系数矩阵： $$R = \begin{pmatrix} R_{11} & R_{12} \\ R_{21} & R_{22} \end{pmatrix}$$ 其中：

• $R_{11} \in \mathbb{R}^{p \times p}$ 为第一组变量自身的互相关阵，对角线全为 1；
• $R_{22} \in \mathbb{R}^{q \times q}$ 为第二组变量自身的互相关阵，对角线全为 1；
• $R_{12} = R_{21}^\top \in \mathbb{R}^{p \times q}$ 为两组变量之间的互相关阵。

二、典型相关变量的代数求解推导

2.1 优化问题构建

设标准化变量的线性组合为典型变量： $$U = a^\top X^* = \sum_{i=1}^p a_i X_i^*, \quad V = b^\top Y^* = \sum_{j=1}^q b_j Y_j^*$$

要求解权向量 $a \in \mathbb{R}^p$ 和 $b \in \mathbb{R}^q$，使得 $U$ 与 $V$ 的相关系数最大化： $$\max_{a, b} \text{Corr}(U, V) = \frac{\text{Cov}(U, V)}{\sqrt{\text{Var}(U)\text{Var}(V)}} = \frac{a^\top R_{12} b}{\sqrt{(a^\top R_{11} a)(b^\top R_{22} b)}}$$

由于线性组合乘以非零常数不改变相关系数，我们引入方差归一化约束条件： $$a^\top R_{11} a = 1, \quad b^\top R_{22} b = 1$$ 此时目标函数转化为在约束条件下最大化协方差： $$\max_{a, b} \quad a^\top R_{12} b$$

2.2 特征方程的建立

构造拉格朗日函数： $$L(a, b, \lambda, \mu) = a^\top R_{12} b - \frac{\lambda}{2}(a^\top R_{11} a - 1) - \frac{\mu}{2}(b^\top R_{22} b - 1)$$

分别对向量 $a$ 和 $b$ 求一阶偏导并令其为零向量： $$\frac{\partial L}{\partial a} = R_{12} b - \lambda R_{11} a = \mathbf{0} \implies R_{12} b = \lambda R_{11} a$$ $$\frac{\partial L}{\partial b} = R_{21} a - \mu R_{22} b = \mathbf{0} \implies R_{21} a = \mu R_{22} b$$

将第一式左乘 $a^\top$，第二式左乘 $b^\top$，结合约束条件 $a^\top R_{11} a = 1$ 与 $b^\top R_{22} b = 1$： $$a^\top R_{12} b = \lambda = \mu = \rho$$ 即拉格朗日乘子 $\lambda = \mu = \rho$，它恰好等于典型相关系数！

假定 $R_{11}$ 与 $R_{22}$ 均为可逆满秩阵，由第二式解得： $$b = \frac{1}{\rho} R_{22}^{-1} R_{21} a$$ 代入第一式中，消去 $b$ 可得： $$R_{12} \left(\frac{1}{\rho} R_{22}^{-1} R_{21} a\right) = \rho R_{11} a$$ 两边同时左乘 $R_{11}^{-1}$，整理得标准特征值方程： $$(R_{11}^{-1} R_{12} R_{22}^{-1} R_{21}) a = \rho^2 a$$

同理，若先消去 $a$，可得关于向量 $b$ 的特征值方程： $$(R_{22}^{-1} R_{21} R_{11}^{-1} R_{12}) b = \rho^2 b$$

2.3 结论与性质

1. 矩阵 $M_1 = R_{11}^{-1} R_{12} R_{22}^{-1} R_{21}$ 与 $M_2 = R_{22}^{-1} R_{21} R_{11}^{-1} R_{12}$ 具有完全相同的非零特征值；
2. 非零特征值的个数为 $k = \min(p, q)$，且所有特征值落在区间 $[0, 1)$ 内；
3. 特征值按降序排列：$\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_k \ge 0$；
4. 第 $i$ 对典型相关系数即为 $\rho_i = \sqrt{\lambda_i}$，对应的特征向量经过归一化后即为典型系数向量 $a_i$ 与 $b_i$。

三、典型变量得分与典型结构矩阵

在完成典型系数的计算后，需要对典型变量代表的统计学含义进行解释。通常借助两类重要矩阵：

3.1 典型变量得分（Canonical Variate Scores）

将标准化样本矩阵 $X^*$ 与 $Y^*$ 分别代入系数向量： $$u_i = X^* a_i, \quad v_i = Y^* b_i$$ 向量 $u_i$ 与 $v_i$ 记录了每个样本在第 $i$ 对典型变量上的得分，可用于后续的低维投影作图与聚类分析。

3.2 典型载荷矩阵（Canonical Loadings）

典型载荷定义为原始标准化变量与典型变量之间的简单相关系数： $$\text{Corr}(X^*, U_i) = \mathbb{E}[X^* U_i] = \mathbb{E}[X^* (a_i^\top X^*)] = \mathbb{E}[X^* (X^*)^\top] a_i = R_{11} a_i$$ $$\text{Corr}(Y^*, V_i) = \mathbb{E}[Y^* V_i] = \mathbb{E}[Y^* (b_i^\top Y^*)] = \mathbb{E}[Y^* (Y^*)^\top] b_i = R_{22} b_i$$ 载荷绝对值越接近 1，表示该原始变量受该典型变量的影响越大，这对解释典型变量的实际物理含义至关重要。

四、Python 代码实战：基于 sklearn 的 CCA

在工程实践中，推荐直接使用 sklearn.cross_decomposition.CCA。下面我们通过模拟两组含有共性潜变量的高维指标，演示从拟合、相关系数计算到典型载荷分析的全流程。

import numpy as np
from sklearn.cross_decomposition import CCA

np.random.seed(42)
n = 120
t1 = np.random.normal(0, 1, size=n)
t2 = np.random.normal(0, 1, size=n)

X = np.column_stack([
    0.8 * t1 + np.random.normal(0, 0.4, size=n),
    0.7 * t1 + np.random.normal(0, 0.5, size=n),
    0.4 * t2 + np.random.normal(0, 0.6, size=n)
])

Y = np.column_stack([
    0.75 * t1 + np.random.normal(0, 0.4, size=n),
    0.6 * t1 + 0.5 * t2 + np.random.normal(0, 0.4, size=n)
])

cca = CCA(n_components=2)
cca.fit(X, Y)
X_c, Y_c = cca.transform(X, Y)

r1 = np.corrcoef(X_c[:, 0], Y_c[:, 0])[0, 1]
r2 = np.corrcoef(X_c[:, 1], Y_c[:, 1])[0, 1]

x_loadings = np.array([np.corrcoef(X[:, j], X_c[:, 0])[0, 1] for j in range(3)])
y_loadings = np.array([np.corrcoef(Y[:, j], Y_c[:, 0])[0, 1] for j in range(2)])

print("第一典型相关对相关系数:", round(float(r1), 4))
print("第二典型相关对相关系数:", round(float(r2), 4))
print("X 组各变量在第一典型变量上的载荷:", np.round(x_loadings, 4).tolist())
print("Y 组各变量在第一典型变量上的载荷:", np.round(y_loadings, 4).tolist())

代码输出与结果解读

运行上述脚本，控制台输出如下统计结果：

第一典型相关对相关系数: 0.8074
第二典型相关对相关系数: 0.4609
X 组各变量在第一典型变量上的载荷: [0.9542, 0.7928, 0.2911]
Y 组各变量在第一典型变量上的载荷: [0.8887, 0.9006]

实战结果解读：

1. 典型相关强度：第一典型相关对的相关系数高达 0.8074，表明两组变量之间存在极强的线性共线性；第二对典型相关系数降至 0.4609，反映了次级变异结构的关联。
2. X 组变量载荷结构：变量 $X_1$ 和 $X_2$ 在第一典型变量上的载荷分别达到 0.9542 和 0.7928，而 $X_3$ 仅为 0.2911。这清晰表明第一典型变量主要由前两个指标驱动。
3. Y 组变量载荷结构：变量 $Y_1$ 与 $Y_2$ 在第一典型变量上的载荷均在 0.88 以上，说明第一典型变量对准确捕获了两组变量受公共潜在信号 $t_1$ 驱动的整体联动特征。

📝 动手练一练

1. 特征值方程推导：已知 $X$ 组有 2 个变量，$Y$ 组有 3 个变量，且两组变量相关系数矩阵非奇异。请问矩阵 $M = R_{11}^{-1} R_{12} R_{22}^{-1} R_{21}$ 的阶数是多少？其非零特征值最多有几个？
2. 典型载荷计算：在某典型相关分析中，第一典型变量权向量为 $a_1 = (0.6, 0.8)^\top$，$X$ 组变量的相关系数矩阵为 $R_{11} = \begin{pmatrix} 1.0 & 0.5 \\ 0.5 & 1.0 \end{pmatrix}$。请计算两个原始变量在第一典型变量上的因子载荷。

本章小结

本节重点讨论了基于相关系数矩阵求解典型相关分析的理论与工程实战：

1. 标准化必要性：多元统计中消除量纲和方差悬殊的关键在于基于样本相关阵开展分析；
2. 代数本质：典型相关系数的平方等价于广义特征值方程的特征值，典型权向量即为对应的特征向量；
3. 典型载荷：通过原始变量与典型变量的互相关（因子载荷），为抽象的典型综合指标赋予具体的业务涵义；
4. Python 工具链：熟练运用 sklearn.cross_decomposition.CCA 即可高效完成模型训练与指标提取。

📋 行动清单

• 用 Python 对两组变量分别实施 Z-score 标准化并构建相关系数矩阵 R 替代原始协方差阵。
• 将 R 分块为 R11、R12、R21、R22 后求解特征值方程，输出典型相关系数与典型系数向量。
• 调用 sklearn 的 CCA 或自编代码计算典型变量得分和因子载荷，检查标准化后变量贡献是否受量纲影响。

—— 小象教研组