privacyLDA.py 文件源码

def run(X,y1,y2):
    # function [E,dd] = privacyLDA(X,y1,y2)
    # %max_W0 tr(W0'*C1*W0) - tr(W0'*C2*W0)
    D,N = X.shape
    y1_unique = np.unique(y1)
    #print y1_unique
    #[y1_unique,~,J1] = unique(y1);
    #ny1 = y1_unique.size

    y2_unique = np.unique(y2)
    #print y2_unique
    #[y2_unique,~,J2] = unique(y2);
    #ny2 = y2_unique.size

    C1 = np.zeros((D,D))
    C2 = np.zeros((D,D))
    mu = X.mean(axis=1).reshape((D,1))
    #print mu.shape

    for k in np.nditer(y1_unique):
        indk = np.where(y1==k)[0]        
        muk = X[:,indk].mean(axis=1).reshape((D,1))
        #muk -= np.kron(np.ones((1,len(indk))),mu)
        #%C1 = C1 + ny1*(muk-mu)*(muk-mu)';
        C1 = C1 + len(indk)*np.dot(muk-mu,(muk-mu).T)

    for k in np.nditer(y2_unique):
        indk = np.where(y2==k)[0]        
        muk = X[:,indk].mean(axis=1).reshape((D,1))
        #muk -= np.kron(np.ones((1,len(indk))),mu)
        #%C1 = C1 + ny1*(muk-mu)*(muk-mu)';
        C2 = C2 + len(indk)*np.dot(muk-mu,(muk-mu).T)

    C1 = C1 + 1e-8*np.trace(C1)*np.eye(D)# 
    C2 = C2 + 1e-8*np.trace(C2)*np.eye(D)#
    C1 = 0.5*(C1+C1.T)#;% + 1E-8*trace(C1)*eye(D); 
    C2 = 0.5*(C2+C2.T)#;% + 1E-8*trace(C2)*eye(D);


    dd,E = scipy.linalg.eigh(C1,C2) # ascending order

    #print dd.shape
    #print E.shape

    assert np.isnan(dd).any()==False
    assert np.iscomplex(dd).any()==False
    #[dd,ind] = sort(diag(dd),'descend'); 
    #print dd
    dd = dd[::-1] #
    E = E[:,::-1]

    E = E/np.tile(np.sqrt((E**2).sum(axis=0,keepdims=True)),(D,1))
    #% They need not be orthogonal.

    #print dd.shape
    #print E.shape

    return (E,dd)