流形学习——MDS & ISOMAP

最新推荐文章于 2024-06-23 15:10:21 发布

Vk_ch

最新推荐文章于 2024-06-23 15:10:21 发布

阅读量2.8k

点赞数 3

分类专栏： anomaly detection 文章标签： mds isomap 流形学习

本文链接：https://blog.csdn.net/inkky/article/details/81128650

版权

anomaly detection 专栏收录该内容

3 篇文章 0 订阅

订阅专栏

流形学习

流形学习

基本思想：在高维空间中，有低维结构。所以在局部具有欧式空间的性质。能用欧氏距离来进行距离计算。

如在三维空间中有二维平面

测地线距离geodesic：沿着地面的路线走，即b图红线的距离。
高维直线距离：虚线部分欧氏距离

所以高维中的测地线距离需要用近邻距离来近似。

1. 矩阵内积

矩阵的转置： $(AB)'=B'A'$
方阵的幂 $A^0=E$
由方阵A的元素所构成的行列式（各元素的位置不变），称为方阵A的行列式，记作|A|或 $detA$
$|\lambda A|=\lambda^n |A|$ （ $\lambda$ 是常数，A的阶数为n)
AB=BA=E。则我们称B是A的逆矩阵 $A^{-1}=B$
方阵不可逆，则行列式为0
若n阶方阵A满足 $A^TA$ =E （即 $A^{-1}$ = $A^T$ ），那称A为正交矩阵
正交变换刚好就是对数据做旋转、翻转操作的

先介绍向量的两种运算，一个行向量乘以一个列向量称作向量的内积，又叫作点积，结果是一个数；

一个列向量乘以一个行向量称作向量的外积，外积是一种特殊的克罗内克积，结果是一个矩阵，

2. 特征值，特征向量，特征分解

两篇文章参考
参考1
参考2

2.1 特征值和特征向量

红线为特征向量 v
A（变换矩阵）为特征向量矩阵
λ 是特征值
$\sum$ 为特征值的对角矩阵
I 为单位矩阵

1、在这种情况下变换仅仅是水平方向乘以因子2和垂直方向乘以因子0.5，使得变换矩阵A定义为

2、在一般情况下，特征向量矩阵A的特征向量 v 满足下列式子

A v = λ v

$Av=\lambda v$
λ为特征向量 v 对应的特征值。我们可以将上式重写为

3、如果一个方阵是不可逆的，这意味着它的行列式必须等于零。因此，要找到A的特征向量，我们只需要解决以下公式：

4、利用判别式 $\Delta=b^2-4ac$ 来找到 λ 的值。需要注意的是大小为NxN的方阵总是具有N个特征值，每一个对应一个特征向量。特征值指定特征向量的大小。

5、将 λ 代入公式8，找到特征向量。

2.2 特征分解

局限：变换的矩阵必须是方阵

特征值分解是将一个矩阵分解为如下形式：

其中，Q是这个矩阵A的特征向量组成的矩阵，Σ是一个对角矩阵，每一个对角线元素就是一个特征值，里面的特征值是由大到小排列的，这些特征值所对应的特征向量就是描述这个矩阵变化方向（从主要的变化到次要的变化排列）。也就是说矩阵A的信息可以由其特征值和特征向量表示。

对于矩阵为高维的情况下，那么这个矩阵就是高维空间下的一个线性变换。可以想象，这个变换也同样有很多的变换方向，我们通过特征值分解得到的前N个特征向量，那么就对应了这个矩阵最主要的N个变化方向。我们利用这前N个变化方向，就可以近似这个矩阵（变换）。

3. 数据降维

3.1 MDS（Multiple Dimensional Scaling）

参考

我们有n个样本，每个样本维度为m。我们的目标是用不同的新的k维向量 ( $k<<m$ ) 替代原来的n个m维向量，使得在新的低维空间中，所有样本相互之间的距离等于（或最大程度接近）原空间中的距离（默认欧氏距离）。

MDS接收的输入是一个距离矩阵D（已经降维后的矩阵）

不能直接转换到原始数据的原因

对点做平移、旋转等变化时，点之间的距离时不变的

3.1.1 实现流程

思路是：从D->B->X

假设： X是一个n×q的矩阵，n为样本数，q是原始的维度；降维后的距离矩阵D

定义：内积矩阵B

B = X X T (n \times n) = (X M) (X M) T (M 是 一 组 正 交 基 ， 对 X 做 正 交 变 换) = X M M T X = X X T (7)

