流形学习（一）LLE 在 MATLAB 中的实现及实例

最新推荐文章于 2024-08-09 08:36:15 发布

waitingwinter

最新推荐文章于 2024-08-09 08:36:15 发布

阅读量7.8k

点赞数 18

分类专栏：深度学习与机器学习文章标签： matlab 算法机器学习人工智能

本文链接：https://blog.csdn.net/waitingwinter/article/details/105467074

版权

深度学习与机器学习专栏收录该内容

8 篇文章 4 订阅

订阅专栏

LLE(Locally linear Embedding)

流形学习，全称流形学习方法(Manifold Learning)，自2000年在著名的科学杂志《Science》被首次提出以来，已成为信息科学领域的研究热点。在理论和应用上，流形学习方法都具有重要的研究意义。假设数据是均匀采样于一个高维欧氏空间中的低维流形，流形学习就是从高维采样数据中恢复低维流形结构，即找到高维空间中的低维流形，并求出相应的嵌入映射，以实现维数约简或者数据可视化。它是从观测到的现象中去寻找事物的本质，找到产生数据的内在规律。（百度百科）

下面我们将介绍三种流形学习及其在瑞士卷上的应用

LLE(Locally linear Embedding)

LLE 方法是发表在science上的文章，是流形学习的开山之作，自然被最先介绍。论文地址为

S. Roweis and L. Saul, “Nonlinear dimensionality reduction by locally linear embedding,” Science, vol. 290, pp. 2323–2326, 2000.

基本思想

给定数据集，通过最近邻等方式构造一个数据图(data graph)。然后在每一个局部区域，高维空间中的样本线性重构关系在低维空间中均得以保持。

Step1 最近邻确定与最佳权重

对于数据集中的任意一个 x_i , 按照欧式距离我们确定出离它最近的 k 个元素，记为 $x_{i_j}$ , 我们计算 x_i 的最优权重表示

由拉格朗日乘子法，可得 W_i 的线性表示解为

$W_i=\frac{(X_i^TX_i)^{-1}e}{e^T(X_i^TX_i)^{-1}e}$

其中

$e =(1,1,\cdots,1)^T\in \mathbb{R}^k,\\ X_i=(x_{i_1},x_{i_2},\cdots,x_{i_k})\in \mathbb{R}^{m\times k},\\ W_i=(w_{i,i_1},w_{i,i_2},\cdots,w_{i,i_k})\in \mathbb{R}^k.$

step2 全局嵌入学习模型

我们希望对于降维后的数据也有这种局部线性重构关系，即

$y_i \approx \sum_{j=1}^k w_{i,i_j}y_{i_j},$

故有此可得目标函数

$\sum_{i=1}^m\Vert y_i -\sum_{j=1}^k w_{i,i_j}y_{i_j}\Vert_2^2=tr(Y(I-W)^T(I-W)Y^T)$

故学习模型为

$\min\limits_{Y}\quad tr(Y(I-W)^T(I-W)Y^T)\quad s.t. \quad YY^T=I.$

step3 求解优化问题

对矩阵 (I-W)^T(I-W ) 进行特征分解

取出该矩阵最小的k+1个特征值对应的特征向量
丢弃特征值零对应的分量全相等的特征向量
即采用第2至第d+1个最小的特征值对应的特征向量组成样本的新的坐标

由于我们侧重于程序实现，故具体证明细节不作展开，只作简要分析，可以参考这篇博客查看具体细节

https://www.cnblogs.com/pinard/p/6266408.html?utm_source=itdadao&utm_medium=referral

算法步骤

MATLAB程序实现

考虑对如下瑞士卷数据进行分类

    N = 2000;
    rand('state',123456789);
    %Gaussian noise
    noise = 0.001*randn(1,N);
    %standard swiss roll data
    tt = (3*pi/2)*(1+2*rand(1,N));   height = 21*rand(1,N);
    X = [(tt+ noise).*cos(tt); height; (tt+ noise).*sin(tt)];

