% iteratively reweighted least-squares method
n = 1000;
m = round(n/2);
A = randn(m,n);
xe= zeros(n,1);
ind = rand(n,1);
ind = find(ind >.9);
xe(ind) = 10*randn(size(ind));
A  = A/norm(A);
be = A*xe;
b  = be + max(abs(be))*1e-5*randn(m,1);
al = 1e-2;
figure(10), plot(1:n,xe,'r', 1:m,b,'b'), legend('solution','data')


%%  ill-posed problem: first-order derivative
clear, close all
n = 100;
m = n;
N = 500;
A = ones(n);
A = tril(A);
xe= zeros(n,1);
ind = rand(n,1);
ind = find(ind >.4);
xe(ind) = 10*randn(size(ind));
A  = A/norm(A);
be = A*xe;
b  = be + max(abs(be))*1e-5*randn(m,1);
al = 1e-5;
figure(10), plot(1:n,xe,'r', 1:m,b,'b'), legend('solution','data')
%}

%%
% iteratively reweighed least-squares
x = zeros(n,1);
res = [];
fcnl= [];
for k = 1:200
    W = diag(2./(abs(x)+1e-10));
    x = (A'*A+al*W)\(A'*b);
    figure(1), plot(1:n,x,'r*',1:n,xe,'k.'), axis square, set(gca,'fontsize',20), title(['IRLS, iter = ' num2str(k)])
    res(k) = norm(A*x-b);
    fcnl(k)= norm(A*x-b)^2/2+al*norm(x,1);
    pause
end

%%
figure(2), plot(1:200,res,'k',1:200,fcnl,'r','linewidth',2), 
legend('residual','functional value')
axis square
set(gca,'fontsize',20)

% fast iterative soft thresholding
%% FIST 
x0 = zeros(n,1);
y  = x0;
t1 = 1;
fcnl_fista = [];
err_fista = [];
res_fista = [];
for k = 1:N
    % update x
    x1 = y - A'*(A*y-b);
    x1 = sign(x1).*max(abs(x1)-al,0);
    % update t
    t2 = (1+sqrt(1+4*t1^2))/2;
    % update y
    y = x1 + (t1-1)/t2*(x1-x0);
    x0 = x1;
    t1 = t2;
    figure(3), plot(1:n,x1,'r.',1:n,xe,'k.'), axis square, set(gca,'fontsize',20), title(['FISTA, iter ' num2str(k)])
    fcnl_fista(k) = norm(A*x1-b)^2/2+al*norm(x1,1);
    err_fista(k) = norm(x1-xe);
    res_fista(k) = norm(A*x-b);
end
x_fista = x;
figure(4), plot(1:N,res_fista,'k',1:N,fcnl_fista,'r',1:N,err_fista,'g','linewidth',2)
legend('residual','functional','error'), axis square, set(gca,'fontsize',20),
axis square