Description of bbnnls

0001 function out = bbnnls(A, b, x0, opt)
0002 % BBNNLS   -- Solve NNLS problems via SBB
0003 %
0004 % WARNING Use at own risk!
0005 % NOTE --- guaranteed convergence phase: *REMOVED* for speedup!!
0006 % NOTE --- To speed up code further, *REMOVE* debugging part
0007 %
0008 %
0009 % function out = bbnnls(A, b, x0, opt)
0010 % Solve a bound-constrained least squares problem,
0011 %    min    0.5*||Ax-b||^2, s.t. x >= 0
0012 %
0013 %
0014 % x0 -- Starting vector (useful for warm-starts).
0015 %
0016 % OPT -- This structure contains important opt that control how the
0017 % optimization procedure runs. To obtain a default structure the user can
0018 % use 'opt = solopt'. Use 'help solopt' to get a description of
0019 % what the individual opt mean.
0020 %
0021 % Most important options to tune as: opt.tolg, opt.maxit
0022 %
0023 %
0024 % OUT contains the solution and other information about the optimization or
0025 % info indicating whether the method succeeded or failed.
0026 %
0027 % See also: solopt, bcls
0028 %
0029 % Version 1.1 (c) 2010 Suvrit Sra, Dongmin Kim
0030 %
0031 % Released under the GNU General Public License
0032 % For details on license and terms please see http://www.gnu.org/copyleft/gpl.html
0033 
0034 
0035     fgx = @(x) funcGrad(A,b, x); % function to compute obj and grad
0036 
0037     % do some initialization for maintaining statistics
0038     out.iter = 0;
0039     out.iterTimes = nan*ones(opt.maxit,1);
0040     out.objTimes  = nan*ones(opt.maxit,1);
0041     out.pgTimes   = nan*ones(opt.maxit,1);
0042     out.trueError = nan*ones(opt.maxit,1);
0043     out.startTime = tic;
0044     out.status = 'Failure';
0045 
0046     % HINT: Very important for overall speed is to have a good x0
0047     out.x      = x0;
0048     out.refx   = x0;
0049     [out.refobj, out.grad]   = fgx(out.x);
0050     out.oldg   = out.grad;
0051     out.refg   = out.oldg;
0052 
0053 
0054     %% Begin the main algorithm
0055     if (opt.verbose)
0056        fprintf('Running: **** SBB-NNLS ****\n\n');
0057        fprintf('Iter   \t     Obj\t\t  ||pg||_inf\t\t ||x-x*||\n');
0058        fprintf('-------------------------------------------------------\n');
0059     end
0060 
0061     objectives = zeros(opt.maxit,1);
0062     %f = figure;
0063     while 1
0064         out.iter = out.iter + 1;
0065 
0066         % HINT: edit checkTermination to determine good criterion for yourself!
0067         [termReason, out.pgTimes(out.iter)] = checkTermination(opt, out);
0068         if (termReason > 0), break; end
0069 
0070         % HINT: computeBBStep is the one to implement most carefully
0071         [step out] = computeBBStep(A, b, out);
0072         out.x = out.x - step * out.grad;
0073         out.oldg = out.grad;
0074         
0075         % HINT: projection step: can replace the 0 by an epsilon to truncate
0076         % values close to 0
0077         out.x(out.x < 0) = 0;
0078 
0079         [out.obj out.grad] = fgx(out.x);
0080         
0081         objectives(out.iter) = out.obj;
0082 %         clf;
0083 %         plot(objectives);
0084 %         title('Objective ||Ax-b||^2');
0085 %         xlabel('Iteration');
0086 %         ylabel('Objective');
0087 %         drawnow;
0088         
0089         % HINT: can remove, as this is just for statistics
0090         out.objTimes (out.iter) = out.obj;
0091         out.iterTimes(out.iter) = toc(out.startTime);
0092         
0093         % HINT: for debugging, to see how result develops if true x* is known
0094         if (opt.truex), out.trueError(out.iter) = norm(opt.xt-out.x); end
0095         if (opt.verbose)
0096             fprintf('%04d\t %E\t%E\t%E\n', out.iter, out.obj, out.pgTimes(out.iter), out.trueError(out.iter)); 
0097         end
0098     end % of while
0099 
0100     %%  Final statistics and wrap up
0101     out.time = toc(out.startTime);
0102     out.status = 'Success';
0103     out.termReason = setTermReason(termReason);
0104 end
0105 
0106 % Compute BB step; for SBB also modifies out.oldg, and this change must be
0107 % passed back to the calling routine, else it will fail!
0108 function [step out] = computeBBStep(A, b, out)
0109     
0110     % HINT: Can tune the x==0 to replace it by an epsilon to do TUNING
0111     gp = find(out.x == 0 & out.grad > 0);
0112     out.oldg(gp) = 0;
0113 
0114     Ag = A*out.oldg; % A*oldg
0115     
0116     % HINT: In my experience, the falling alternating steps perform better
0117     if (mod(out.iter, 2) == 0)
0118         step = (out.oldg' * out.oldg) / (Ag' * Ag);
0119     else
0120         numer = Ag' * Ag;
0121         Ag = A'*Ag; %
0122         Ag(gp) = 0;
0123         step = numer / (Ag' * Ag);
0124     end
0125 end
0126 
0127 % compute obj function and gradient --- requires good implementation of A*x
0128 % and A'*y for appropriate x and y
0129 function [f g] = funcGrad(A, b, x)
0130     Ax = A*x - b;
0131     f = 0.5*norm(Ax)^2;
0132     if (nargout > 1)
0133         g = A'*Ax;
0134     end
0135 end
0136 
0137 % check various termination criteria; return norm of pg
0138 % the strictest is norm of pg
0139 % HINT: for speedup, use maybe just opt.tolo or some other criterion that
0140 % you like.
0141 function [v pg] = checkTermination(options, out)
0142     % pgnorm limit -- need to check this first of all
0143     gp = find( (out.x ~= 0 | out.grad < 0));
0144 
0145     pg = norm(out.grad(gp), 'inf');
0146     if (pg < options.tolg), v=8; return; end
0147 
0148     % First check if we are doing termination based on running time
0149     if (options.time_limit)
0150         out.time = etime(clock, out.start_time);
0151         if (out.time >= options.maxtime)
0152             v = 1;
0153             return;
0154         end
0155     end
0156 
0157     % Now check if we are doing break by tolx
0158     if (options.use_tolx)
0159         if (norm(out.x-out.oldx)/norm(out.oldx) < options.tolx)
0160             v = 2;
0161             return;
0162         end
0163     end
0164 
0165     % Are we doing break by tolo (tol obj val)
0166     if (options.use_tolo && out.iter > 2)
0167         delta = abs(out.objTimes(out.iter-1)-out.objTimes(out.iter-2));
0168         if (delta < options.tolo)
0169             v = 3;
0170             return;
0171         end
0172     end
0173 
0174     % Finally the plain old check if max iter has been achieved
0175     if (out.iter >= options.maxit)
0176         v = 4;
0177         return;
0178     end
0179 
0180     % KKT violation
0181     if (options.use_kkt)
0182         if abs(out.x' * out.grad) <= options.tolk
0183             v = 7;
0184             return;
0185         end
0186     end
0187 
0188 
0189     % All is ok...
0190     v = 0;
0191 end
0192 
0193 %% Prints status
0194 function showStatus(out, options)
0195     if (options.verbose)
0196         fprintf('.');
0197         if (mod(out.iter, 30) == 0)
0198             fprintf('\n');
0199         end
0200     end
0201 end
0202 
0203 % String representation of termination
0204 function r = setTermReason(t)
0205     switch t
0206       case 1
0207         r = 'Exceeded time limit';
0208       case 2
0209         r = 'Relative change in x small enough';
0210       case 3
0211         r = 'Relative change in objvalue small enough';
0212       case 4
0213         r = 'Maximum number of iterations reached';
0214       case 5
0215         r = '|x_t+1 - x_t|=0 or |grad_t+1 - grad_t| < 1e-9';
0216       case 6
0217         r = 'Line search failed';
0218       case 7
0219         r = '|x^T * grad| < opt.pbb_gradient_norm';
0220       case 8
0221         r = '|| grad ||_inf < opt.tolg';
0222       case 100
0223         r = 'The active set converged';
0224       otherwise
0225         r = 'Undefined';
0226     end
0227 end
0228
bbnnls

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
