# Why is Python's scikit-learn LDA not working correctly and how does it compute LDA via SVD?

I was using the Linear Discriminant Analysis (LDA) from the scikit-learn machine learning library (Python) for dimensionality reduction and was a little bit curious about the results. I am wondering now what the LDA in scikit-learn is doing so that the results look different from, e.g., a manual approach or an LDA done in R. It would be great if someone could give me some insights here.

What's basically most concerning is that the scikit-plot shows a correlation between the two variables where there should be a correlation 0.

For a test, I used the Iris dataset and the first 2 linear discriminants looked like this:

IMG-1. LDA via scikit-learn

This is basically consistent with the results I found in the scikit-learn documentation here.

Now, I went through the LDA step by step and got a different projection. I tried different approaches in order to find out what was going on:

IMG-2. LDA on raw data (no centering, no standardization)

And here would be the step-by-step approach if I standardized (z-score normalization; unit variance) the data first. I did the same thing with mean-centering only, which should lead to the same relative projection image (and which it indeed did).

IMG-3. Step-by-step LDA after mean-centering, or standardization

IMG-4. LDA in R (default settings)

LDA in IMG-3 where I centered the data (which would be the preferred approach) looks also exactly the same as the one that I found in a Post by someone who did the LDA in R

## Code for reference

I did not want to paste all the code here, but I have uploaded it as an IPython notebook here broken down into the several steps I used (see below) for the LDA projection.

1. Step 1: Computing the d-dimensional mean vectors $$\mathbf m_i = \frac{1}{n_i} \sum\limits_{\mathbf x \in D_i}^n \; \mathbf x_k$$
2. Step 2: Computing the Scatter Matrices

2.1 The within-class scatter matrix $S_W$ is computed by the following equation:
$$S_W = \sum\limits_{i=1}^{c} S_i = \sum\limits_{i=1}^{c} \sum\limits_{\mathbf x \in D_i}^n (\mathbf x - \mathbf m_i)\;(\mathbf x - \mathbf m_i)^T$$

2.2 The between-class scatter matrix $S_B$ is computed by the following equation:
$$S_B = \sum\limits_{i=1}^{c} n_i (\mathbf m_i - \mathbf m) (\mathbf m_i - \mathbf m)^T$$ where $\mathbf m$ is the overall mean.

3. Step 3. Solving the generalized eigenvalue problem for the matrix $S_{W}^{-1}S_B$

3.1. Sorting the eigenvectors by decreasing eigenvalues

3.2. Choosing k eigenvectors with the largest eigenvalues. Combining the two eigenvectors with the highest eigenvalues to construct our $d \times k$-dimensional eigenvector matrix $\mathbf W$

4. Step 5: Transforming the samples onto the new subspace $$\mathbf y = \mathbf W^T \times \mathbf x.$$

• I haven't gone through to look for the differences, but you can see exactly what scikit-learn is doing in the source. – Dougal Jul 28 '14 at 19:52
• It looks like they are also standardizing (centering and then scaling via division by the standard deviation). This, I would expect a result similar to the one in my 3rd plot (and the R) plot...hmm – user39663 Jul 28 '14 at 20:15
• Weird: the plot you obtained with scikit (and the one they show in their documentation) does not make sense. LDA always yields projections that have correlation zero, but obviously there is a very strong correlation between scikit's projections on discriminant axes 1 and 2. Something is clearly wrong there. – amoeba says Reinstate Monica Jul 28 '14 at 21:52
• @ameoba Yes, I think so too. What's also weird is that the same plot I am showing for scikit is in the example documentation: scikit-learn.org/stable/auto_examples/decomposition/… That makes me think that my usage of scikit it correct, but that there is something odd about the LDA function – user39663 Jul 28 '14 at 21:59
• @SebastianRaschka: Yes, I noticed. It is weird indeed. However, notice that the first of your own (non-scikit) LDA plots also shows non-zero correlation and hence something must be wrong with it as well. Did you centre the data? Projection on the second axis does not seem to have zero mean. – amoeba says Reinstate Monica Jul 28 '14 at 23:04

Update: Thanks to this discussion, scikit-learn was updated and works correctly now. Its LDA source code can be found here. The original issue was due to a minor bug (see this github discussion) and my answer was actually not pointing at it correctly (apologies for any confusion caused). As all of that does not matter anymore (bug is fixed), I edited my answer to focus on how LDA can be solved via SVD, which is the default algorithm in scikit-learn.

After defining within- and between-class scatter matrices $\boldsymbol \Sigma_W$ and $\boldsymbol \Sigma_B$, the standard LDA calculation, as pointed out in your question, is to take eigenvectors of $\boldsymbol \Sigma_W^{-1} \boldsymbol \Sigma_B$ as discriminant axes (see e.g. here). The same axes, however, can be computed in a slightly different way, exploiting a whitening matrix:

