Daniel@0: function [x, options, flog, pointlog] = conjgrad(f, x, options, gradf, ... Daniel@0: varargin) Daniel@0: %CONJGRAD Conjugate gradients optimization. Daniel@0: % Daniel@0: % Description Daniel@0: % [X, OPTIONS, FLOG, POINTLOG] = CONJGRAD(F, X, OPTIONS, GRADF) uses a Daniel@0: % conjugate gradients algorithm to find the minimum of the function Daniel@0: % F(X) whose gradient is given by GRADF(X). Here X is a row vector and Daniel@0: % F returns a scalar value. The point at which F has a local minimum Daniel@0: % is returned as X. The function value at that point is returned in Daniel@0: % OPTIONS(8). A log of the function values after each cycle is Daniel@0: % (optionally) returned in FLOG, and a log of the points visited is Daniel@0: % (optionally) returned in POINTLOG. Daniel@0: % Daniel@0: % CONJGRAD(F, X, OPTIONS, GRADF, P1, P2, ...) allows additional Daniel@0: % arguments to be passed to F() and GRADF(). Daniel@0: % Daniel@0: % The optional parameters have the following interpretations. Daniel@0: % Daniel@0: % OPTIONS(1) is set to 1 to display error values; also logs error Daniel@0: % values in the return argument ERRLOG, and the points visited in the Daniel@0: % return argument POINTSLOG. If OPTIONS(1) is set to 0, then only Daniel@0: % warning messages are displayed. If OPTIONS(1) is -1, then nothing is Daniel@0: % displayed. Daniel@0: % Daniel@0: % OPTIONS(2) is a measure of the absolute precision required for the Daniel@0: % value of X at the solution. If the absolute difference between the Daniel@0: % values of X between two successive steps is less than OPTIONS(2), Daniel@0: % then this condition is satisfied. Daniel@0: % Daniel@0: % OPTIONS(3) is a measure of the precision required of the objective Daniel@0: % function at the solution. If the absolute difference between the Daniel@0: % objective function values between two successive steps is less than Daniel@0: % OPTIONS(3), then this condition is satisfied. Both this and the Daniel@0: % previous condition must be satisfied for termination. Daniel@0: % Daniel@0: % OPTIONS(9) is set to 1 to check the user defined gradient function. Daniel@0: % Daniel@0: % OPTIONS(10) returns the total number of function evaluations Daniel@0: % (including those in any line searches). Daniel@0: % Daniel@0: % OPTIONS(11) returns the total number of gradient evaluations. Daniel@0: % Daniel@0: % OPTIONS(14) is the maximum number of iterations; default 100. Daniel@0: % Daniel@0: % OPTIONS(15) is the precision in parameter space of the line search; Daniel@0: % default 1E-4. Daniel@0: % Daniel@0: % See also Daniel@0: % GRADDESC, LINEMIN, MINBRACK, QUASINEW, SCG Daniel@0: % Daniel@0: Daniel@0: % Copyright (c) Ian T Nabney (1996-2001) Daniel@0: Daniel@0: % Set up the options. Daniel@0: if length(options) < 18 Daniel@0: error('Options vector too short') Daniel@0: end Daniel@0: Daniel@0: if(options(14)) Daniel@0: niters = options(14); Daniel@0: else Daniel@0: niters = 100; Daniel@0: end Daniel@0: Daniel@0: % Set up options for line search Daniel@0: line_options = foptions; Daniel@0: % Need a precise line search for success Daniel@0: if options(15) > 0 Daniel@0: line_options(2) = options(15); Daniel@0: else Daniel@0: line_options(2) = 1e-4; Daniel@0: end Daniel@0: Daniel@0: display = options(1); Daniel@0: Daniel@0: % Next two lines allow conjgrad to work with expression strings Daniel@0: f = fcnchk(f, length(varargin)); Daniel@0: gradf = fcnchk(gradf, length(varargin)); Daniel@0: Daniel@0: % Check gradients Daniel@0: if (options(9)) Daniel@0: feval('gradchek', x, f, gradf, varargin{:}); Daniel@0: end Daniel@0: Daniel@0: options(10) = 0; Daniel@0: options(11) = 0; Daniel@0: nparams = length(x); Daniel@0: fnew = feval(f, x, varargin{:}); Daniel@0: options(10) = options(10) + 1; Daniel@0: gradnew = feval(gradf, x, varargin{:}); Daniel@0: options(11) = options(11) + 1; Daniel@0: d = -gradnew; % Initial search direction Daniel@0: br_min = 0; Daniel@0: br_max = 1.0; % Initial value for maximum distance to search along Daniel@0: tol = sqrt(eps); Daniel@0: Daniel@0: j = 1; Daniel@0: if nargout >= 3 Daniel@0: flog(j, :) = fnew; Daniel@0: if nargout == 4 Daniel@0: pointlog(j, :) = x; Daniel@0: end Daniel@0: end Daniel@0: Daniel@0: while (j <= niters) Daniel@0: Daniel@0: xold = x; Daniel@0: fold = fnew; Daniel@0: gradold = gradnew; Daniel@0: Daniel@0: gg = gradold*gradold'; Daniel@0: if (gg == 0.0) Daniel@0: % If the gradient is zero then we are done. Daniel@0: options(8) = fnew; Daniel@0: return; Daniel@0: end Daniel@0: Daniel@0: % This shouldn't occur, but rest of code depends on d being downhill Daniel@0: if (gradnew*d' > 0) Daniel@0: d = -d; Daniel@0: if options(1) >= 0 Daniel@0: warning('search direction uphill in conjgrad'); Daniel@0: end Daniel@0: end Daniel@0: Daniel@0: line_sd = d./norm(d); Daniel@0: [lmin, line_options] = feval('linemin', f, xold, line_sd, fold, ... Daniel@0: line_options, varargin{:}); Daniel@0: options(10) = options(10) + line_options(10); Daniel@0: options(11) = options(11) + line_options(11); Daniel@0: % Set x and fnew to be the actual search point we have found Daniel@0: x = xold + lmin * line_sd; Daniel@0: fnew = line_options(8); Daniel@0: Daniel@0: % Check for termination Daniel@0: if (max(abs(x - xold)) < options(2) & max(abs(fnew - fold)) < options(3)) Daniel@0: options(8) = fnew; Daniel@0: return; Daniel@0: end Daniel@0: Daniel@0: gradnew = feval(gradf, x, varargin{:}); Daniel@0: options(11) = options(11) + 1; Daniel@0: Daniel@0: % Use Polak-Ribiere formula to update search direction Daniel@0: gamma = ((gradnew - gradold)*(gradnew)')/gg; Daniel@0: d = (d .* gamma) - gradnew; Daniel@0: Daniel@0: if (display > 0) Daniel@0: fprintf(1, 'Cycle %4d Function %11.6f\n', j, line_options(8)); Daniel@0: end Daniel@0: Daniel@0: j = j + 1; Daniel@0: if nargout >= 3 Daniel@0: flog(j, :) = fnew; Daniel@0: if nargout == 4 Daniel@0: pointlog(j, :) = x; Daniel@0: end Daniel@0: end Daniel@0: end Daniel@0: Daniel@0: % If we get here, then we haven't terminated in the given number of Daniel@0: % iterations. Daniel@0: Daniel@0: options(8) = fold; Daniel@0: if (options(1) >= 0) Daniel@0: disp(maxitmess); Daniel@0: end