$\begin{equation} \begin{aligned} B &= XX^T \space\space\space\space(n\times n) \\ &= (XM)(XM)^T \space\space\space\space(M是一组正交基，对X做正交变换)\\ &= XMM^TX \\ &= XX^T \end{aligned} \end{equation}$
所以如果我们想通过B反算出X，肯定是没法得到真正的XX, 而是它的任意一种正交变换后的结果。

B的每个元素：

b i j = \sum k = 1 q x i k x j k ， i 是 行 ， j 是 列

$b_{ij} = \sum_{k=1}^{q}x_{ik}x_{jk}， i 是行，j是列$
D的每个元素：

d 2 i j = (x i - x j) 2 = \sum k = 1 q (x i k - x j k) 2 (求 每 个 向 量 的 距 离) = \sum k = 1 q x 2 i k + x 2 j k - 2 x i k x j k = b i i + b j j - 2 b i j (8)

$\begin{equation} \begin{aligned} d_{ij}^2 &= (x_i-x_j)^2 \\ &= \sum_{k=1}^{q}(x_{ik}-x_{jk})^2 \space\space\space\space(求每个向量的距离)\\ &= \sum_{k=1}^{q}x_{ik}^2+x_{jk}^2-2x_{ik}x_{jk}\\ &=b_{ii}+b_{jj}-2b_{ij} \end{aligned} \end{equation}$
根据公式3，我们可以得到

b i j = - 1 2 (d 2 i j - b i i - b j j)

$b_{ij} = -\frac12(d_{ij}^2-b_{ii}-b_{jj})$
trick去中心化：

数据的中心点平移到原点 $\sum_{i=1}^nx_{ik} = 0, for \space all \space k =1..q$

即B的任意行(row)之和以及任意列(column)之和都为0了。所以我们可以得到

b i i b j j 2 T n = T n + 1 n \sum j = 1 n d 2 i j = T n + 1 n \sum i = 1 n d 2 i j = 1 n 2 \sum i = 1 n \sum j = 1 n d 2 i j (9) (10) (11)

$\begin{align} b_{ii} &= \frac{T}n+\frac1n\sum_{j=1}^nd_{ij}^2\\ b_{jj} &= \frac{T}n+\frac1n\sum_{i=1}^nd_{ij}^2\\ \frac{2T}{n} &= \frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^nd_{ij}^2 \end{align}$
特别重要的公式（在代码中应用）

b i j = - 1 2 (d 2 i j - b i i - b j j) = - 1 2 (d 2 i j - 1 n \sum j = 1 n d 2 i j - 1 n \sum i = 1 n d 2 i j + 2 T n) = - 1 2 (d 2 i j - 1 n \sum j = 1 n d 2 i j - 1 n \sum i = 1 n d 2 i j + 1 n 2 \sum i = 1 n \sum j = 1 n d 2 i j) = - 1 2 (d 2 i j - d 2 i \cdot - d 2 \cdot j + d 2 \cdot \cdot) (12)

$\begin{equation} \begin{aligned} b_{ij} &= -\frac12(d_{ij}^2-b_{ii}-b_{jj})\\ &= -\frac12(d_{ij}^2-\frac1n\sum_{j=1}^nd_{ij}^2-\frac1n\sum_{i=1}^nd_{ij}^2+\frac{2T}{n})\\ &= -\frac12(d_{ij}^2-\frac1n\sum_{j=1}^nd_{ij}^2-\frac1n\sum_{i=1}^nd_{ij}^2+\frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^nd_{ij}^2)\\ &= -\frac12(d_{ij}^2-d_{i\cdot}^2-d_{\cdot j}^2+d_{\cdot\cdot}^2) \end{aligned} \end{equation}$
可以看到

d2i⋅ d i ⋅ 2 $d_{i\cdot}^2$ 是

D2 D 2 $D^2$ 行均值；

d2⋅j d ⋅ j 2 $d_{\cdot j}^2$ 是列均值；

d2⋅⋅ d ⋅ ⋅ 2 $d_{\cdot\cdot}^2$ 是矩阵的均值。

这样我们就可以通过矩阵D得到矩阵B了。因为B是对称的矩阵，所以可以通过特征分解得到：

B = V Λ V - 1 = V Λ V T

