Why is the Scaling Matrix in LDA unnormalized?

Question

I was carrying out LDA (linear Discriminant Analysis) and noticed that the Scaling matrix produced by R is not normalized. Here is an example:

(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 

Normalizing the scaling matrix:

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 

Why is the scaling matrix not normalized?

Notice that if we try to fit lda manually:

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

Notice that the results Obtained by directly computing the scaling matrix is normalized. Of course the difference between the manually computed scaling matrix and the one produced by lda is the scale factor.

But this has an impact on the posterior probabilities.

Looking at the R code, i noticed they are doing svd twice rather than doing eigen decomposition. Tried analyzing the R code and after hours came up with the following:

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

Question: This is exactly like the unnormalized Scaling. But what is the intuition behind it?

Main QUESTION: in LDA Why is the scaling matrix not normalized? and How can I obtain this unnormalized scaling matrix using eigen decomposition? (Note that I am more interested with eigen decomposition since the eigenvectors are the solutions to the lda problem - But i do not mind having the SVD approach)

The lda objective: $$\max_{w} J(w) = \max_{w} \frac{w'S_{between}w}{w'S_{within}w}$$

EDIT:

I have found out that the formula used in R is the same used in python. Am I missing something regarding lda? Python too does give the unnormalized vectors.

if interested in the python code check here

Answer 1

Please study both my links in the comments. The 1st gives the LDA extraction algorithm formulas, and the 2nd (in my answer there) contains the demonstration of LDA on iris dataset according to such computations. I am not using R's term "scaling matrix", finding it a bit vague.

Per the algorithm expressed in my 1st link, LDA produces eigenvalues (which we recalculate into canonical correlations) and eigenvectors$^1$. Eigenvectors - we can multiply them by the constant $\sqrt{N-k}$ to obtain the (raw, unstandardized) discriminant coefficients $\bf C$ (called by you "unnormalized scaling matrix") which are the coefficients to compute discriminant scores with. Alternatively, we could normalize the eigenvectors to column SS=1, to obtain the matrix $\bf C_n$ (you call it "normalized scaling matrix") usable, for visual purposes, as the (oblique) rotation matrix of the variables into the discriminants.

So, you may use either way to get the values (scores) of the discriminants ($\bf X$ being the centered analyzed variables): one via discriminat coefficients: $\bf XC$, another via normalized eigenvectors: $\bf XC_n$. The first set of scores have the property that their pooled within-class covariance matrix is the identity matrix. The second set of scores are the perpendicular projections of the data points onto the discriminants as axes drawn in the space of the variables. By both sets of scores, the discriminants are uncorrelated variates. But, as drawn in the space of the original variables, they are not orthogonal axes.

By tradition, the first way to compute discr. scores is used most often by packages.

[As far as I'm aware, R and Python do not follow the fast algorithm I described under my link. They probably use equivalent, slower yet computationally more stable, variant utilizing SVD instead of eigendecomposition. Again, I'm mentioning that issue in my linked answer.]

---

$^1$ Because $\bf{S_w^{-1} S_b}$ of LDA isn't symmetric, its eigenvectors does not have to be an orthonormal matrix, unlike eigenvectors of PCA.