clc;
clear all;
close all;
%用mds对瑞士卷降维
% d=[0,411,213,219,296,397;...
% 411,0,204,203,120,152;...
% 213,204,0,73,136,245;...
% 219,203,73,0,90,191;...
% 296,120,136,90,0,109;...
% 397,152,245,191,109,0];
%瑞士卷的生成图
N=400;
t=(3*pi/2)*(1+2*rand(1,N));
s=21*rand(1,N);
X=[t.*cos(t);s;t.*sin(t)];
plot3(X(1,:),X(2,:),X(3,:),'.')
X=X';
[m,n]=size(X);
D=zeros(m,m);
for i=1:m
for j=i:m
D(i,j)=norm(X(i,:)-X(j,:));
D(j,i)=D(i,j);
end
end
d=D;
n=size(d,1);
B=zeros(n,n);
for i=1:n
for j=1:n
B(i,j)=-0.5*(d(i,j)^2 -1/n*d(i,:)*d(i,:)' -1/n*d(:,j)'*d(:,j) +1/n^2*sum(sum(d.^2)));
end
end
[V,D]=eig(B);
A=sort(diag(D),'descend');%把对角矩阵变成列向量
V1=[V(:,1),V(:,2)]*[A(1,1),0;0,A(2,1)];
figure;
plot(V1(:,1),V1(:,2),'.')