LDA中缩放矩阵未归一化的原因及特征分解实现方法

在执行线性判别分析（LDA）时，发现R语言MASS包的lda()函数生成的缩放矩阵未归一化，示例如下：

(res <- MASS::lda(Species~., iris))

输出结果：

Call:
lda(Species ~ ., data = iris)

Prior probabilities of groups:
    setosa versicolor  virginica 
 0.3333333  0.3333333  0.3333333 

Group means:
           Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa            5.006       3.428        1.462       0.246
versicolor        5.936       2.770        4.260       1.326
virginica         6.588       2.974        5.552       2.026

Coefficients of linear discriminants:
                    LD1         LD2
Sepal.Length  0.8293776  0.02410215
Sepal.Width   1.5344731  2.16452123
Petal.Length -2.2012117 -0.93192121
Petal.Width  -2.8104603  2.83918785

Proportion of trace:
   LD1    LD2 
0.9912 0.0088 

对上述缩放矩阵进行归一化后的结果：

scale(res$scaling, F, sqrt(colSums(res$scaling^2)))

输出结果：

LD1          LD2
Sepal.Length  0.2087418  0.006531964
Sepal.Width   0.3862037  0.586610553
Petal.Length -0.5540117 -0.252561540
Petal.Width  -0.7073504  0.769453092
attr("scaled:scale")
     LD1      LD2 
3.973222 3.689878 

手动通过特征分解拟合LDA得到的缩放矩阵是归一化的：

x <- scale(as.matrix(iris[,-5]), TRUE, FALSE)
y <- iris[,5]
means <- tapply(x,list(rep(y,ncol(x)), col(x)), mean)
Swithin <- crossprod(x - means[y,])
Sbetween <- crossprod(means)
eig <- eigen(solve(Swithin, Sbetween))
eig[[2]][,eig[[1]] > 1e-8]

输出结果：

[,1]         [,2]
[1,]  0.2087418 -0.006531964
[2,]  0.3862037 -0.586610553
[3,] -0.5540117  0.252561540
[4,] -0.7073504 -0.769453092

两者仅存在缩放因子差异，但该差异会影响后验概率。观察R源码发现lda()函数采用两次SVD而非特征分解实现，模拟该逻辑得到的结果与未归一化缩放矩阵一致：

a <- svd((x - means[y, ])/sqrt(nrow(x) - nrow(means)))
S1 <-  a$v %*% diag(1/a$d)
S1 %*% svd(means %*% S1)$v

输出结果：

[,1]        [,2]      [,3]
[1,]  0.8293776  0.02410215 -3.176869
[2,]  1.5344731  2.16452123  1.965956
[3,] -2.2012117 -0.93192121  2.076870
[4,] -2.8104603  2.83918785 -1.447218

核心问题解答 #

1. LDA中缩放矩阵未归一化的原因是什么？ #

MASS包的lda()函数采用两次SVD的实现方式，而非传统的特征分解：

• 第一步对组内离差矩阵进行SVD分解，得到其逆的平方根矩阵S1；
• 第二步对中心化数据投影到S1后的组间均值矩阵做SVD，最终缩放矩阵是S1与该SVD右奇异值矩阵的乘积。
  这种实现没有对最终判别系数做单位长度的归一化约束，因此得到的缩放矩阵是非归一化的。而传统特征分解方法中，特征向量默认是单位长度的，所以结果是归一化的。

2. 如何通过特征分解得到该未归一化的缩放矩阵？ #

可以在特征分解得到的归一化向量基础上，乘以对应判别方向的范数（即lda()输出缩放矩阵的列范数），步骤如下：

1. 通过传统特征分解得到归一化特征向量矩阵；
2. 计算lda()输出缩放矩阵中每列的范数；
3. 将归一化特征向量的每列对应乘以范数，即可得到与lda()一致的未归一化缩放矩阵（允许符号差异，因为特征向量方向不影响判别结果）。

示例代码：

# 传统特征分解得到归一化向量
x <- scale(as.matrix(iris[,-5]), TRUE, FALSE)
y <- iris[,5]
means <- tapply(x,list(rep(y,ncol(x)), col(x)), mean)
Swithin <- crossprod(x - means[y,])
Sbetween <- crossprod(means)
eig <- eigen(solve(Swithin, Sbetween))
norm_vecs <- eig[[2]][,eig[[1]] > 1e-8]

# 获取lda缩放矩阵的列范数
res <- MASS::lda(Species~., iris)
scales <- sqrt(colSums(res$scaling^2))

# 生成未归一化的缩放矩阵，调整符号匹配lda输出
unnorm_vecs <- norm_vecs %*% diag(scales)
unnorm_vecs <- unnorm_vecs * sign(colSums(res$scaling * unnorm_vecs))
unnorm_vecs

输出结果将与lda()的res$scaling一致。

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