nikcleju@46: function [x,k,l] = fastProjection( U, S, V, y, b, epsilon, lambda0, DISP )
nikcleju@46: % [x,niter,lambda] = fastProjection(U, S, V, y, b, epsilon, [lambda0], [DISP] )
nikcleju@46: %
nikcleju@46: % minimizes || x - y ||
nikcleju@46: %   such that || Ax - b || <= epsilon
nikcleju@46: %
nikcleju@46: % where USV' = A (i.e the SVD of A)
nikcleju@46: %
nikcleju@46: % The optional input "lambda0" is a guess for the Lagrange parameter
nikcleju@46: %
nikcleju@46: % Warning: for speed, does not calculate A(y) to see if x = y is feasible
nikcleju@46: %
nikcleju@46: % NESTA Version 1.1
nikcleju@46: %   See also Core_Nesterov
nikcleju@46: 
nikcleju@46: % Written by Stephen Becker, September 2009, srbecker@caltech.edu
nikcleju@46: 
nikcleju@46: DEBUG = true;
nikcleju@46: if nargin < 8
nikcleju@46:     DISP = false;
nikcleju@46: end
nikcleju@46: % -- Parameters for Newton's method --
nikcleju@46: MAXIT = 70;
nikcleju@46: TOL = 1e-8 * epsilon;
nikcleju@46: % TOL = max( TOL, 10*eps );  % we can't do better than machine precision
nikcleju@46: 
nikcleju@46: m = size(U,1);
nikcleju@46: n = size(V,1);
nikcleju@46: mn = min([m,n]);
nikcleju@46: if numel(S) > mn^2, S = diag(diag(S)); end  % S should be a small square matrix
nikcleju@46: r = size(S);
nikcleju@46: if size(U,2) > r, U = U(:,1:r); end
nikcleju@46: if size(V,2) > r, V = V(:,1:r); end
nikcleju@46: 
nikcleju@46: s = diag(S);
nikcleju@46: s2 = s.^2;
nikcleju@46: 
nikcleju@46: % What we want to do:
nikcleju@46: %   b = b - A*y;
nikcleju@46: %   bb = U'*b;
nikcleju@46: 
nikcleju@46: % if A doesn't have full row rank, then b may not be in the range
nikcleju@46: if size(U,1) > size(U,2)
nikcleju@46:     bRange = U*(U'*b);
nikcleju@46:     bNull = b - bRange;
nikcleju@46:     epsilon = sqrt( epsilon^2 - norm(bNull)^2 );
nikcleju@46: end
nikcleju@46: b = U'*b - S*(V'*y);  % parenthesis is very important!  This is expensive.
nikcleju@46:     
nikcleju@46: % b2 = b.^2;
nikcleju@46: b2 = abs(b).^2;  % for complex data
nikcleju@46: bs2 = b2.*s2;
nikcleju@46: epsilon2 = epsilon^2;
nikcleju@46: 
nikcleju@46: % The following routine need to be fast
nikcleju@46: % For efficiency (at cost of transparency), we are writing the calculations
nikcleju@46: % in a way that minimize number of operations.  The functions "f"
nikcleju@46: % and "fp" represent f and its derivative.
nikcleju@46: 
nikcleju@46: % f = @(lambda) sum( b2 .*(1-lambda*s2).^(-2) ) - epsilon^2;
nikcleju@46: % fp = @(lambda) 2*sum( bs2 .*(1-lambda*s2).^(-3) );
nikcleju@46: if nargin < 7, lambda0 = 0; end
nikcleju@46: l = lambda0; oldff = 0;
nikcleju@46: one = ones(m,1);
nikcleju@46: alpha = 1;      % take full Newton steps
nikcleju@46: for k = 1:MAXIT
nikcleju@46:     % make f(l) and fp(l) as efficient as possible:
nikcleju@46:     ls = one./(one-l*s2);
nikcleju@46:     ls2 = ls.^2;
nikcleju@46:     ls3 = ls2.*ls;
nikcleju@46:     ff = b2.'*ls2; % should be .', not ', even for complex data
nikcleju@46:     ff = ff - epsilon2;
nikcleju@46:     fpl = 2*( bs2.'*ls3 );  % should be .', not ', even for complex data
nikcleju@46: %     ff = f(l);    % this is a little slower
nikcleju@46: %     fpl = fp(l);  % this is a little slower
nikcleju@46:     d = -ff/fpl;
nikcleju@46:     if DISP, fprintf('%2d, lambda is %5.2f, f(lambda) is %.2e, f''(lambda) is %.2e\n',...
nikcleju@46:             k,l,ff,fpl ); end
nikcleju@46:     if abs(ff) < TOL, break; end        % stopping criteria
nikcleju@46:     l_old = l;
nikcleju@46:     if k>2 && ( abs(ff) > 10*abs(oldff+100) ) %|| abs(d) > 1e13 )
nikcleju@46:         l = 0; alpha = 1/2;  
nikcleju@46: %         oldff = f(0);
nikcleju@46:         oldff = sum(b2); oldff = oldff - epsilon2;
nikcleju@46:         if DISP, disp('restarting'); end
nikcleju@46:     else
nikcleju@46:         if alpha < 1, alpha = (alpha+1)/2; end
nikcleju@46:         l = l + alpha*d;
nikcleju@46:         oldff = ff;
nikcleju@46:         if l > 0
nikcleju@46:             l = 0;  % shouldn't be positive
nikcleju@46:             oldff = sum(b2);  oldff = oldff - epsilon2;
nikcleju@46:         end
nikcleju@46:     end
nikcleju@46:     if l_old == l && l == 0
nikcleju@46:         if DISP, disp('Making no progress; x = y is probably feasible'); end
nikcleju@46:         break;
nikcleju@46:     end
nikcleju@46: end
nikcleju@46: % if k == MAXIT && DEBUG, disp('maxed out iterations'); end
nikcleju@46: if l < 0
nikcleju@46:     xhat = -l*s.*b./( 1 - l*s2 );
nikcleju@46:     x = V*xhat + y;
nikcleju@46: else
nikcleju@46:     % y is already feasible, so no need to project
nikcleju@46:     l = 0;
nikcleju@46:     x = y;
nikcleju@46: end