privacyPLS.py 文件源码

def run(X,y1,y2,dmax=None):

    D,N = X.shape
    y1_unique,J1 = np.unique(y1, return_inverse=True)
    ny1 = y1_unique.size

    y2_unique,J2 = np.unique(y2, return_inverse=True)
    ny2 = y2_unique.size

    Y1 = np.zeros((ny1,N))
    Y2 = np.zeros((ny2,N))
    Y1[J1,range(N)] = 1.
    Y2[J2,range(N)] = 1.

    XY2 = np.dot(X,Y2.T) # D x ny2
    XY2Y2X = np.dot(XY2,XY2.T) # D x D
    XX = np.dot(X,X.T) # D x D

    P = np.zeros((D,0))
    Sj = np.dot(X,Y1.T) #  D x ny1

    if dmax==None:
        dmax = D

    for d in range(dmax):
        if d>0: 
            invPP = np.linalg.pinv(np.dot(P.T,P))
            Sj -= np.dot(np.dot(np.dot(P,invPP),P.T),Sj) 

        C = np.dot(Sj,Sj.T) - XY2Y2X
        C = 0.5*(C+C.T)
        dd,E = scipy.linalg.eigh(C,eigvals=(D-1,D-1)) # ascending order

        assert np.isnan(dd).any()==False
        assert np.iscomplex(dd).any()==False
        #dd = dd[::-1] #
        #E = E[:,::-1]
        wj = E#E[:,0] # D x 1
        pj = np.dot(XX,wj) / np.dot(np.dot(wj.T,XX),wj) # D x 1
        P = np.hstack((P,pj.reshape((D,1)))) #  D x d


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