Varying results when calculating scatter matrices for LDA

Question

I'm following a Linear Discriminant Analysis tutorial from here for dimensionality reduction. After working through the tutorial (did the PCA part, too), I shortened the code using sklearn modules where applicable and verified it on the Iris data set (same code, same result), a synthetic data set (with make_classification) and the sklearn-digits dataset.

However, then I tried the exact same code on a complete different (unfortunately non-public) data set that contains spectra recordings of two classes. The LDA crashes at the eigenvector verification part, where the $\lambda \mathbf{v}$ is supposed to be almost equal to $S_W^{-1} S_B \mathbf{v}$ (with $\lambda$ being the eigenvalue and $\mathbf{v}$ the corresponding eigenvector; $S_W$ and $S_B$ are the in/between-class scatter matrices). The first vector to be wrong seems at random positions, meaning each run it's a different vector that's causing this error.

I suspect it's related to rounding during calculations, since I get complex eigenvectors. For the PCA I just discarded the complex part (I think I read it somewhere in this forum), but this approach does not seem to work with LDA. Has anybody encountered similar problems or knows what's wrong?

Following is my code for the analysis, which is more or less the same as in the tutorial. I'm using the manual approach, since I'm interested in how many linear discriminants are needed to describe my data. (I'm not sure how to do this with sklearn's LDA.)

def LDAnalysis_manual(X, y):
    n_features = X.shape[1]
    n_classes = len(np.unique(y))

    print("Mean vectors...")
    mean_vectors = []
    for cl in range(n_classes):
        mean_vectors.append(np.mean(X[y == cl], axis=0))
        # print("Mean vector class {}: {}".format(cl, mean_vectors[cl - 1]))

    print("In-class scatter matrix...")
    S_W = np.zeros((n_features, n_features))
    for cl, mv in zip(range(1, n_classes), mean_vectors):
        class_sc_mat = np.zeros((n_features, n_features))  # each class' scatter matrix
        for row in X[y == cl]:
            row, mv = row.reshape(n_features, 1), mv.reshape(n_features, 1)  # column vectors
            class_sc_mat += (row - mv).dot((row - mv).T)
        S_W += class_sc_mat  # sum class scatter matrices

    overall_mean = np.mean(X, axis=0)
    print("Between-class scatter matrix...")
    S_B = np.zeros((n_features, n_features))

    for i, mean_vec in enumerate(mean_vectors):
        n = X[y == i + 1].shape[0]
        mean_vec = mean_vec.reshape(n_features, 1)  # make column vector
        overall_mean = overall_mean.reshape(n_features, 1)
        S_B += n * (mean_vec - overall_mean).dot((mean_vec - overall_mean).T)

    eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))

    print("Eigenvector test")
    for i in range(len(eig_vals)):
        print("\r{:3}".format(i), end=" ")
        sys.stdout.flush()

        eigv = eig_vecs[:, i].reshape(n_features, 1)
        np.testing.assert_array_almost_equal(np.linalg.inv(S_W).dot(S_B).dot(eigv).real,
                                             (eig_vals[i] * eigv).real,
                                             decimal=6, err_msg='', verbose=True)
    __log.debug("\nAll values ok.")

    eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:, i]) for i in
                 range(len(eig_vals))]  # make list of value & vector tuples
    eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)  # Sort tuple-list from high to low

    __log.info("\nEigenvalues (decending):")
    for i in eig_pairs:
        __log.info(i[0])

    tot = sum(eig_vals)
    var_exp = [(i / tot) for i in sorted(eig_vals, reverse=True)]
    cum_var_exp = np.cumsum(var_exp)
    cum_var_exp = cum_var_exp.real
    plot(len(var_exp), var_exp, cum_var_exp)

    idx_98 = next(idx for idx, val in enumerate(cum_var_exp) if val > .98)
    return idx_98 + 1

Answer 1

