Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/netlab3.3/conjgrad.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [x, options, flog, pointlog] = conjgrad(f, x, options, gradf, ...
wolffd@0 2 varargin)
wolffd@0 3 %CONJGRAD Conjugate gradients optimization.
wolffd@0 4 %
wolffd@0 5 % Description
wolffd@0 6 % [X, OPTIONS, FLOG, POINTLOG] = CONJGRAD(F, X, OPTIONS, GRADF) uses a
wolffd@0 7 % conjugate gradients algorithm to find the minimum of the function
wolffd@0 8 % F(X) whose gradient is given by GRADF(X). Here X is a row vector and
wolffd@0 9 % F returns a scalar value. The point at which F has a local minimum
wolffd@0 10 % is returned as X. The function value at that point is returned in
wolffd@0 11 % OPTIONS(8). A log of the function values after each cycle is
wolffd@0 12 % (optionally) returned in FLOG, and a log of the points visited is
wolffd@0 13 % (optionally) returned in POINTLOG.
wolffd@0 14 %
wolffd@0 15 % CONJGRAD(F, X, OPTIONS, GRADF, P1, P2, ...) allows additional
wolffd@0 16 % arguments to be passed to F() and GRADF().
wolffd@0 17 %
wolffd@0 18 % The optional parameters have the following interpretations.
wolffd@0 19 %
wolffd@0 20 % OPTIONS(1) is set to 1 to display error values; also logs error
wolffd@0 21 % values in the return argument ERRLOG, and the points visited in the
wolffd@0 22 % return argument POINTSLOG. If OPTIONS(1) is set to 0, then only
wolffd@0 23 % warning messages are displayed. If OPTIONS(1) is -1, then nothing is
wolffd@0 24 % displayed.
wolffd@0 25 %
wolffd@0 26 % OPTIONS(2) is a measure of the absolute precision required for the
wolffd@0 27 % value of X at the solution. If the absolute difference between the
wolffd@0 28 % values of X between two successive steps is less than OPTIONS(2),
wolffd@0 29 % then this condition is satisfied.
wolffd@0 30 %
wolffd@0 31 % OPTIONS(3) is a measure of the precision required of the objective
wolffd@0 32 % function at the solution. If the absolute difference between the
wolffd@0 33 % objective function values between two successive steps is less than
wolffd@0 34 % OPTIONS(3), then this condition is satisfied. Both this and the
wolffd@0 35 % previous condition must be satisfied for termination.
wolffd@0 36 %
wolffd@0 37 % OPTIONS(9) is set to 1 to check the user defined gradient function.
wolffd@0 38 %
wolffd@0 39 % OPTIONS(10) returns the total number of function evaluations
wolffd@0 40 % (including those in any line searches).
wolffd@0 41 %
wolffd@0 42 % OPTIONS(11) returns the total number of gradient evaluations.
wolffd@0 43 %
wolffd@0 44 % OPTIONS(14) is the maximum number of iterations; default 100.
wolffd@0 45 %
wolffd@0 46 % OPTIONS(15) is the precision in parameter space of the line search;
wolffd@0 47 % default 1E-4.
wolffd@0 48 %
wolffd@0 49 % See also
wolffd@0 50 % GRADDESC, LINEMIN, MINBRACK, QUASINEW, SCG
wolffd@0 51 %
wolffd@0 52
wolffd@0 53 % Copyright (c) Ian T Nabney (1996-2001)
wolffd@0 54
wolffd@0 55 % Set up the options.
wolffd@0 56 if length(options) < 18
wolffd@0 57 error('Options vector too short')
wolffd@0 58 end
wolffd@0 59
wolffd@0 60 if(options(14))
wolffd@0 61 niters = options(14);
wolffd@0 62 else
wolffd@0 63 niters = 100;
wolffd@0 64 end
wolffd@0 65
wolffd@0 66 % Set up options for line search
wolffd@0 67 line_options = foptions;
wolffd@0 68 % Need a precise line search for success
wolffd@0 69 if options(15) > 0
wolffd@0 70 line_options(2) = options(15);
wolffd@0 71 else
wolffd@0 72 line_options(2) = 1e-4;
wolffd@0 73 end
wolffd@0 74
wolffd@0 75 display = options(1);
wolffd@0 76
wolffd@0 77 % Next two lines allow conjgrad to work with expression strings
wolffd@0 78 f = fcnchk(f, length(varargin));
wolffd@0 79 gradf = fcnchk(gradf, length(varargin));
wolffd@0 80
wolffd@0 81 % Check gradients
wolffd@0 82 if (options(9))
wolffd@0 83 feval('gradchek', x, f, gradf, varargin{:});
wolffd@0 84 end
wolffd@0 85
wolffd@0 86 options(10) = 0;
wolffd@0 87 options(11) = 0;
wolffd@0 88 nparams = length(x);
wolffd@0 89 fnew = feval(f, x, varargin{:});
wolffd@0 90 options(10) = options(10) + 1;
wolffd@0 91 gradnew = feval(gradf, x, varargin{:});
wolffd@0 92 options(11) = options(11) + 1;
wolffd@0 93 d = -gradnew; % Initial search direction
wolffd@0 94 br_min = 0;
wolffd@0 95 br_max = 1.0; % Initial value for maximum distance to search along
wolffd@0 96 tol = sqrt(eps);
wolffd@0 97
wolffd@0 98 j = 1;
wolffd@0 99 if nargout >= 3
wolffd@0 100 flog(j, :) = fnew;
wolffd@0 101 if nargout == 4
wolffd@0 102 pointlog(j, :) = x;
wolffd@0 103 end
wolffd@0 104 end
wolffd@0 105
wolffd@0 106 while (j <= niters)
wolffd@0 107
wolffd@0 108 xold = x;
wolffd@0 109 fold = fnew;
wolffd@0 110 gradold = gradnew;
wolffd@0 111
wolffd@0 112 gg = gradold*gradold';
wolffd@0 113 if (gg == 0.0)
wolffd@0 114 % If the gradient is zero then we are done.
wolffd@0 115 options(8) = fnew;
wolffd@0 116 return;
wolffd@0 117 end
wolffd@0 118
wolffd@0 119 % This shouldn't occur, but rest of code depends on d being downhill
wolffd@0 120 if (gradnew*d' > 0)
wolffd@0 121 d = -d;
wolffd@0 122 if options(1) >= 0
wolffd@0 123 warning('search direction uphill in conjgrad');
wolffd@0 124 end
wolffd@0 125 end
wolffd@0 126
wolffd@0 127 line_sd = d./norm(d);
wolffd@0 128 [lmin, line_options] = feval('linemin', f, xold, line_sd, fold, ...
wolffd@0 129 line_options, varargin{:});
wolffd@0 130 options(10) = options(10) + line_options(10);
wolffd@0 131 options(11) = options(11) + line_options(11);
wolffd@0 132 % Set x and fnew to be the actual search point we have found
wolffd@0 133 x = xold + lmin * line_sd;
wolffd@0 134 fnew = line_options(8);
wolffd@0 135
wolffd@0 136 % Check for termination
wolffd@0 137 if (max(abs(x - xold)) < options(2) & max(abs(fnew - fold)) < options(3))
wolffd@0 138 options(8) = fnew;
wolffd@0 139 return;
wolffd@0 140 end
wolffd@0 141
wolffd@0 142 gradnew = feval(gradf, x, varargin{:});
wolffd@0 143 options(11) = options(11) + 1;
wolffd@0 144
wolffd@0 145 % Use Polak-Ribiere formula to update search direction
wolffd@0 146 gamma = ((gradnew - gradold)*(gradnew)')/gg;
wolffd@0 147 d = (d .* gamma) - gradnew;
wolffd@0 148
wolffd@0 149 if (display > 0)
wolffd@0 150 fprintf(1, 'Cycle %4d Function %11.6f\n', j, line_options(8));
wolffd@0 151 end
wolffd@0 152
wolffd@0 153 j = j + 1;
wolffd@0 154 if nargout >= 3
wolffd@0 155 flog(j, :) = fnew;
wolffd@0 156 if nargout == 4
wolffd@0 157 pointlog(j, :) = x;
wolffd@0 158 end
wolffd@0 159 end
wolffd@0 160 end
wolffd@0 161
wolffd@0 162 % If we get here, then we haven't terminated in the given number of
wolffd@0 163 % iterations.
wolffd@0 164
wolffd@0 165 options(8) = fold;
wolffd@0 166 if (options(1) >= 0)
wolffd@0 167 disp(maxitmess);
wolffd@0 168 end