function [X , S , r]= MultiLS(A,b,band,sparsity,supp,seqmax,ep)
%% locally optimize the object in a random order
%% input: A:        sensing matrix
%%        b:        measurement
%%        band:     radius of searching band in Local Optimization
%%        sparsity: # elements in supprec
%%        supp:     input support to be locally optimized
%%        seqmax:   set seqmax = 1
%%        ep:       noise level Ax + e = b where ||e||_2 \le ep. By
%%                  default, ep = 1e-12.
%% output: X:       output vector
%%         S:       output support which has been locally optimized
%%         r:       A*X-b

if(nargin<6)
    seqmax = 1;
    ep     = 1e-12;
end

%% apply OLS to local search 
[m,N] = size(A);

Cell = cell(sparsity,1);
for j = 1:sparsity
      Cell{j,1} = max(1,supp(j)-band):min(supp(j)+band,N);
end

%% first round
for seq = 1 : seqmax
   S = supp;
   order = randperm(sparsity);
   for iter = 1 : sparsity
       now = order(iter);
       Size = length(Cell{now,1});
       Residue = zeros(Size,1);
       for k = 1 : Size
           temsupport = setdiff(S,S(now));
           temsupport = union(temsupport,Cell{now,1}(k));
           Matrix = A(:,temsupport);
           x = lscov(Matrix,b);
           Residue(k) = norm(Matrix*x-b)/norm(b);
       end
       [minres , newindex] = min(Residue);
       newindex = Cell{now,1}(newindex);
       S(now) = newindex;
   end
   % update supp
   supp = S;
   %error
   A_S = A(:,S);
   X = zeros(N,1);
   x = lscov(A_S,b);
   X(S) = x;
   r = A_S*x-b;
   res = norm(r);
   % stop if residue is small and supp doesn't change
   %fprintf('seq = %6.0f , relative residue = %e \n',seq,relativeres)
   if (res < ep) %%|| isempty(setdiff(S,supp))
       break;
   end
end
        
