LTSA 改進版本，和經典LTSA類似

function Y = iterative_LTSA(data,d,K,max_iterate_time,errors)
[D,N] = size(data);  % D是輸入樣本點的維數. N是輸入樣本個數
indexs=zeros(N,K);

% Step 0:  Neighborhood Index（k-d method）
%if nargin<4                   %如果輸入參數個數<4

    NI = cell(1,N);            %産生一個1*N的空單元數組
    if D>N
        a = sum(data.*data);  %  .* 用于矩陣中對應元素的相乘,并且豎向相加
        dist2 = sqrt(repmat(a',[1 N]) + repmat(a,[N 1]) - 2*(data'*data)); %計算距離
        for i=1:N
            % Determine ki nearest neighbors of x_j
            [dist_sort,J] = sort(dist2(:,i));   %sort 對dist2 進行排序，dist_sort 傳回對每列排序後的結果
            %J傳回的是排序後矩陣中每個數在矩陣未排序前每列中的位置（即索引值）
            indexs(i,:)=J(1:K);
            Ii = J(1:K);             %i個近鄰點
            NI{i} = Ii;                 %放進單元數組N中
        end;
    else                                %   m<=N
        for i=1:N
            % Determine ki nearest neighbors of x_j
            x = data(:,i);      %data中的第i列所有元素存到存到x中
            dist2 = sum((data-repmat(x,[1 N])).^2,1);    
            [dist_sort,J] = sort(dist2);  
            indexs(i,:)=J(1:K);
            Ii = J(1:K);  
            NI{i} = Ii;
        end;
    end;

C_K=eye(K)-ones(K,K)/K;
current_Y=zeros(d,N);
next_Y=zeros(d,N);
L=zeros(N,N);
S_i=zeros(N,K);
L_i=zeros(N,K);
Y0=zeros(d,N); 
Y1=zeros(d,N);

for i=1:N 
    % Compute the d largest right singular eigenvectors of the cenantered matrix
    Ii =NI{i}; 
    %   data(:,Ii)表示data矩陣的所有行都是由矩陣Ii上的元素組成
    %   （這是矩陣的尋訪）
    Xi = data(:,Ii)-repmat(mean(data(:,Ii),2),[1,K]);  %中心化樣本矩陣，得到hatX
    %W = Xi'*Xi;                     %協方差矩陣（PSD）
    [Ui,Si,Vi] = svd(Xi);
    Om = Si * Vi' ;        
    Vi_=Vi(:,1:d);
    Y0(:,indexs(i,:)) =Om(1:d,:) ; %選取1：d行的全部元素，這裡得到Y0
    temp=C_K*(eye(K)-Vi_*Vi_');
    L_i(:,:)=0;
    L_i(indexs(i,:),:)=temp;     
    L=L+L_i*L_i';
end

L1=(L+L')/2;
L2=sparse(L1);
options.disp = 0;                   %
options.isreal = 1;                 %B為實矩陣
options.issym = 1;                  %矩陣B對稱
[V_vectors,e_values] = eigs(L2,d+1,0,options);          %D為d個最大特征值對角陣，U的列向量為對應的特征向量
%     lambda = diag(e_values);                       %将特征值取出來構成一個列向量
%     [lambda_s,J] = sort(abs(lambda));        %将特征值構成的列向量按升序排列
%     U = U(:,J);                               %将J放入矩陣U(U wei 全局對齊矩陣)
%     lambda = lambda(J);                             
%     T1 = U(:,2:d+1)';    
e_value=diag(e_values);
[e_value,sort_i]=sort(abs(e_value));
Y1=V_vectors(:,sort_i(2:d+1))';

current_Y=Y1;
for j=2:max_iterate_time
    L(:,:)=0;      
    for i=1:N
        Y_i=current_Y(:,indexs(i,:));
        Y_i=Y_i-repmat(mean(Y_i,2),1,K);
        
        [Ui,Si,Vi]=svd(Y_i);
%         Ti=Si;
%         Ti(find(Si~=0)) = 1./Si(find(Si~=0));
%         temp=C_K*(eye(K)-Vi*Ti'*Si*Vi');
        Vi_=Vi(:,1:d);
        temp=C_K*(eye(K)-Vi_*Vi_');
               
        L_i(:,:)=0;
        L_i(indexs(i,:),:)=temp;     
        L=L+L_i*L_i';
    end
   
    L1=(L+L')/2;
    L2=sparse(L1);
    options.disp = 0;                   %
    options.isreal = 1;                 %B為實矩陣
    options.issym = 1;                  %矩陣B對稱
    [V_vectors,e_values] = eigs(L2,d+1,0,options);          %D為d個最大特征值對角陣，U的列向量為對應的特征向量
%     lambda = diag(e_values);                       %将特征值取出來構成一個列向量
%     [lambda_s,J] = sort(abs(lambda));        %将特征值構成的列向量按升序排列
%     U = U(:,J);                               %将J放入矩陣U(U wei 全局對齊矩陣)
%     lambda = lambda(J);                             
%     T1 = U(:,2:d+1)';    
    e_value=diag(e_values);
    [e_value ,sort_i]=sort(abs(e_value));
    next_Y=V_vectors(:,sort_i(2:d+1))';
    error = trace(next_Y*L*next_Y');    
    current_Y=next_Y;
    fprintf('times:\t%d\t\terror:\t%d\n',j,abs(error));   
   if error<errors
        break;
    end  
end
Y=current_Y;