1. Compute $\boldsymbol \Sigma_W^{-1/2}$. This is a whitening transformation with respect to the pooled within-class covariance (see my linked answer for details).

Note that if you have eigen-decomposition $\boldsymbol \Sigma_W = \mathbf{U}\mathbf{S}\mathbf{U}^\top$, then $\boldsymbol \Sigma_W^{-1/2}=\mathbf{U}\mathbf{S}^{-1/2}\mathbf{U}^\top$. Note also that one compute the same by doing SVD of pooled within-class data: $\mathbf{X}_W = \mathbf{U} \mathbf{L} \mathbf{V}^\top \Rightarrow \boldsymbol\Sigma_W^{-1/2}=\mathbf{U}\mathbf{L}^{-1}\mathbf{U}^\top$.

2. Find eigenvectors of $\boldsymbol \Sigma_W^{-1/2} \boldsymbol \Sigma_B \boldsymbol \Sigma_W^{-1/2}$, let us call them $\mathbf{A}^*$.

Again, note that one can compute it by doing SVD of between-class data $\mathbf{X}_B$, transformed with $\boldsymbol \Sigma_W^{-1/2}$, i.e. between-class data whitened with respect to the within-class covariance.

3. The discriminant axes $\mathbf A$ will be given by $\boldsymbol \Sigma_W^{-1/2} \mathbf{A}^*$, i.e. by the principal axes of transformed data, transformed again.

Indeed, if $\mathbf a^*$ is an eigenvector of the above matrix, then $$\boldsymbol \Sigma_W^{-1/2} \boldsymbol \Sigma_B \boldsymbol \Sigma_W^{-1/2}\mathbf a^* = \lambda \mathbf a^*,$$ and multiplying from the left by $\boldsymbol \Sigma_W^{-1/2}$ and defining $\mathbf a = \boldsymbol \Sigma_W^{-1/2}\mathbf a^*$, we immediately obtain: $$\boldsymbol \Sigma_W^{-1} \boldsymbol \Sigma_B \mathbf a = \lambda \mathbf a.$$

In summary, LDA is equivalent to whitening the matrix of class means with respect to within-class covariance, doing PCA on the class means, and back-transforming the resulting principal axes into the original (unwhitened) space.

This is pointed out e.g. in The Elements of Statistical Learning, section 4.3.3. In scikit-learn this is the default way to compute LDA because SVD of a data matrix is numerically more stable than eigen-decomposition of its covariance matrix.

Note that one can use any whitening transformation instead of $\boldsymbol \Sigma_W^{-1/2}$ and everything will still work exactly the same. In scikit-learn $\mathbf{L}^{-1}\mathbf{U}^\top$ is used (instead of $\mathbf{U}\mathbf{L}^{-1}\mathbf{U}^\top$), and it works just fine (contrary to what was originally written in my answer).

• Thanks for this nice answer. I appreciate it that you took the time to write it up that nicely. Maybe you could mention it in the discussion on GitHub; I am sure that would be of help to fix the LDA in the next version of sci-kit – user39663 Aug 5 '14 at 15:07
• @SebastianRaschka: I don't have an account on GitHub. But if you want, you can give there a link to this thread. – amoeba says Reinstate Monica Aug 5 '14 at 15:17
• @amoeba: Textbooks usually describe the LDA as you did - an eigenvalue decomposition of $\boldsymbol \Sigma_W^{-1} \boldsymbol \Sigma_B$. Curiously, a number of LDA implementations I know take a different approach. Their axes are the vectors to the class means transformed with $\boldsymbol \Sigma_W^{-1}$. Your LDA solution is an orthonormal basis of these vectors. Scikit-learn's LDA gives the same results as these implementations, so I don't think there is actually an error. – kazemakase Aug 6 '14 at 7:38
• For reference, here are the implementations I was talking about: sourceforge.net/p/mlpy/code/ci/default/tree/mlpy/da.py#l24 github.com/sccn/BCILAB/blob/master/code/machine_learning/… mathworks.com/matlabcentral/fileexchange/… – kazemakase Aug 6 '14 at 7:41
• @kazemakase: Well, of course if there are only two classes, then $\boldsymbol \Sigma_B$ has rank 1, and everything simplifies a lot, as the only eigenvector of $\boldsymbol \Sigma_W^{-1}\boldsymbol \Sigma_B$ is given by $\boldsymbol \Sigma_W^{-1}(\boldsymbol\mu_1 - \boldsymbol\mu_2)$, where $\boldsymbol\mu_i$ are class means. I guess that is what you meant before? This is nicely covered e.g. in Bishop's ML textbook, section 4.1.4. But the generalization to more classes requires eigen-analysis (Ibid., 4.1.6). Also, scikit's code (that we are discussing here!) does use svd, twice actually. – amoeba says Reinstate Monica Aug 6 '14 at 13:45

Just to close this question, the discussed issue with the LDA has been fixed in scikit-learn 0.15.2.