局部保留投影(Locality Preserving Projection, LPP)

Locality Preserving Projections.

paper：Locality Preserving Projection

局部保留投影(Locality Preserving Projection, LPP)算法是一种线性降维方法。对于高维空间中的样本$X \in \Bbb{R}^{N \times D}$，LPP寻找一个投影矩阵$A \in \Bbb{R}^{D \times d}$，从而构造低维空间的降维结果$Y=XA \in \Bbb{R}^{N \times d}$。

LPP在降维时考虑了样本点之间的相对位置关系。构造一个关系矩阵$W \in \Bbb{R}^{N \times N}$衡量任意两个样本点之间的距离。关系矩阵$W$构造如下：

\[W_{ij} = \begin{cases} e^{-\frac{||x_i-x_j||^2}{t}}, & j\neq i, j \in \text{KNN}(i) \\ 0, & \text{otherwise} \end{cases}\]

直观地，当点$j$与点$i$距离比较近时，使用径向基函数衡量两者的相对位置关系；并且对点$i$只考虑最接近的$k$个点(使用kNN实现)。定义$W_{ii}=0$是为了后续推导方便。

关系矩阵$W$的构造过程为：

def cal_pairwise_dist(X):
    #返回任意两个点之间的欧式距离
    N,D = np.shape(X)
    tile_xi = np.tile(np.expand_dims(X,1),[1,N,1]) # [N,N,D]
    tile_xj = np.tile(np.expand_dims(X,axis=0),[N,1,1]) # [N,N,D]
    dist = np.sum((tile_xi-tile_xj)**2,axis=-1)  # [N,N]
    return dist
    
def rbf(dist, t = 1.0):
    #径向基函数
    return np.exp(-(dist/t))

def cal_rbf_dist(data, n_neighbors = 10, t = 1):
    #计算关系矩阵
    dist = cal_pairwise_dist(data)
    N = dist.shape[0]
    rbf_dist = rbf(dist, t)
    W = np.zeros([N, N])
    for i in range(N):
        #跳过与样本点与自身的关系
        index_ = np.argsort(dist[i])[1:1 + n_neighbors]
        W[i, index_] = rbf_dist[i, index_]
        W[index_, i] = rbf_dist[index_, i]
    return W

LPP的优化目标函数如下：

\[\sum_{i,j} (y_i-y_j)^2W_{ij}\]

如果样本点$j$与点$i$的距离比较接近，则$W_{ij}$比较大，通过最小化上述目标使得降维后$y_j$与$y_i$比较接近。如果样本点$j$与点$i$的距离比较远，则$W_{ij}=0$，不再限制降维后$y_j$与$y_i$的关系。

目标函数可进一步写作：

\[\begin{aligned} \sum_{i,j} (y_i-y_j)^2W_{ij} &= \sum_{i,j} (y_i^TW_{ij}y_i+y_j^TW_{ij}y_j-2y_i^TW_{ij}y_j) \\ &= 2\sum_{i,j} y_i^TW_{ij}y_i-2\sum_{i,j}y_i^TW_{ij}y_j \\ &= 2\sum_{i} y_i^T(\sum_{j}W_{ij})y_i-2\sum_{i,j}y_i^TW_{ij}y_j \\ &(\text{记}D_{ii} = \sum_{j}W_{ij}) \\ &= 2(\sum_{i} y_i^TD_{ii}y_i-\sum_{i,j}y_i^TW_{ij}y_j) \end{aligned}\]

记$D=\text{diag}(D_{ii})\in \Bbb{R}^{N \times N}$。首先考虑降维维度是$1$的情况，即$A=a\in \Bbb{R}^{D \times 1}$，则目标函数写作矩阵形式：

\[Y^T(D-W)Y =Y^TLY= a^TX^TLXa\]

其中$L=D-W\in \Bbb{R}^{N \times N}$是样本的拉普拉斯矩阵。

引入约束$Y^TDY=a^TX^TDXa=1$，通过将$Y^TDY$固定为常数，使得当样本的相对关系比较接近时($D$比较接近)，降维结果也比较接近($Y$比较接近)。

则LPP的优化问题写作：

\[\begin{aligned} \mathop{\min}_a \quad & a^TX^TLXa \\ \text{s.t.} \quad & a^TX^TDXa=1 \end{aligned}\]

构造拉格朗日函数：

\[L=a^TX^TLXa + \lambda(1-a^TX^TDXa)\]

对上式求极值，得：

\[2X^TLXa -2\lambda X^TDXa = 0\]

整理得：

\[(X^TDX)^{-1}X^TLXa = \lambda a\]

因此$a$为矩阵$(X^TDX)^{-1}X^TLX$的最小特征值对应的特征向量。

一般地，当$d>1$时，降维矩阵$A$为矩阵$(X^TDX)^{-1}X^TLX$的$d$个最小的特征值对应的特征向量组成的矩阵。

LPP的实现过程如下：

def lpp(X, n_dims = 2, n_neighbors = 30, t = 1.0):
    N = X.shape[0]
    W = cal_rbf_dist(X, n_neighbors, t)

    #计算矩阵D和L
    D = np.zeros_like(W)
    for i in range(N):
        D[i,i] = np.sum(W[i])
    L = D - W

    #计算矩阵特征值
    XTDX = np.dot(np.dot(X.T, D), X)
    XTLX = np.dot(np.dot(X.T, L), X)
    eig_val, eig_vec = np.linalg.eig(np.dot(np.linalg.pinv(XDXT), XLXT))

    #特征值从小到大排序
    sort_index_ = np.argsort(np.abs(eig_val))
    eig_val = eig_val[sort_index_]

    #舍弃太小的特征值
    j = 0
    while eig_val[j] < 1e-6:
        j+=1
    sort_index_ = sort_index_[j:j+n_dims]
    eig_val_picked = eig_val[j:j+n_dims]

    #根据特征值对应的特征向量构造降维矩阵A
    A = eig_vec[:, sort_index_]

    #降维
    Y = np.dot(X, A)
    return Y