1. The LDA crashes for the exact reason you suspected. You have complex eigenvalues. If you use np.linalg.eigh, which was designed to decompose Hermetian matrices, you will always get real eigenvalues. np.linalg.eig can decompose nonsymetric square matrices, but, as you've suspected, it can produce complex eigenvalues. In short, np.linalg.eigh is more stable, and I would suggest using it for both PCA and LDA.
2. Dropping the complex part of the eigenvalues may have been acceptable in your specific example, but in practice it should be avoided. Depending on the size of the complex part of the number, it can significantly change the result. For example think of the multiplication of two complex conjugates. $(3+ .1i)*(3-.1i)=9-.01=9.01$ comared to $9$ when droping the complex part is a relatively safe, but $(3-2i)*(3+2i)=13$ compared to $9$ is a significant miscalculation. Using the above method for the eigen-decomposition will prevent this situation from arising.
3. Remember that one of the assumptions of LDA is that the features are normally distributed and independent of each other. Try running print('Class label distribution: %s' % np.bincount(y_train)[1:]). If the counts are not close to being equal, you've violated the first asumption of LDA, and the within-class scatter matrix must be scaled, in short divide each count by the number of class samples $N_i$. By doing this it should be obvious that computing the normalized scatter matrix is the same as computing the covariance matrix $\Sigma_i$. $$\Sigma_i=\frac{1}{N_i}S_W=\frac{1}{N_i}(x-m_i)(x-m_i)^T $$
4. Make sure your scaling your features before you do your PCA/LDA.
5. If the above doesn't fix your eigenvector verification step I suspect the problem is that the eigenvectors are scalled differently. Remember from your linear algebra class that a single eigenvalue, $\lambda_i$, has infinitely many eigenvectors, each being a scalar multiple of the others. $v_i=[1,2,3]$ and $v_i=[2,4,6]$ can both be the eigenvector of $\lambda_i$. So while you may get different values when calculating values at any given step after the decomposition, the end result should be the same.

Below is a template I use for LDA data compression. It assumes that you've split your data into a training and test set, the feature space has been properly scaled, and there are three classes in your label vector (you can adjust accordingly). It plots the individual and cumulative "discriminability" of each linear discriminant and then relies on the lda package in sklearn to transform the feature space using the number of discriminants you intend on using (here I chose to use the first 2). It also scales the within class scatter matrices by default.

LINEAR DISCRIMINANT ANALYSIS

calculate mean vectors

mean_vecs = []
for label in range(1, 4):
    mean_vecs.append(np.mean(X_train_std[y_train==label], axis=0))
    print('MV %s: %s\n' %(label, mean_vecs[label-1]))

calculate within-class scatter matrix

d = X_train_std.shape[1]
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
    class_scatter = np.cov(X_train_std[y_train==label].T)
    S_W += class_scatter
print('Scaled within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))

calculate between-class scatter matrix

mean_overall = np.mean(x_train_std, axis=0)
S_B = np.zeros((d, d))
for i, mean_vec in enumerate(mean_vec):
    n = X_train_std[y_train==i+1, :].shape[0]
    mean_vec = mean_vec.reshape(d, 1)
    mean_overall = mean_overall.reshape(d, 1)
    S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)
print('Between-class scatter matrix: %sx%s' % (S_B.shape[0], S_B.shape[1]))

eigen decomposition

eigen_vals, eigen_vecs = np.linalg.eigh(np.linalg.inv(S_W).dot(S_B))
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]
eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)
print('Eigendecomposition: \nEigenvalues in decreasing order:\n')
for eigen_val in eigen_pairs:
    print(eigen_val[0])

plot discriminablithy and select number of linear discriminants

tot = sum(eigen_vals.real)
discr = [(i / tot) for i in sorted(eigen_vals.real, reverse=True)]
cum_discr = np.cumsum(discr)

plt.bar(range(1, 14), discr, alpha=0.5, align='center',
        label='individual "discriminability"')
plt.step(range(1, 14), cum_discr, where='mid',
         label='cumulative "discriminability"')
plt.ylabel('"discriminability" ratio')
plt.xlabel('Linear Discriminants')
plt.ylim([-0.1, 1.1])
plt.legend(loc='best')
plt.tight_layout()
plt.show()

from sklearn.lda import LDA

lda = LDA(n_components=2)
x_train_lda = lda.fit_transform(X_train_std)
x_test_lda = lda.transform(X_test_std)
print('Features projected onto %d-dimensional LD subspace' % lda.n_components)

Answer 2

Note that using np.linalg.eigh will produce wrong results since np.linalg.inv(S_W).dot(S_B) is not Hermitian. It still should have real eigenvalues and eigenvectors, and nonzero imaginary parts are rounding errors. It should therefore be safe to just use the real parts, but you can ad a check if the imaginary parts are indeed small.