计算 LDA 的散布矩阵时的不同结果

数据挖掘 Python scikit-学习判别分析

2021-10-12 19:37:59

我正在按照这里的线性判别分析教程进行降维。在完成教程（也做了 PCA 部分）之后，我在适用的情况下使用 sklearn 模块缩短了代码，并在 Iris 数据集（相同的代码，相同的结果）、合成数据集（使用make_classification）和 sklearn-数字数据集。

但是，然后我在包含两类光谱记录的完全不同（不幸的是非公开）数据集上尝试了完全相同的代码。LDA 在特征向量验证部分崩溃，其中 $\lambda \mathbf{v}$ 应该几乎等于 $S_W^{-1} S_B \mathbf{v}$ （和 $\lambda$ 是特征值和 $\mathbf{v}$ 对应的特征向量； $S_W$ 和 $S_B$ 是类内/类间散布矩阵）。第一个错误的向量似乎在随机位置，这意味着每次运行都是导致此错误的不同向量。

我怀疑它与计算期间的舍入有关，因为我得到了复杂的特征向量。对于 PCA，我只是丢弃了复杂的部分（我想我在这个论坛的某个地方读过它），但这种方法似乎不适用于 LDA。有没有人遇到过类似的问题或知道出了什么问题？

以下是我的分析代码，与教程中的大致相同。我正在使用手动方法，因为我对需要多少线性判别式来描述我的数据感兴趣。（我不确定如何使用 sklearn 的 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

2个回答

LDA 因您怀疑的确切原因而崩溃。你有复杂的特征值。如果您使用np.linalg.eigh旨在分解 Hermetian 矩阵的，您将始终获得真正的特征值。 np.linalg.eig可以分解非对称方阵，但是，正如您所怀疑的，它可以产生复杂的特征值。简而言之，np.linalg.eigh它更稳定，我建议将它用于 PCA 和 LDA。
在您的具体示例中，删除特征值的复杂部分可能是可以接受的，但在实践中应该避免。根据复数部分的大小，它可以显着改变结果。例如考虑两个复共轭的乘法。 $(3+ .1i)*(3-.1i)=9-.01=9.01$ 来到 $9$ 当丢弃复杂的部分是相对安全的，但是 $(3-2i)*(3+2i)=13$ 相比 $9$ 是一个严重的误判。使用上述方法进行特征分解将防止出现这种情况。
请记住，LDA 的假设之一是特征是正态分布的并且彼此独立。尝试运行print('Class label distribution: %s' % np.bincount(y_train)[1:])。如果计数不接近相等，则违反了 LDA 的第一个假设，并且必须缩放类内散布矩阵，简而言之，将每个计数除以类样本的数量 $N_i$ . 通过这样做，很明显计算归一化散布矩阵与计算协方差矩阵相同 $\Sigma_i$ . $Σ_{一世} = \frac{1}{ñ_{一世}} {小号}_{W} = \frac{1}{ñ_{一世}} (X - 米_{一世}) (X - 米_{一世})^{吨}$ $\Sigma_i=\frac{1}{N_i}S_W=\frac{1}{N_i}(x-m_i)(x-m_i)^T$
确保在进行 PCA/LDA 之前缩放功能。
如果上述方法不能解决您的特征向量验证步骤，我怀疑问题在于特征向量的命名方式不同。从你的线性代数课中记住一个单一的特征值， $\lambda_i$ ，有无限多个特征向量，每个都是其他的标量倍数。 $v_i=[1,2,3]$ 和 $v_i=[2,4,6]$ 都可以是特征向量 $\lambda_i$ . 因此，虽然在分解后的任何给定步骤计算值时可能会得到不同的值，但最终结果应该是相同的。

下面是我用于 LDA 数据压缩的模板。它假设您已将数据拆分为训练集和测试集，特征空间已正确缩放，并且标签向量中有三个类（您可以进行相应调整）。它绘制了每个线性判别式的个体和累积“可辨别性”，然后依靠lda包sklearn使用您打算使用的判别式数量来转换特征空间（这里我选择使用前 2 个）。默认情况下，它还缩放类内散布矩阵。

线性判别分析

计算平均向量

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]))

计算类内散布矩阵

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]))

计算类间散布矩阵

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_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])

绘制可辨别性并选择线性判别式的数量

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)

请注意，使用 np.linalg.eigh 会产生错误的结果，因为 np.linalg.inv(S_W).dot(S_B) 不是 Hermitian。它仍然应该有实特征值和特征向量，非零虚部是舍入误差。因此，只使用实部应该是安全的，但是您可以检查虚部是否确实很小。

其它你可能感兴趣的问题

上一篇我应该信任哪个功能选择？下一篇在 R 中对具有多个级别的列进行虚拟编码