将其可视化

    %show the picture
    point_size = 20;
    figure(1)
    cla
    scatter3(X(1,:),X(2,:),X(3,:), point_size,tt,'filled');
    view([12 12]); grid off; axis off; hold on;
    axis off;
    axis equal;
    drawnow;

有

我们取最近邻个数为12 ，降维为2维，即

k = 12;
d = 2;
X = X';

计算每个元素对应的最近邻元素，并按照

来计算 W 的值，此处lambda是为了防止矩阵奇异

[m, ~] = size(X); % d < ~
lambda = 1e-10;
W = zeros(m); 
e = ones(k,1);
for i =1 : m
    xx = repmat(X(i, :), m, 1);
    diff = xx - X;
    dist = sum(diff.* diff, 2);
    [~, pos] = sort(dist);
    index = pos(1 : k + 1)';
    index(index == i) = [];
    w_numerator = (X(index, :) * X(index, :)' + lambda * eye(k)) \ e;
    w_denominator = e' * w_numerator;
    w = w_numerator / w_denominator;
    W(i, index) = w;
end

令 A =(I-W)^T(I-W) , 则问题的解为 A的第2到d+1个小特征值对应的特征向量，即

W =sparse(W);
I = eye(m);
A = (I - W)' * (I - W);
[eigenvector, eigenvalue] = eig(A);
eigenvalue = diag(eigenvalue);
[~,pos] = sort(eigenvalue);
index = pos(1: d+1);
tran = eigenvector(:, index);
p =sum(tran.*tran);
j = find(p == min(p));
tran(:, j) = [];
X = tran;

最后显示降维结果

figure(2)
scatter(X(:,1),X(:,2),point_size, tt,'filled');
grid off; axis off;

即

总的MATLAB代码

function X = LLE(X, k, d)
%% minifold learning: locally linear embedding
% k : neighbors
% x : data set 
% d : low dimension;
if nargin < 1
    N = 2000;
    rand('state',123456789);
    %Gaussian noise
    noise = 0.001*randn(1,N);
    %standard swiss roll data
    tt = (3*pi/2)*(1+2*rand(1,N));   height = 21*rand(1,N);
    X = [(tt+ noise).*cos(tt); height; (tt+ noise).*sin(tt)];
    
    %show the picture
    point_size = 20;
    figure(1)
    cla
    scatter3(X(1,:),X(2,:),X(3,:), point_size,tt,'filled');
    view([12 12]); grid off; axis off; hold on;
    axis off;
    axis equal;
    drawnow;
    X = X';
    k = 12;
    d = 2;
end
[m, ~] = size(X); % d < ~
lambda = 1e-10;
W = zeros(m); 
e = ones(k,1);
for i =1 : m
    xx = repmat(X(i, :), m, 1);
    diff = xx - X;
    dist = sum(diff.* diff, 2);
    [~, pos] = sort(dist);
    index = pos(1 : k + 1)';
    index(index == i) = [];
    w_numerator = (X(index, :) * X(index, :)' + lambda * eye(k)) \ e;
    w_denominator = e' * w_numerator;
    w = w_numerator / w_denominator;
    W(i, index) = w;
end
W =sparse(W);
I = eye(m);
A = (I - W)' * (I - W);
[eigenvector, eigenvalue] = eig(A);
eigenvalue = diag(eigenvalue);
[~,pos] = sort(eigenvalue);
index = pos(1: d+1);
tran = eigenvector(:, index);
p =sum(tran.*tran);
j = find(p == min(p));
tran(:, j) = [];
X = tran;
figure(2)
scatter(X(:,1),X(:,2),point_size, tt,'filled');
grid off; axis off;

下一篇我们将介绍 Isomap 算法的思想及实现，详见

https://blog.csdn.net/waitingwinter/article/details/105469511