ivan@155: function [unhat,er] = mm1(Phi,x,u0,to,lambda,maxIT,eps,map)
ivan@155: %% Iterative Soft Thresholding (with optional debiasing)
ivan@155: %
ivan@155: % Phi = Normalized Dictionary
ivan@155: % x = Signal(x). This can be a vector or a matrix
ivan@155: % u0 = Initial guess for the coefficients
ivan@155: % to = 1/(step size) . It is larger than spectral norm of dictionary Phi
ivan@155: % lambda = Lagrangian multiplier. (regulates shrinkage)
ivan@155: % eps = Stopping criterion for iterative softthresholding and MM dictionary update
ivan@155: % map = Debiasing. 0 = No, 1 = Yes
ivan@155: % unhat = Updated coefficients
ivan@155: % er = Objective cost
ivan@155: %%
ivan@155: cont = 1;
ivan@155: in = 1;
ivan@155: % un = zeros(size(u0,1),size(u0,2));    
ivan@155: un = u0;
ivan@155: c1 = (1/to^2)*Phi'*x;
ivan@155: c2 = (1/to^2)*(Phi'*Phi);
ivan@155: %%%%
ivan@155: while (cont && (in<=maxIT))
ivan@155:     unold = un;
ivan@155:     %%%%%% Soft Thresholding %%%%%%%
ivan@155:     alphap = (un + c1 - c2*un);
ivan@155:     un = (alphap-(lambda/(2*to^2))*sign(alphap)).*(abs(alphap)>=(lambda/(2*to^2)));
ivan@155:     in = in+1;
ivan@155:     cont = sum(sum((unold-un).^2))>eps;
ivan@155: end
ivan@155: %%%%%%%%%%
ivan@155: if map == 1,
ivan@155:     %% Mapping on the selected space %%%%
ivan@155:     [uN,uM] = size(un);
ivan@155:     unhat = zeros(uN,uM);
ivan@155:     for l = 1:uM,
ivan@155:         unz = (abs(un(:,l))>0);
ivan@155:         M = diag(unz);
ivan@155:         PhiNew = Phi*M;
ivan@155:         PhiS = PhiNew(:,unz);
ivan@155:         unt = inv(PhiS'*PhiS+.0001*eye(sum(unz)))*PhiS'*x(:,l);
ivan@155:         unhat(unz,l) = unt;
ivan@155:     end    
ivan@155: else
ivan@155:     unhat = un;
ivan@155: end
ivan@155: %%% Cost function calculation
ivan@155: if map == 1,
ivan@155:     er = sum(sum((Phi*unhat-x).^2))+lambda*(sum(sum(abs(unhat)>0))); %% l_0 Cost function     
ivan@155: else   
ivan@155:     er = sum(sum((Phi*unhat-x).^2))+lambda*(sum(sum(abs(unhat))));
ivan@155: end