Varying results when calculating scatter matrices for LDA
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.
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.
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 $$
Make sure your scaling your features before you do your PCA/LDA.
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)