function s=SL0(A, x, sigma_min, sigma_decrease_factor, mu_0, L, A_pinv, true_s)
%
% SL0(A, x, sigma_min, sigma_decrease_factor, mu_0, L, A_pinv, true_s)
%
% Returns the sparsest vector s which satisfies underdetermined system of
% linear equations A*s=x, using Smoothed L0 (SL0) algorithm. Note that
% the matrix A should be a 'wide' matrix (more columns than rows). The
% number of the rows of matrix A should be equal to the length of the
% column vector x.
%
% The first 3 arguments should necessarily be provided by the user. The
% other parameters have defult values calculated within the function, or
% may be provided by the user.
%
% Sequence of Sigma (sigma_min and sigma_decrease_factor):
% This is a decreasing geometric sequence of positive numbers:
% - The first element of the sequence of sigma is calculated
% automatically. The last element is given by 'sigma_min', and the
% change factor for decreasing sigma is given by 'sigma_decrease_factor'.
% - The default value of 'sigma_decrease_factor' is 0.5. Larger value
% gives better results for less sparse sources, but it uses more steps
% on sigma to reach sigma_min, and hence it requires higher
% computational cost.
% - There is no default value for 'sigma_min', and it should be
% provided by the user (depending on his/her estimated source noise
% level, or his/her desired accuracy). By `noise' we mean here the
% noise in the sources, that is, the energy of the inactive elements of
% 's'. For example, by the noiseless case, we mean the inactive
% elements of 's' are exactly equal to zero. As a rule of tumb, for the
% noisy case, sigma_min should be about 2 to 4 times of the standard
% deviation of this noise. For the noiseless case, smaller 'sigma_min'
% results in better estimation of the sparsest solution, and hence its
% value is determined by the desired accuracy.
%
% mu_0:
% The value of mu_0 scales the sequence of mu. For each vlue of
% sigma, the value of mu is chosen via mu=mu_0*sigma^2. Note that this
% value effects Convergence.
% The default value is mu_0=2 (see the paper).
%
% L:
% number of iterations of the internal (steepest ascent) loop. The
% default value is L=3.
%
% A_pinv:
% is the pseudo-inverse of matrix A defined by A_pinv=A'*inv(A*A').
% If it is not provided, it will be calculated within the function. If
% you use this function for solving x(t)=A s(t) for different values of
% 't', it would be a good idea to calculate A_pinv outside the function
% to prevent its re-calculation for each 't'.
%
% true_s:
% is the true value of the sparse solution. This argument is for
% simulation purposes. If it is provided by the user, then the function
% will calculate the SNR of the estimation for each value of sigma and
% it provides a progress report.
%
% Authors: Massoud Babaie-Zadeh and Hossein Mohimani
% Version: 1.3
% Last modified: 4 August 2008.
%
%
% Web-page:
% ------------------
% http://ee.sharif.ir/~SLzero
%
% Code History:
%--------------
% Version 1.2: Adding some more comments in the help section
%
% Version 1.1: 4 August 2008
% - Using MATLAB's pseudo inverse function to generalize for the case
% the matrix A is not full-rank.
%
% Version 1.0 (first official version): 4 July 2008.
%
% First non-official version and algorithm development: Summer 2006
if nargin < 4
sigma_decrease_factor = 0.5;
A_pinv = pinv(A);
mu_0 = 2;
L = 3;
ShowProgress = logical(0);
elseif nargin == 4
A_pinv = pinv(A);
mu_0 = 2;
L = 3;
ShowProgress = logical(0);
elseif nargin == 5
A_pinv = pinv(A);
L = 3;
ShowProgress = logical(0);
elseif nargin == 6
A_pinv = pinv(A);
ShowProgress = logical(0);
elseif nargin == 7
ShowProgress = logical(0);
elseif nargin == 8
ShowProgress = logical(1);
else
error('Error in calling SL0 function');
end
% Initialization
%s = A\x;
%size(x)
size(A_pinv)
s = A_pinv*x;
sigma = 2*max(abs(s));
% Main Loop
while sigma>sigma_min
for i=1:L
delta = OurDelta(s,sigma);
s = s - mu_0*delta;
s = s - A_pinv*(A*s-x); % Projection
end
if ShowProgress
fprintf(' sigma=%f, SNR=%f\n',sigma,estimate_SNR(s,true_s))
end
sigma = sigma * sigma_decrease_factor;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function delta=OurDelta(s,sigma)
delta = s.*exp(-(s.*conj(s))/sigma^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function SNR=estimate_SNR(estim_s,true_s)
err = true_s - estim_s;
SNR = 10*log10(sum(true_s.*conj(true_s))/sum(err.*conj(err)));