# Understanding MATLAB's fmincg optimization function

I'm researching numerical optimization. Recently I've come across a variant of a conjugate gradient method named fmincg.

The function is written in MATLAB and is used in the famous Andrew Ng's course on Machine Learning on Coursera. According to the copyright notice, fmincg was written by Carl Edward Rasmussen.

You can find the full text of fmincg on the MATLAB file exchange.

I'm trying to figure out the algorithm. The structure is quite straightforward as it closely follows the reference text I'm currently using (Numerical Optimization, Nocedal & Wright).

However I'm having a harder time figuring out the details of the line search procedure. Specifically, here's the code for the line search loop:

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


I would like to understand the following:

• before the loop, the point 3 is initialized as: f3 = f1; d3 = d1; z3 = -z1; and the comment says Save point 0. But why does he set z3 = -z1?
(in case someone's wondering, no, it's - apparently - not a mistake because 1. the algo works and 2. I've tried running it with z3 = z1 and it breaks)
• how the author is using the Wolfe conditions and its constants?
• interpolation / extrapolation. Where do the equations he's using come from? E.g.:

A = 6*(f2-f3)/z3+3*(d2+d3); B = 3*(f3-f2)-z3*(d3+2*d2); z2 = -d2*z3*z3/(B+sqrt(B*B-A*d2*z3*z3));

• what's going on in the if-else checks that involve the limit variable? Note that before entering the while 1, he set limit = -1

EDIT:
added one more question to the bullet list

• I think the Wolfe conditions are (the negation of) the loop test (f2 > f1+z1*RHO*d1) | (d2 > -SIG*d1), but I can't readily suggest anything for your other questions. – Kirill Jan 3 '17 at 17:23
• yes, that was the easy one :) What's not very clear to me is why he makes the other check elseif d2 > SIG*d1, success = 1; break; which I cannot relate to the literature – blackgreen Jan 3 '17 at 17:43
• Yeah :) I didn't see the second one, I think it's the curvature condition with absolute signs (en.wikipedia.org/wiki/…) (because d2 < -SIG*d1 from the first if-guard) and means it's right to terminate the line search. If I may say, whoever wrote this had a page full of notes in front of them, and the code isn't written to be read without those notes. – Kirill Jan 3 '17 at 17:51
• I think it's z3 = -z1 because it says later in a comment that "z3 is now relative to the location of z2", so point 1 is at coordinate $-z_1$ relative to point 2 at $z_1$. – Kirill Jan 3 '17 at 17:59
• That's a good point. It is actually checking the absolute value of d2 as part of the strong Wolfe-Powell condition. This should more or less solve question number 2! – blackgreen Jan 3 '17 at 18:00