Ax=A'*x;
for i=1:m,
Ax(:,i)=Ax(:,i)*mu(i);
end
tresh=max(Ax.^2-repmat(RHO,n,1),0);
f(iter)=sum(sum(tresh)); %cost function
if f(iter)==0, %the sparsity parameter is too high: all entries of the loading vector are zero.
break
else
grad=zeros(p,m);
for i=1:m,
pattern=find(tresh(:,i)'>0);
grad(:,i)=A(:,pattern)*Ax(pattern,i)*mu(i); %gradient
end
[U,S,V]=svd(grad,'econ');
x=U*V';
end