% Section 6.3.1, Figure 6.7
% Original by Lieven Vandenberghe
% Adapted for CVX Argyris Zymnis - 10/2005
%
% Solves
%        minimize   ||A*x-b||_2
%        subject to card(x) <= k
%
% where card(x) denotes the number of nonzero elements in x,
% by first solving (for some value of alpha close to ||x_ln||_1)
%        minimize   ||A*x-b||_2
%        subject to ||x||_1 <= alpha
%
% and iteratively decreasing alpha so as to get card(x) = k
% The sparsity pattern is then fixed in A and b and
%        minimize   ||A*x-b||_2
%
% is solved

rand('state',0);

m = 10;
n = 20;

A = randn(m,n);
b = A*[randn(round(m/2),1); zeros(n-round(m/2),1)];
b = b + 0.1*norm(b)*randn(m,1);

if (1) %%%%%%%%%%%%

%
% min.  ||Ax-b||_2
% s.t.  ||x||_1 <= alpha

residuals_heur = [norm(b)];
xln = A'*((A*A')\b);
lnorm = norm(xln,1);
nopts = 100;
alphas = linspace(0,lnorm,nopts);
residuals_heur = [norm(b)];
card_heur = [0];

for k=2:(nopts-1)
alpha = alphas(k);

cvx_begin quiet
variable x(n)
minimize(norm(A*x-b))
subject to
norm(x,1) <= alpha;
cvx_end

x(find(abs(x) < 1e-3*max(abs(x)))) = 0;
ind = find(abs(x));
sparsity = length(ind);
fprintf(1,'Current sparsity pattern k = %d \n',sparsity);
x = zeros(n,1);  x(ind) = A(:,ind)\b;
card_heur = [card_heur, sparsity];
residuals_heur = [residuals_heur, norm(A*x-b)];
end;

obj1 = norm(b)
obj2 = [0];

i=1;
for k=1:m-1
if ~isempty(find(card_heur == k))
obj2(i+1) = k;
obj1(i+1) = min(residuals_heur(find(card_heur ==k)));
i=i+1;
end;
end;
obj2(i) = m;  obj1(i) = 0;

end; %%%%%%%%%%%%%%%%%%%

if (1) %%%%%%%%%%%%%

bestx = zeros(n,m);
bestres = zeros(1,m);

for k=1:m-1
k
% enumerate sparsity patterns with exactly k nonzeros
bestres(k) = Inf;
ind = 1:k
nocases = 1;
done = 0;
while ~done
done = 1;
for i=0:k-1
if (ind(k-i) < n-i),
ind(k-i:k) = ind(k-i)+[1:i+1];
done = 0;
break;
end;
end;
if done, break; end;
x = zeros(n,1);
x(ind) = A(:,ind)\b;
if (norm(A*x-b) < bestres(k)),
bestres(k) = norm(A*x-b);
bestx(:,k) = x;
end;
nocases = nocases + 1;
end;
nocases
factorial(n)/(factorial(n-k)*factorial(k))
end;

x = A\b;
bestres(m) = norm(A*x-b);
bestres = [norm(b) bestres];

end; %%%%%%%%%

figure
hold off
obj1dbl =[];
obj2dbl =[];
for i=1:length(obj1)-1
obj1dbl = [obj1dbl, obj1(i), obj1(i)];
obj2dbl = [obj2dbl, obj2(i), obj2(i+1)];
end;
obj1dbl = [obj1dbl, obj1(length(obj1))];
obj2dbl = [obj2dbl, obj2(length(obj1))];

bestobj1 = bestres;
bestobj2 = [0:1:m];
bestobj1dbl =[];
bestobj2dbl =[];
for i=1:length(bestobj1)-1
bestobj1dbl = [bestobj1dbl, bestobj1(i), bestobj1(i)];
bestobj2dbl = [bestobj2dbl, bestobj2(i), bestobj2(i+1)];
end;
bestobj1dbl = [bestobj1dbl, bestobj1(length(bestobj1))];
bestobj2dbl = [bestobj2dbl, bestobj2(length(bestobj1))];

plot(obj1dbl,obj2dbl,'-', bestobj1dbl, bestobj2dbl,'--');
hold on
plot(obj1,obj2,'o', bestobj1, bestobj2,'o');
axis([0 ceil(2*norm(b))/2 0 m+1])
xlabel('x');
ylabel('y');
hold off

%print -deps sparse_regressor_global_helv.eps
%save regressor_results
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 1
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 2
Current sparsity pattern k = 3
Current sparsity pattern k = 3
Current sparsity pattern k = 4
Current sparsity pattern k = 4
Current sparsity pattern k = 4
Current sparsity pattern k = 4
Current sparsity pattern k = 4
Current sparsity pattern k = 4
Current sparsity pattern k = 4
Current sparsity pattern k = 4
Current sparsity pattern k = 5
Current sparsity pattern k = 5
Current sparsity pattern k = 5
Current sparsity pattern k = 5
Current sparsity pattern k = 5
Current sparsity pattern k = 5
Current sparsity pattern k = 5
Current sparsity pattern k = 6
Current sparsity pattern k = 7
Current sparsity pattern k = 8
Current sparsity pattern k = 7
Current sparsity pattern k = 7
Current sparsity pattern k = 7
Current sparsity pattern k = 9
Current sparsity pattern k = 10
Current sparsity pattern k = 10
Current sparsity pattern k = 10
Current sparsity pattern k = 10
Current sparsity pattern k = 10
Current sparsity pattern k = 10
Current sparsity pattern k = 10
Current sparsity pattern k = 9
Current sparsity pattern k = 9
Current sparsity pattern k = 8
Current sparsity pattern k = 8
Current sparsity pattern k = 8
Current sparsity pattern k = 8
Current sparsity pattern k = 8
Current sparsity pattern k = 9
Current sparsity pattern k = 10
Current sparsity pattern k = 12
Current sparsity pattern k = 14
Current sparsity pattern k = 15
Current sparsity pattern k = 16
Current sparsity pattern k = 19
Current sparsity pattern k = 18
Current sparsity pattern k = 18
Current sparsity pattern k = 18
Current sparsity pattern k = 19
Current sparsity pattern k = 19
Current sparsity pattern k = 19
Current sparsity pattern k = 19
Current sparsity pattern k = 18
Current sparsity pattern k = 18
Current sparsity pattern k = 19
Current sparsity pattern k = 19
Current sparsity pattern k = 20
Current sparsity pattern k = 19
Current sparsity pattern k = 20
Current sparsity pattern k = 20
Current sparsity pattern k = 20
Current sparsity pattern k = 19
Current sparsity pattern k = 20
Current sparsity pattern k = 19
Current sparsity pattern k = 19
Current sparsity pattern k = 20
Current sparsity pattern k = 20
Current sparsity pattern k = 20
Current sparsity pattern k = 20

obj1 =

2.8954

k =

1

ind =

1

nocases =

20

ans =

20

k =

2

ind =

1     2

nocases =

190

ans =

190

k =

3

ind =

1     2     3

nocases =

1140

ans =

1140

k =

4

ind =

1     2     3     4

nocases =

4845

ans =

4845

k =

5

ind =

1     2     3     4     5

nocases =

15504

ans =

15504

k =

6

ind =

1     2     3     4     5     6

nocases =

38760

ans =

38760

k =

7

ind =

1     2     3     4     5     6     7

nocases =

77520

ans =

77520

k =

8

ind =

1     2     3     4     5     6     7     8

nocases =

125970

ans =

125970

k =

9

ind =

1     2     3     4     5     6     7     8     9

nocases =

167960

ans =

167960