$\begin{aligned} B &= V\Lambda V^{-1}\\ &= V\Lambda V^T \end{aligned}$
D会是一个对称实矩阵，此时得到的B刚好会有q个非0的特征值，也就是说B的秩等于q，如果我们想还原X，就选择前q个特征值和特征向量；如果想要达到降维的目的，就选择制定的p个(p

3.1.2 pseudocode

3.1.3 code(python)

def mds(D, q):
    """
    将 D 转化成内积矩阵 B：
    bij=-0.5(dij^2-di.^2-d.j^2+d..^2)
    将 B 转化成 X：
    X=特征向量*sqrt(特征值)

    :param D:距离矩阵
    :param q:降维的目标
    :return:降维后的目标矩阵
    """
    D = np.asarray(D)
    DSquare = D ** 2  # d^2
    totalMean = np.mean(DSquare)  # d..^2
    columnMean = np.mean(DSquare, axis=0)  # d.j^2
    rowMean = np.mean(DSquare, axis=1)  # di.^2
    B = np.zeros(DSquare.shape)
    for i in range(B.shape[0]):
        for j in range(B.shape[1]):
            B[i][j] = -0.5 * (DSquare[i][j] - rowMean[i] - columnMean[j] + totalMean)
    eigVal, eigVec = np.linalg.eig(B)  # 找到特征值和特征向量
    print(eigVal, eigVec)
    X = np.dot(eigVec[:, :q], np.sqrt(np.diag(eigVal[:q])))

    return X

3.2 isomap（等度量映射）

A Global Geometric Framework for Nonlinear Dimensionality Reduction. Joshua B. Tenenbaum, Vin de Silva, John C. Langford. 2000.

3.2.1 核心

构造点之间的距离

对于每个点，基于欧氏距离找出其邻近点，（邻近点之间存在连接，非邻近点不存在连接），建立邻近连接图。

3.2.2 过程

1.为每个数据点确定邻居 ，将它们连接起来构造带权邻接图G（复杂度： $O(DN^2)$ ）

1.1 通过kNN(k-Nearest Neighbor)找到点的k个最近邻

1.2 把半径ϵ内的所有点作为邻居

2.通过计算同中各点之间的最短路径Dijkstra ，作为点之间的距离dij放入距离矩阵D（复杂度： $O(DN^2)$ ）

3.将D传给经典的MDS算法，得到降维后的结果。（复杂度： $O(qN^2)$ ）

3.2.3 pesudocode

# -*- coding: utf-8 -*-
# @Time    : 2018/7/11 14:37
# @Author  : Inkky
# @Email   : yingyang_chen@163.com
'''
isomap算法主要流程
1、找k近邻的点形成近邻图
2、找最短路径（Floyd或者Dijkstra）形成近邻矩阵
3、调用MDS降维

N ：训练数据点数
D： 近邻矩阵
k ：最近邻居数
q ：输出尺寸
'''

import numpy as np
from sklearn import datasets
from MDS import mds


import matplotlib.pyplot as plt


def distancematirx(data):
    n = len(data)
    dis_matrix = np.zeros([n, n], np.float32)
    for i in range(n):
        for j in range(n):
            dis_matrix[i, j] = np.linalg.norm(data[i] - data[j])
    return dis_matrix


def dijkstra(data, start):
    n = len(data)
    col = data[start].copy()  # 一维数组 dis 来存储 1 号顶点到其余各个顶点的初始路程
    for i in range(n):
        tmp = np.argpartition(col, 1)[1]  # 找到第二小的值，最小为自身
        dis = data[start, tmp]  # 找到最短的长度
        for j in range(n):  # 更新距离矩阵
            if data[start, j] > dis + data[i, j]:
                data[start, j] = dis + data[i, j]
                data[j, start] = data[start, j]
        col[tmp] = float('inf')
    return data


def isomap(rawdata, target, k):
    inf = float('inf')
    n = len(rawdata)
    if k >= n:
        raise ValueError('K is too large')
    mat_matrix = distancematirx(rawdata)  # 计算原来矩阵的距离
    print(mat_matrix)

    # 找到k近邻
    k_near = np.ones([n, n], np.float32) * inf
    for i in range(n):
        temp = np.argpartition(mat_matrix[i], k)[:k + 1]
        print(temp)
        k_near[i][temp] = mat_matrix[i][temp]
    print('找到k近邻\n', k_near)

    # 调用最短路径短发计算任意两样本点之间的距离
    for i in range(n):
        dist = dijkstra(k_near, i)
    print('调用最短路径短发计算任意两样本点之间的距离\n', dist)

    # 调用mds降维
    return mds(dist, target)


if __name__ == '__main__':
    digits = datasets.load_digits(n_class=4)
    X = digits.data[:10,:10] #只能为偶数
    y = digits.target
    n_samples, n_features = X.shape
    print(n_samples, n_features)


    print('开始降维.....')
    # D = np.array([[1, 2, 3, 4], [2, 1, 5, 6], [3, 5, 1, 7], [4, 6, 7, 1], [3, 4, 5, 6], [6, 1, 2, 3]])  # test data
    outcome = isomap(X, 2, 4)
    print('降维完成\n')
    print(outcome)