function lyp_mean = lyapunov(X, num_neighbors, a, TS, num_steps)
    D_embed = size(X, 2);
    num_points = size(X,1);
    n=0;
    Q=eye(D_embed);
    lyp_est=[];
    d = [];
    jump = a;
    start = 1;

    for i=start:jump:num_points-a
        for k=1:num_points-a
            if (k~=i)
                d(k)=norm(X(i,:)-X(k,:));
            end;
        end;
        d(i)=max(d)+100;
%                    d(find(d<4.0)) = max(d)+100;
        [Y,I]=sort(d);
        I=I(1:num_neighbors);
        %d(I(num_neighbors));
        for k=1:num_neighbors
            B_i(k,:)=X(I(k),:)-X(i,:);
            B_i_plus_a(k,:)=X(I(k)+a,:)-X(i+a,:);
        end;
        T=(pinv(B_i)*B_i_plus_a)';

        [Q,R]=qr(T*Q);

        if (sum((abs(diag(R))>0))==D_embed)
            lyp_est=[lyp_est; log(abs(diag(R)))'/(TS*a)];
            n=n+1;
        else
            R
            disp('In lyapunov.m: trouble');
        end;

        if n == num_steps
             break;
        end;
    end;
    
    lyp_mean = mean(lyp_est);
        