# ml-class / ex3 / fmincg.m

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175``` ```function [X, fX, i] = fmincg(f, X, options, P1, P2, P3, P4, P5) % Minimize a continuous differentialble multivariate function. Starting point % is given by "X" (D by 1), and the function named in the string "f", must % return a function value and a vector of partial derivatives. The Polack- % Ribiere flavour of conjugate gradients is used to compute search directions, % and a line search using quadratic and cubic polynomial approximations and the % Wolfe-Powell stopping criteria is used together with the slope ratio method % for guessing initial step sizes. Additionally a bunch of checks are made to % make sure that exploration is taking place and that extrapolation will not % be unboundedly large. The "length" gives the length of the run: if it is % positive, it gives the maximum number of line searches, if negative its % absolute gives the maximum allowed number of function evaluations. You can % (optionally) give "length" a second component, which will indicate the % reduction in function value to be expected in the first line-search (defaults % to 1.0). The function returns when either its length is up, or if no further % progress can be made (ie, we are at a minimum, or so close that due to % numerical problems, we cannot get any closer). If the function terminates % within a few iterations, it could be an indication that the function value % and derivatives are not consistent (ie, there may be a bug in the % implementation of your "f" function). The function returns the found % solution "X", a vector of function values "fX" indicating the progress made % and "i" the number of iterations (line searches or function evaluations, % depending on the sign of "length") used. % % Usage: [X, fX, i] = fmincg(f, X, options, P1, P2, P3, P4, P5) % % See also: checkgrad % % Copyright (C) 2001 and 2002 by Carl Edward Rasmussen. Date 2002-02-13 % % % (C) Copyright 1999, 2000 & 2001, Carl Edward Rasmussen % % Permission is granted for anyone to copy, use, or modify these % programs and accompanying documents for purposes of research or % education, provided this copyright notice is retained, and note is % made of any changes that have been made. % % These programs and documents are distributed without any warranty, % express or implied. As the programs were written for research % purposes only, they have not been tested to the degree that would be % advisable in any important application. All use of these programs is % entirely at the user's own risk. % % [ml-class] Changes Made: % 1) Function name and argument specifications % 2) Output display % % Read options if exist('options', 'var') && ~isempty(options) && isfield(options, 'MaxIter') length = options.MaxIter; else length = 100; end RHO = 0.01; % a bunch of constants for line searches SIG = 0.5; % RHO and SIG are the constants in the Wolfe-Powell conditions INT = 0.1; % don't reevaluate within 0.1 of the limit of the current bracket EXT = 3.0; % extrapolate maximum 3 times the current bracket MAX = 20; % max 20 function evaluations per line search RATIO = 100; % maximum allowed slope ratio argstr = ['feval(f, X']; % compose string used to call function for i = 1:(nargin - 3) argstr = [argstr, ',P', int2str(i)]; end argstr = [argstr, ')']; if max(size(length)) == 2, red=length(2); length=length(1); else red=1; end S=['Iteration ']; i = 0; % zero the run length counter ls_failed = 0; % no previous line search has failed fX = []; [f1 df1] = eval(argstr); % get function value and gradient i = i + (length<0); % count epochs?! s = -df1; % search direction is steepest d1 = -s'*s; % this is the slope z1 = red/(1-d1); % initial step is red/(|s|+1) while i < abs(length) % while not finished i = i + (length>0); % count iterations?! X0 = X; f0 = f1; df0 = df1; % make a copy of current values X = X + z1*s; % begin line search [f2 df2] = eval(argstr); i = i + (length<0); % count epochs?! d2 = df2'*s; f3 = f1; d3 = d1; z3 = -z1; % initialize point 3 equal to point 1 if length>0, M = MAX; else M = min(MAX, -length-i); end success = 0; limit = -1; % initialize quanteties while 1 while ((f2 > f1+z1*RHO*d1) | (d2 > -SIG*d1)) & (M > 0) limit = z1; % tighten the bracket if f2 > f1 z2 = z3 - (0.5*d3*z3*z3)/(d3*z3+f2-f3); % quadratic fit else A = 6*(f2-f3)/z3+3*(d2+d3); % cubic fit B = 3*(f3-f2)-z3*(d3+2*d2); z2 = (sqrt(B*B-A*d2*z3*z3)-B)/A; % numerical error possible - ok! end if isnan(z2) | isinf(z2) z2 = z3/2; % if we had a numerical problem then bisect end z2 = max(min(z2, INT*z3),(1-INT)*z3); % don't accept too close to limits z1 = z1 + z2; % update the step X = X + z2*s; [f2 df2] = eval(argstr); M = M - 1; i = i + (length<0); % count epochs?! d2 = df2'*s; z3 = z3-z2; % z3 is now relative to the location of z2 end if f2 > f1+z1*RHO*d1 | d2 > -SIG*d1 break; % this is a failure elseif d2 > SIG*d1 success = 1; break; % success elseif M == 0 break; % failure end A = 6*(f2-f3)/z3+3*(d2+d3); % make cubic extrapolation B = 3*(f3-f2)-z3*(d3+2*d2); z2 = -d2*z3*z3/(B+sqrt(B*B-A*d2*z3*z3)); % num. error possible - ok! if ~isreal(z2) | isnan(z2) | isinf(z2) | z2 < 0 % num prob or wrong sign? if limit < -0.5 % if we have no upper limit z2 = z1 * (EXT-1); % the extrapolate the maximum amount else z2 = (limit-z1)/2; % otherwise bisect end elseif (limit > -0.5) & (z2+z1 > limit) % extraplation beyond max? z2 = (limit-z1)/2; % bisect elseif (limit < -0.5) & (z2+z1 > z1*EXT) % extrapolation beyond limit z2 = z1*(EXT-1.0); % set to extrapolation limit elseif z2 < -z3*INT z2 = -z3*INT; elseif (limit > -0.5) & (z2 < (limit-z1)*(1.0-INT)) % too close to limit? z2 = (limit-z1)*(1.0-INT); end f3 = f2; d3 = d2; z3 = -z2; % set point 3 equal to point 2 z1 = z1 + z2; X = X + z2*s; % update current estimates [f2 df2] = eval(argstr); M = M - 1; i = i + (length<0); % count epochs?! d2 = df2'*s; end % end of line search if success % if line search succeeded f1 = f2; fX = [fX' f1]'; fprintf('%s %4i | Cost: %4.6e\r', S, i, f1); s = (df2'*df2-df1'*df2)/(df1'*df1)*s - df2; % Polack-Ribiere direction tmp = df1; df1 = df2; df2 = tmp; % swap derivatives d2 = df1'*s; if d2 > 0 % new slope must be negative s = -df1; % otherwise use steepest direction d2 = -s'*s; end z1 = z1 * min(RATIO, d1/(d2-realmin)); % slope ratio but max RATIO d1 = d2; ls_failed = 0; % this line search did not fail else X = X0; f1 = f0; df1 = df0; % restore point from before failed line search if ls_failed | i > abs(length) % line search failed twice in a row break; % or we ran out of time, so we give up end tmp = df1; df1 = df2; df2 = tmp; % swap derivatives s = -df1; % try steepest d1 = -s'*s; z1 = 1/(1-d1); ls_failed = 1; % this line search failed end if exist('OCTAVE_VERSION') fflush(stdout); end end fprintf('\n'); ```