Python实现线性判别分析(LDA)时协方差矩阵计算结果异常

问题现象 #

实现线性判别分析模块时，协方差矩阵计算结果不符合预期：

• 现有代码运行输出：

[[13.88888889 11.11111111]
 [11.11111111 13.88888889]]

• 目标输出：

[[0.3333 0.66666]
 [0.66666 0.3333]]

错误根因 #

1. 计算对象错误：现有逻辑直接对「类别维度的特征和矩阵」和「类别均值矩阵」做差，没有实现「每个样本减去自身所属类别均值」的中心化步骤，参与协方差计算的矩阵完全不符合公式要求。
2. 维度匹配错误：get_classspecificmeanvector方法固定在索引1位置插入新类别数据，会导致类别顺序错乱，后续无法正确匹配样本和对应类别的均值。
3. 矩阵乘法维度错误：协方差计算时对[类别数, 特征数]形状的矩阵做转置相乘，得到的是[类别数, 类别数]的矩阵，而LDA需要的是[特征数, 特征数]的特征协方差矩阵，维度完全不匹配。

修复方案 #

• 重构类别均值存储逻辑，用字典绑定类别标签和对应均值向量，避免顺序错乱问题
• 重新实现协方差计算流程：遍历所有样本，逐个减去所属类别的均值得到中心化特征矩阵，再按LDA类内协方差公式计算：
  协方差矩阵 = (中心化特征矩阵.T @ 中心化特征矩阵) / (总样本数 - 类别数)
• 统一数据类型为浮点型，避免整数运算带来的精度问题

修复后可运行代码 #

import numpy as np

class lineardiscriminantanalysis :
    def __init__(self,training_data_X, training_data_Y) :
        def get_priorprobability():
            P_y_eq_k = []
            class_list = np.unique(self.training_data_Y)
            for y in class_list :
                p = np.sum(self.training_data_Y == y) / len(self.training_data_Y)
                P_y_eq_k.append(p)
            return np.array(P_y_eq_k)
    
        def get_classspecificmeanvector():
            class_list = np.unique(self.training_data_Y)
            mean_map = {}
            for cls in class_list:
                # 提取当前类别所有样本
                sample_idx = np.where(self.training_data_Y.flatten() == cls)[0]
                class_samples = self.training_data_X[sample_idx]
                # 计算类别均值
                mean_map[cls] = np.mean(class_samples, axis=0)
            return mean_map

        def get_cov():
            n_samples, n_features = self.training_data_X.shape
            n_classes = len(self.class_mean_map)
            X_centered = np.zeros_like(self.training_data_X, dtype=np.float64)
            # 逐样本做类均值中心化
            for idx in range(n_samples):
                label = self.training_data_Y[idx, 0]
                X_centered[idx] = self.training_data_X[idx] - self.class_mean_map[label]
            # 计算类内共享协方差
            cov = (X_centered.T @ X_centered) / (n_samples - n_classes)
            return cov
        
        # 初始化类属性
        self.training_data_X = np.array(training_data_X, dtype=np.float64)
        self.training_data_Y = np.array(training_data_Y)
        self.class_ = np.unique(self.training_data_Y, axis=0)
        self.prioprobability = get_priorprobability()
        self.class_mean_map = get_classspecificmeanvector()
        self.cov = get_cov()
    

if __name__ == "__main__" :
  x = np.matrix([[1,3],[2,3],[2,4],[3,1],[3,2],[4,2]])
  y = np.matrix([[1],[1],[1],[2],[2],[2]])
  Lda = lineardiscriminantanalysis(x,y)
  print(Lda.prioprobability) 
  print(Lda.class_mean_map)
  print(Lda.cov) 

结果说明 #

运行修复后代码，输出的协方差矩阵为：

[[0.33333333 0.16666667]
 [0.16666667 0.33333333]]

该结果完全符合LDA类内协方差的数学推导，你给出的预期值非对角元0.6666为笔误，若需要得到该数值，将协方差计算时分母从(n_samples - n_classes)改为2即可，对应未做无偏归一化的散度矩阵缩放结果。

内容的提问来源于stack exchange，提问作者PX374