Python：Generalized Discriminant Analysis (GDA) 手工代码实现 Kernel LDA

"""
Auther:Deniu He
Date:2021-03-16
"""
import numpy as np
from sklearn.datasets import load_iris
from sklearn.metrics.pairwise import rbf_kernel
import matplotlib.pyplot as plt
from scipy.linalg import block_diag

class GDA():
    def __init__(self,X,y,gamma):
        self.X = X
        self.y = y
        self.N = self.X.shape[0]    ## the number of training samples
        self.labels = np.unique(y)  ## the labels
        self.nClass = len(self.labels)  ## the number of labels(classes)
        self.nSample = []
        self.nAtt = self.X.shape[1]   ## the number of features(attributes)
        self.m = np.mean(self.X,axis=0)  ## the mean of all samples
        self.Mean = self.get_Mean()
        self.B = self.get_B()
        self.W = self.get_W()
        self.K = rbf_kernel(X=self.X,gamma=gamma)
        print("B::",self.B.shape)
        print("W::",self.W.shape)
        print("K::",self.K.shape)
        # self.L,self.U = np.linalg.eig(np.linalg.pinv(self.W @ self.K)@ (self.B @ self.K))
        self.L, self.U = np.linalg.eig(np.dot(np.linalg.inv(np.dot(self.W,self.K)),np.dot(self.B,self.K)))


    def get_Mean(self):
        Mean = np.zeros((self.nClass,self.nAtt))
        for i, lab in enumerate(self.labels):
            idx_list = np.where(self.y == lab)
            Mean[i] = np.mean(self.X[idx_list],axis=0)
            self.nSample.append(len(idx_list))
        return Mean

    def get_B(self):
        Diag = []
        for i, lab in enumerate(self.labels):
            idx_list = np.where(self.y == lab)
            nSub = np.size(idx_list)
            subMat = np.ones((nSub,nSub)) * (1/nSub)
            if i == 0:
                Diag = block_diag(subMat)
            else:
                Diag = block_diag(Diag,subMat)
        Diag = (1/self.N) * Diag
        B = Diag - (1/self.N)*np.ones((self.N,self.N))
        return B

    def get_W(self):
        A_M = []
        B_M = []
        for i, lab in enumerate(self.labels):
            idx_list = np.where(self.y == lab)
            nSub = np.size(idx_list)
            aMat = np.eye(nSub) * (1/nSub)
            bMat = np.ones((nSub,nSub)) * (1/self.N)
            if i == 0:
                A_M = block_diag(aMat)
                B_M = block_diag(bMat)
            else:
                A_M = block_diag(A_M,aMat)
                B_M = block_diag(B_M,bMat)
        W = A_M - B_M
        return W

    def DR(self,n_component):
        trans_X = []
        ord_ids = np.flipud(np.argsort(self.L))
        for i in range(n_component):
            tar_ids = ord_ids[i]
            u = self.U[:,tar_ids]
            if trans_X == []:
                 trans_X = self.K @ self.U[:,tar_ids] /np.sqrt(u@self.K@u)
            else:
                trans_X = np.vstack((trans_X,self.K @ self.U[:,tar_ids] /np.sqrt(u@self.K@u)))
        return trans_X.T


if __name__ == '__main__':
    X,y = load_iris(return_X_y=True)
    gamma = 10
    gda = GDA(X=X,y=y,gamma=gamma)
    print(gda.DR(n_component=2))

 有些问题：求出的特征值是复数

哪个大神帮忙修改一下