Code covered by the BSD License  

Highlights from
Sparse Randomized Kaczmarz for Multiple Measurement Vectors.

Sparse Randomized Kaczmarz for Multiple Measurement Vectors.

by

 

Implementation of sparse randomized Kaczmarz algorithm to handle multiple measurements.

Xhat=srkMMV(A,B,J,estSupp,numsol)
function Xhat=srkMMV(A,B,J,estSupp,numsol)
%This function implements our proposed modification in sparse randomize
%algorithm(see reference in myreadme.txt) to handle applications with multiple measurement vectors. 
% Input :
% A is m by n matrix
% B is n by L observation matrix
% J is no. times you want to sweep through all the rows of A
% estSupp : it is the initial estimate of no. of non-zero rows
%          it should be around twice of actual value for best results
%numsol : it takes two values 'all' or 'last'. It just tells how many
%         take the value 'all' if you want estimated 'X' after every sweep
%         or set the value 'last' if you just want the last sweep result.
%          'last' is obviously the best result.
%Output: 
%Xhat : it is the estimated value of actual matrix X.

[m,n]=size(A);
[~,L]=size(B);
%J represents total no. of sweeps 
X=zeros (n,L); % it is the reconstructed signal
B1=A.^2;
p=sum(B1,2)/sum(B1(:)); %Probability of each row
idx=randsample(1:m,J*m,true,p);
Xfinal=zeros(n,L,J);
w=zeros(1,n);

tt=1;
for j=1:J*m %sweeps
       %select the rows based on l2-norm
       [~,id]=sort(sqrt(sum(X.^2,2)),'descend'); 
        maxid=max(estSupp,n-j+1);
        S=id(1:maxid);
        %set the weighting criterion
        w(S)=1;
        Sc=setdiff(1:n,S);
        w(Sc)=1/sqrt(j);         
        a=w.*A(idx(j),:);
        %Repeat the projection L times 
        Xtemp=zeros(n,L);
        for i=1:L
          b=B(idx(j),i);
          x=X(:,i);         
          Xtemp(:,i)= ((b-a*x)/(a*a'))*a';
        end
        X=X+Xtemp;
      %store the values after every sweep
    if mod(j,m)==0
        Xfinal(:,:,tt)=X;
        tt=tt+1;
    end   
end

if strcmp(numsol,'all')
    Xhat=Xfinal;
elseif strcmp(numsol,'last')
    Xhat=Xfinal(:,:,J);
end


%err

Contact us