function [X , S , r]= LS(A,b,band,sparsity,supp)
%% locally optimize the object one by one according to the order of the input
%% support(supp)
%% 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
%% output: X: output vector
%% S: output support which has been locally optimized
%% r: A*X-b
[m,N] = size(A);

%% Neighborhood
Cell = cell(sparsity,1);
for j = 1:sparsity
BERegion = supp(j)-band : supp(j)+band;
if BERegion(1)<0.5
Ind = BERegion<0.5;
BERegion = BERegion.*(1-Ind)+(BERegion+N).*Ind;
elseif BERegion(end)>(N+0.5)
Ind = BERegion>(N+0.5);
BERegion = BERegion.*(1-Ind)+(BERegion-N).*Ind;
end
Cell{j,1} = BERegion;
end

%% Locally optimize in the neighborhood of each object
for iter = 1 : sparsity
Size = length(Cell{iter,1});
Residue = zeros(Size,1);
for k = 1 : Size
temsupport = setdiff(supp,supp(iter));
temsupport = union(temsupport,Cell{iter,1}(k));
Matrix = A(:,temsupport);
x = lscov(Matrix,b);
Residue(k) = norm(Matrix*x-b)/norm(b);
end
[minres , newindex] = min(Residue);
newindex = Cell{iter,1}(newindex);
supp(iter) = newindex;
end

S = supp;
A_S = A(:,S);
X = zeros(N,1);
x = lscov(A_S,b);
X(S)= x;
r = A_S*x-b;