function [x,f,exitflag,output] = minFunc(funObj,x0,options,varargin)

% minFunc(funObj,x0,options,varargin)

%

% Unconstrained optimizer using a line search strategy

%

% Uses an interface very similar to fminunc

%   (it doesn't support all of the optimization toolbox options,

%       but supports many other options).

%

% It computes descent directions using one of ('Method'):

%   - 'sd': Steepest Descent

%       (no previous information used, not recommended)

%   - 'csd': Cyclic Steepest Descent

%       (uses previous step length for a fixed length cycle)

%   - 'bb': Barzilai and Borwein Gradient

%       (uses only previous step)

%   - 'cg': Non-Linear Conjugate Gradient

%       (uses only previous step and a vector beta)

%   - 'scg': Scaled Non-Linear Conjugate Gradient

%       (uses previous step and a vector beta, 

%           and Hessian-vector products to initialize line search)

%   - 'pcg': Preconditionined Non-Linear Conjugate Gradient

%       (uses only previous step and a vector beta, preconditioned version)

%   - 'lbfgs': Quasi-Newton with Limited-Memory BFGS Updating

%       (default: uses a predetermined nunber of previous steps to form a 

%           low-rank Hessian approximation)

%   - 'newton0': Hessian-Free Newton

%       (numerically computes Hessian-Vector products)

%   - 'pnewton0': Preconditioned Hessian-Free Newton 

%       (numerically computes Hessian-Vector products, preconditioned

%       version)

%   - 'qnewton': Quasi-Newton Hessian approximation

%       (uses dense Hessian approximation)

%   - 'mnewton': Newton's method with Hessian calculation after every

%   user-specified number of iterations

%       (needs user-supplied Hessian matrix)

%   - 'newton': Newton's method with Hessian calculation every iteration

%       (needs user-supplied Hessian matrix)

%   - 'tensor': Tensor

%       (needs user-supplied Hessian matrix and Tensor of 3rd partial derivatives)

%

% Several line search strategies are available for finding a step length satisfying

%   the termination criteria ('LS'):

%   - 0: Backtrack w/ Step Size Halving

%   - 1: Backtrack w/ Quadratic/Cubic Interpolation from new function values

%   - 2: Backtrack w/ Cubic Interpolation from new function + gradient

%   values (default for 'bb' and 'sd')

%   - 3: Bracketing w/ Step Size Doubling and Bisection

%   - 4: Bracketing w/ Cubic Interpolation/Extrapolation with function +

%   gradient values (default for all except 'bb' and 'sd')

%   - 5: Bracketing w/ Mixed Quadratic/Cubic Interpolation/Extrapolation

%   - 6: Use Matlab Optimization Toolbox's line search

%           (requires Matlab's linesearch.m to be added to the path)

%

%   Above, the first three find a point satisfying the Armijo conditions,

%   while the last four search for find a point satisfying the Wolfe

%   conditions.  If the objective function overflows, it is recommended

%   to use one of the first 3.

%   The first three can be used to perform a non-monotone

%   linesearch by changing the option 'Fref'.

%

% Several strategies for choosing the initial step size are avaiable ('LS_init'):

%   - 0: Always try an initial step length of 1 (default for all except 'cg' and 'sd')

%       (t = 1)

%   - 1: Use a step similar to the previous step (default for 'cg' and 'sd')

%       (t = t_old*min(2,g'd/g_old'd_old))

%   - 2: Quadratic Initialization using previous function value and new

%   function value/gradient (use this if steps tend to be very long)

%       (t = min(1,2*(f-f_old)/g))

%   - 3: The minimum between 1 and twice the previous step length

%       (t = min(1,2*t)

%   - 4: The scaled conjugate gradient step length (may accelerate

%   conjugate gradient methods, but requires a Hessian-vector product)

%       (t = g'd/d'Hd)

%

% Inputs:

%   funObj is a function handle

%   x0 is a starting vector;

%   options is a struct containing parameters

%  (defaults are used for non-existent or blank fields)

%   all other arguments are passed to funObj

%

% Outputs:

%   x is the minimum value found

%   f is the function value at the minimum found

%   exitflag returns an exit condition

%   output returns a structure with other information

%

% Supported Input Options

%   Display - Level of display [ off | final | (iter) | full | excessive ]

%   MaxFunEvals - Maximum number of function evaluations allowed (1000)

%   MaxIter - Maximum number of iterations allowed (500)

%   TolFun - Termination tolerance on the first-order optimality (1e-5)

%   TolX - Termination tolerance on progress in terms of function/parameter changes (1e-9)

%   Method - [ sd | csd | bb | cg | scg | pcg | {lbfgs} | newton0 | pnewton0 |

%       qnewton | mnewton | newton | tensor ]

%   c1 - Sufficient Decrease for Armijo condition (1e-4)

%   c2 - Curvature Decrease for Wolfe conditions (.2 for cg methods, .9 otherwise)

%   LS_init - Line Search Initialization -see above (2 for cg/sd, 4 for scg, 0 otherwise)

%   LS - Line Search type -see above (2 for bb, 4 otherwise)

%   Fref - Setting this to a positive integer greater than 1

%       will use non-monotone Armijo objective in the line search.

%       (20 for bb, 10 for csd, 1 for all others)

%   numDiff - compute derivative numerically

%       (default: 0) (this option has a different effect for 'newton', see below)

%   useComplex - if 1, use complex differentials when computing numerical derivatives

%       to get very accurate values (default: 0)

%   DerivativeCheck - if 'on', computes derivatives numerically at initial

%       point and compares to user-supplied derivative (default: 'off')

%   outputFcn - function to run after each iteration (default: []).  It

%       should have the following interface:

%       outputFcn(x,infoStruct,state,varargin{:})

%   useMex - where applicable, use mex files to speed things up (default: 1)

%

% Method-specific input options:

%   newton:

%       HessianModify - type of Hessian modification for direct solvers to

%       use if the Hessian is not positive definite (default: 0)

%           0: Minimum Euclidean norm s.t. eigenvalues sufficiently large

%           (requires eigenvalues on iterations where matrix is not pd)

%           1: Start with (1/2)*||A||_F and increment until Cholesky succeeds

%           (an approximation to method 0, does not require eigenvalues)

%           2: Modified LDL factorization

%           (only 1 generalized Cholesky factorization done and no eigenvalues required)

%           3: Modified Spectral Decomposition

%           (requires eigenvalues)

%           4: Modified Symmetric Indefinite Factorization

%           5: Uses the eigenvector of the smallest eigenvalue as negative

%           curvature direction

%       cgSolve - use conjugate gradient instead of direct solver (default: 0)

%           0: Direct Solver

%           1: Conjugate Gradient

%           2: Conjugate Gradient with Diagonal Preconditioner

%           3: Conjugate Gradient with LBFGS Preconditioner

%           x: Conjugate Graident with Symmetric Successive Over Relaxation

%           Preconditioner with parameter x

%               (where x is a real number in the range [0,2])

%           x: Conjugate Gradient with Incomplete Cholesky Preconditioner

%           with drop tolerance -x

%               (where x is a real negative number)

%       numDiff - compute Hessian numerically

%                 (default: 0, done with complex differentials if useComplex = 1)

%       LS_saveHessiancomp - when on, only computes the Hessian at the

%       first and last iteration of the line search (default: 1)

%   mnewton:

%       HessianIter - number of iterations to use same Hessian (default: 5)

%   qnewton:

%       initialHessType - scale initial Hessian approximation (default: 1)

%       qnUpdate - type of quasi-Newton update (default: 3):

%           0: BFGS

%           1: SR1 (when it is positive-definite, otherwise BFGS)

%           2: Hoshino

%           3: Self-Scaling BFGS

%           4: Oren's Self-Scaling Variable Metric method 

%           5: McCormick-Huang asymmetric update

%       Damped - use damped BFGS update (default: 1)

%   newton0/pnewton0:

%       HvFunc - user-supplied function that returns Hessian-vector products

%           (by default, these are computed numerically using autoHv)

%           HvFunc should have the following interface: HvFunc(v,x,varargin{:})

%       useComplex - use a complex perturbation to get high accuracy

%           Hessian-vector products (default: 0)

%           (the increased accuracy can make the method much more efficient,

%               but gradient code must properly support complex inputs)

%       useNegCurv - a negative curvature direction is used as the descent

%           direction if one is encountered during the cg iterations

%           (default: 1)

%       precFunc (for pnewton0 only) - user-supplied preconditioner

%           (by default, an L-BFGS preconditioner is used)

%           precFunc should have the following interfact:

%           precFunc(v,x,varargin{:})

%   lbfgs:

%       Corr - number of corrections to store in memory (default: 100)

%           (higher numbers converge faster but use more memory)

%       Damped - use damped update (default: 0)

%   pcg:

%       cgUpdate - type of update (default: 2)

%   cg/scg/pcg:

%       cgUpdate - type of update (default for cg/scg: 2, default for pcg: 1)

%           0: Fletcher Reeves

%           1: Polak-Ribiere

%           2: Hestenes-Stiefel (not supported for pcg)

%           3: Gilbert-Nocedal

%       HvFunc (for scg only)- user-supplied function that returns Hessian-vector 

%           products

%           (by default, these are computed numerically using autoHv)

%           HvFunc should have the following interface:

%           HvFunc(v,x,varargin{:})

%       precFunc (for pcg only) - user-supplied preconditioner

%           (by default, an L-BFGS preconditioner is used)

%           precFunc should have the following interfact:

%           precFunc(v,x,varargin{:})

%   bb:

%       bbType - type of bb step (default: 1)

%           0: min_alpha ||delta_x - alpha delta_g||_2

%           1: min_alpha ||alpha delta_x - delta_g||_2

%           2: Conic BB

%           3: Gradient method with retards

%   csd:

%       cycle - length of cycle (default: 3)

%

% Supported Output Options

%   iterations - number of iterations taken

%   funcCount - number of function evaluations

%   algorithm - algorithm used

%   firstorderopt - first-order optimality

%   message - exit message

%   trace.funccount - function evaluations after each iteration

%   trace.fval - function value after each iteration

%

% Author: Mark Schmidt (2006)

% Web: http://www.cs.ubc.ca/~schmidtm

%

% Sources (in order of how much the source material contributes):

%   J. Nocedal and S.J. Wright.  1999.  "Numerical Optimization".  Springer Verlag.

%   R. Fletcher.  1987.  "Practical Methods of Optimization".  Wiley.

%   J. Demmel.  1997.  "Applied Linear Algebra.  SIAM.

%   R. Barret, M. Berry, T. Chan, J. Demmel, J. Dongarra, V. Eijkhout, R.

%   Pozo, C. Romine, and H. Van der Vost.  1994.  "Templates for the Solution of

%   Linear Systems: Building Blocks for Iterative Methods".  SIAM.

%   J. More and D. Thuente.  "Line search algorithms with guaranteed

%   sufficient decrease".  ACM Trans. Math. Softw. vol 20, 286-307, 1994.

%   M. Raydan.  "The Barzilai and Borwein gradient method for the large

%   scale unconstrained minimization problem".  SIAM J. Optim., 7, 26-33,

%   (1997).

%   "Mathematical Optimization".  The Computational Science Education

%   Project.  1995.

%   C. Kelley.  1999.  "Iterative Methods for Optimization".  Frontiers in

%   Applied Mathematics.  SIAM.

if nargin < 3

    options = [];

end

% Get Parameters

[verbose,verboseI,debug,doPlot,maxFunEvals,maxIter,tolFun,tolX,method,...

    corrections,c1,c2,LS_init,LS,cgSolve,qnUpdate,cgUpdate,initialHessType,...

    HessianModify,Fref,useComplex,numDiff,LS_saveHessianComp,...

    DerivativeCheck,Damped,HvFunc,bbType,cycle,...

    HessianIter,outputFcn,useMex,useNegCurv,precFunc] = ...

    minFunc_processInputOptions(options);

if isfield(options, 'logfile')

    logfile = options.logfile;

else

    logfile = [];

end

% Constants

SD = 0;

CSD = 1;

BB = 2;

CG = 3;

PCG = 4;

LBFGS = 5;

QNEWTON = 6;

NEWTON0 = 7;

NEWTON = 8;

TENSOR = 9;

% Initialize

p = length(x0);

d = zeros(p,1);

x = x0;

t = 1;

% If necessary, form numerical differentiation functions

funEvalMultiplier = 1;

if numDiff && method ~= TENSOR

    varargin(3:end+2) = varargin(1:end);

    varargin{1} = useComplex;

    varargin{2} = funObj;

    if method ~= NEWTON

        if debug

            if useComplex

                fprintf('Using complex differentials for gradient computation\n');

            else

                fprintf('Using finite differences for gradient computation\n');

            end

        end

        funObj = @autoGrad;

    else

        if debug

            if useComplex

                fprintf('Using complex differentials for gradient computation\n');

            else

                fprintf('Using finite differences for gradient computation\n');

            end

        end

        funObj = @autoHess;

    end

    if method == NEWTON0 && useComplex == 1

        if debug

            fprintf('Turning off the use of complex differentials\n');

        end

        useComplex = 0;

    end

    if useComplex

        funEvalMultiplier = p;

    else

        funEvalMultiplier = p+1;

    end

end

% Evaluate Initial Point

if method < NEWTON

    [f,g] = feval(funObj, x, varargin{:});

else

    [f,g,H] = feval(funObj, x, varargin{:});

    computeHessian = 1;

end

funEvals = 1;

if strcmp(DerivativeCheck,'on')

    if numDiff

        fprintf('Can not do derivative checking when numDiff is 1\n');

    end

    % Check provided gradient/hessian function using numerical derivatives

    fprintf('Checking Gradient:\n');

    [f2,g2] = autoGrad(x,useComplex,funObj,varargin{:});

    fprintf('Max difference between user and numerical gradient: %f\n',max(abs(g-g2)));

    if max(abs(g-g2)) > 1e-4

        fprintf('User NumDif:\n');

        [g g2]

        diff = abs(g-g2)

        pause;

    end

    if method >= NEWTON

        fprintf('Check Hessian:\n');

        [f2,g2,H2] = autoHess(x,useComplex,funObj,varargin{:});

        fprintf('Max difference between user and numerical hessian: %f\n',max(abs(H(:)-H2(:))));

        if max(abs(H(:)-H2(:))) > 1e-4

            H

            H2

            diff = abs(H-H2)

            pause;

        end

    end

end

% Output Log

if verboseI

    fprintf('%10s %10s %15s %15s %15s\n','Iteration','FunEvals','Step Length','Function Val','Opt Cond');

end

if logfile

    fid = fopen(logfile, 'a');

    if (fid > 0)

        fprintf(fid, '-- %10s %10s %15s %15s %15s\n','Iteration','FunEvals','Step Length','Function Val','Opt Cond');

        fclose(fid);

    end

end

% Output Function

if ~isempty(outputFcn)

    callOutput(outputFcn,x,'init',0,funEvals,f,[],[],g,[],sum(abs(g)),varargin{:});

end

% Initialize Trace

trace.fval = f;

trace.funcCount = funEvals;

% Check optimality of initial point

if sum(abs(g)) <= tolFun

    exitflag=1;

    msg = 'Optimality Condition below TolFun';

    if verbose

        fprintf('%s\n',msg);

    end

    if nargout > 3

        output = struct('iterations',0,'funcCount',1,...

            'algorithm',method,'firstorderopt',sum(abs(g)),'message',msg,'trace',trace);

    end

    return;

end

% Perform up to a maximum of 'maxIter' descent steps:

for i = 1:maxIter

    % ****************** COMPUTE DESCENT DIRECTION *****************

    switch method

        case SD % Steepest Descent

            d = -g;

        case CSD % Cyclic Steepest Descent

            if mod(i,cycle) == 1 % Use Steepest Descent

                alpha = 1;

                LS_init = 2;

                LS = 4; % Precise Line Search

            elseif mod(i,cycle) == mod(1+1,cycle) % Use Previous Step

                alpha = t;

                LS_init = 0;

                LS = 2; % Non-monotonic line search

            end

            d = -alpha*g;

        case BB % Steepest Descent with Barzilai and Borwein Step Length

            if i == 1

                d = -g;

            else

                y = g-g_old;

                s = t*d;

                if bbType == 0

                    yy = y'*y;

                    alpha = (s'*y)/(yy);

                    if alpha <= 1e-10 || alpha > 1e10

                        alpha = 1;

                    end

                elseif bbType == 1

                    sy = s'*y;

                    alpha = (s'*s)/sy;

                    if alpha <= 1e-10 || alpha > 1e10

                        alpha = 1;

                    end

                elseif bbType == 2 % Conic Interpolation ('Modified BB')

                    sy = s'*y;

                    ss = s'*s;

                    alpha = ss/sy;

                    if alpha <= 1e-10 || alpha > 1e10

                        alpha = 1;

                    end

                    alphaConic = ss/(6*(myF_old - f) + 4*g'*s + 2*g_old'*s);

                    if alphaConic > .001*alpha && alphaConic < 1000*alpha

                        alpha = alphaConic;

                    end

                elseif bbType == 3 % Gradient Method with retards (bb type 1, random selection of previous step)

                    sy = s'*y;

                    alpha = (s'*s)/sy;

                    if alpha <= 1e-10 || alpha > 1e10

                        alpha = 1;

                    end

                    v(1+mod(i-2,5)) = alpha;

                    alpha = v(ceil(rand*length(v)));

                end

                d = -alpha*g;

            end

            g_old = g;

            myF_old = f;

        case CG % Non-Linear Conjugate Gradient

            if i == 1

                d = -g; % Initially use steepest descent direction

            else

                gtgo = g'*g_old;

                gotgo = g_old'*g_old;

                if cgUpdate == 0

                    % Fletcher-Reeves

                    beta = (g'*g)/(gotgo);

                elseif cgUpdate == 1

                    % Polak-Ribiere

                    beta = (g'*(g-g_old)) /(gotgo);

                elseif cgUpdate == 2

                    % Hestenes-Stiefel

                    beta = (g'*(g-g_old))/((g-g_old)'*d);

                else

                    % Gilbert-Nocedal

                    beta_FR = (g'*(g-g_old)) /(gotgo);

                    beta_PR = (g'*g-gtgo)/(gotgo);

                    beta = max(-beta_FR,min(beta_PR,beta_FR));

                end

                d = -g + beta*d;

                % Restart if not a direction of sufficient descent

                if g'*d > -tolX

                    if debug

                        fprintf('Restarting CG\n');

                    end

                    beta = 0;

                    d = -g;

                end

                % Old restart rule:

                %if beta < 0 || abs(gtgo)/(gotgo) >= 0.1 || g'*d >= 0

            end

            g_old = g;

        case PCG % Preconditioned Non-Linear Conjugate Gradient

            % Apply preconditioner to negative gradient

            if isempty(precFunc)

                % Use L-BFGS Preconditioner

                if i == 1

                    old_dirs = zeros(length(g),0);

                    old_stps = zeros(length(g),0);

                    Hdiag = 1;

                    s = -g;

                else

                    [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);

                    if useMex

                        s = lbfgsC(-g,old_dirs,old_stps,Hdiag);

                    else

                        s = lbfgs(-g,old_dirs,old_stps,Hdiag);

                    end

                end

            else % User-supplied preconditioner

                s = precFunc(-g,x,varargin{:});

            end

            if i == 1

                d = s;

            else

                if cgUpdate == 0

                    % Preconditioned Fletcher-Reeves

                    beta = (g'*s)/(g_old'*s_old);

                elseif cgUpdate < 3

                    % Preconditioned Polak-Ribiere

                    beta = (g'*(s-s_old))/(g_old'*s_old);

                else

                    % Preconditioned Gilbert-Nocedal

                    beta_FR = (g'*s)/(g_old'*s_old);

                    beta_PR = (g'*(s-s_old))/(g_old'*s_old);

                    beta = max(-beta_FR,min(beta_PR,beta_FR));

                end

                d = s + beta*d;

                if g'*d > -tolX

                    if debug

                        fprintf('Restarting CG\n');

                    end

                    beta = 0;

                    d = s;

                end

            end

            g_old = g;

            s_old = s;

        case LBFGS % L-BFGS

            % Update the direction and step sizes

            if i == 1

                d = -g; % Initially use steepest descent direction

                old_dirs = zeros(length(g),0);

                old_stps = zeros(length(d),0);

                Hdiag = 1;

            else

                if Damped

                    [old_dirs,old_stps,Hdiag] = dampedUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);

                else

                    [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);

                end

                if useMex

                    d = lbfgsC(-g,old_dirs,old_stps,Hdiag);

                else

                    d = lbfgs(-g,old_dirs,old_stps,Hdiag);

                end

            end

            g_old = g;

        case QNEWTON % Use quasi-Newton Hessian approximation

            if i == 1

                d = -g;

            else

                % Compute difference vectors

                y = g-g_old;

                s = t*d;

                if i == 2

                    % Make initial Hessian approximation

                    if initialHessType == 0

                        % Identity

                        if qnUpdate <= 1

                            R = eye(length(g));

                        else

                            H = eye(length(g));

                        end

                    else

                        % Scaled Identity

                        if debug

                            fprintf('Scaling Initial Hessian Approximation\n');

                        end

                        if qnUpdate <= 1

                            % Use Cholesky of Hessian approximation

                            R = sqrt((y'*y)/(y'*s))*eye(length(g));

                        else

                            % Use Inverse of Hessian approximation

                            H = eye(length(g))*(y'*s)/(y'*y);

                        end

                    end

                end

                if qnUpdate == 0 % Use BFGS updates

                    Bs = R'*(R*s);

                    if Damped

                        eta = .02;

                        if y'*s < eta*s'*Bs

                            if debug

                                fprintf('Damped Update\n');

                            end

                            theta = min(max(0,((1-eta)*s'*Bs)/(s'*Bs - y'*s)),1);

                            y = theta*y + (1-theta)*Bs;

                        end

                        R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');

                    else

                        if y'*s > 1e-10

                            R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');

                        else

                            if debug

                                fprintf('Skipping Update\n');

                            end

                        end

                    end

                elseif qnUpdate == 1 % Perform SR1 Update if it maintains positive-definiteness

                    Bs = R'*(R*s);

                    ymBs = y-Bs;

                    if abs(s'*ymBs) >= norm(s)*norm(ymBs)*1e-8 && (s-((R\(R'\y))))'*y > 1e-10

                        R = cholupdate(R,-ymBs/sqrt(ymBs'*s),'-');

                    else

                        if debug

                            fprintf('SR1 not positive-definite, doing BFGS Update\n');

                        end

                        if Damped

                            eta = .02;

                            if y'*s < eta*s'*Bs

                                if debug

                                    fprintf('Damped Update\n');

                                end

                                theta = min(max(0,((1-eta)*s'*Bs)/(s'*Bs - y'*s)),1);

                                y = theta*y + (1-theta)*Bs;

                            end

                            R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');

                        else

                            if y'*s > 1e-10

                                R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');

                            else

                                if debug

                                    fprintf('Skipping Update\n');

                                end

                            end

                        end

                    end

                elseif qnUpdate == 2 % Use Hoshino update

                    v = sqrt(y'*H*y)*(s/(s'*y) - (H*y)/(y'*H*y));

                    phi = 1/(1 + (y'*H*y)/(s'*y));

                    H = H + (s*s')/(s'*y) - (H*y*y'*H)/(y'*H*y) + phi*v*v';

                elseif qnUpdate == 3 % Self-Scaling BFGS update

                    ys = y'*s;

                    Hy = H*y;

                    yHy = y'*Hy;

                    gamma = ys/yHy;

                    v = sqrt(yHy)*(s/ys - Hy/yHy);

                    H = gamma*(H - Hy*Hy'/yHy + v*v') + (s*s')/ys;

                elseif qnUpdate == 4 % Oren's Self-Scaling Variable Metric update

                    % Oren's method

                    if (s'*y)/(y'*H*y) > 1

                        phi = 1; % BFGS

                        omega = 0;

                    elseif (s'*(H\s))/(s'*y) < 1

                        phi = 0; % DFP

                        omega = 1;

                    else

                        phi = (s'*y)*(y'*H*y-s'*y)/((s'*(H\s))*(y'*H*y)-(s'*y)^2);

                        omega = phi;

                    end

                    gamma = (1-omega)*(s'*y)/(y'*H*y) + omega*(s'*(H\s))/(s'*y);

                    v = sqrt(y'*H*y)*(s/(s'*y) - (H*y)/(y'*H*y));

                    H = gamma*(H - (H*y*y'*H)/(y'*H*y) + phi*v*v') + (s*s')/(s'*y);

                elseif qnUpdate == 5 % McCormick-Huang asymmetric update

                    theta = 1;

                    phi = 0;

                    psi = 1;

                    omega = 0;

                    t1 = s*(theta*s + phi*H'*y)';

                    t2 = (theta*s + phi*H'*y)'*y;

                    t3 = H*y*(psi*s + omega*H'*y)';

                    t4 = (psi*s + omega*H'*y)'*y;

                    H = H + t1/t2 - t3/t4;

                end

                if qnUpdate <= 1

                    d = -R\(R'\g);

                else

                    d = -H*g;

                end

            end

            g_old = g;

        case NEWTON0 % Hessian-Free Newton

            cgMaxIter = min(p,maxFunEvals-funEvals);

            cgForce = min(0.5,sqrt(norm(g)))*norm(g);

            % Set-up preconditioner

            precondFunc = [];

            precondArgs = [];

            if cgSolve == 1

                if isempty(precFunc) % Apply L-BFGS preconditioner

                    if i == 1

                        old_dirs = zeros(length(g),0);

                        old_stps = zeros(length(g),0);

                        Hdiag = 1;

                    else

                        [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);

                        if useMex

                            precondFunc = @lbfgsC;

                        else

                            precondFunc = @lbfgs;

                        end

                        precondArgs = {old_dirs,old_stps,Hdiag};

                    end

                    g_old = g;

                else

                    % Apply user-defined preconditioner

                    precondFunc = precFunc;

                    precondArgs = {x,varargin{:}};

                end

            end

            % Solve Newton system using cg and hessian-vector products

            if isempty(HvFunc)

                % No user-supplied Hessian-vector function,

                % use automatic differentiation

                HvFun = @autoHv;

                HvArgs = {x,g,useComplex,funObj,varargin{:}};

            else

                % Use user-supplid Hessian-vector function

                HvFun = HvFunc;

                HvArgs = {x,varargin{:}};

            end

            if useNegCurv

                [d,cgIter,cgRes,negCurv] = conjGrad([],-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs,HvFun,HvArgs);

            else

                [d,cgIter,cgRes] = conjGrad([],-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs,HvFun,HvArgs);

            end

            funEvals = funEvals+cgIter;

            if debug

                fprintf('newtonCG stopped on iteration %d w/ residual %.5e\n',cgIter,cgRes);

            end

            if useNegCurv

                if ~isempty(negCurv)

                    %if debug

                    fprintf('Using negative curvature direction\n');

                    %end

                    d = negCurv/norm(negCurv);

                    d = d/sum(abs(g));

                end

            end

        case NEWTON % Newton search direction

            if cgSolve == 0

                if HessianModify == 0

                    % Attempt to perform a Cholesky factorization of the Hessian

                    [R,posDef] = chol(H);

                    % If the Cholesky factorization was successful, then the Hessian is

                    % positive definite, solve the system

                    if posDef == 0

                        d = -R\(R'\g);

                    else

                        % otherwise, adjust the Hessian to be positive definite based on the

                        % minimum eigenvalue, and solve with QR

                        % (expensive, we don't want to do this very much)

                        if debug

                            fprintf('Adjusting Hessian\n');

                        end

                        H = H + eye(length(g)) * max(0,1e-12 - min(real(eig(H))));

                        d = -H\g;

                    end

                elseif HessianModify == 1

                    % Modified Incomplete Cholesky

                    R = mcholinc(H,debug);

                    d = -R\(R'\g);

                elseif HessianModify == 2

                    % Modified Generalized Cholesky

                    if useMex

                        [L D perm] = mcholC(H);

                    else

                        [L D perm] = mchol(H);

                    end

                    d(perm) = -L' \ ((D.^-1).*(L \ g(perm)));

                elseif HessianModify == 3

                    % Modified Spectral Decomposition

                    [V,D] = eig((H+H')/2);

                    D = diag(D);

                    D = max(abs(D),max(max(abs(D)),1)*1e-12);

                    d = -V*((V'*g)./D);

                elseif HessianModify == 4

                    % Modified Symmetric Indefinite Factorization

                    [L,D,perm] = ldl(H,'vector');

                    [blockPos junk] = find(triu(D,1));

                    for diagInd = setdiff(setdiff(1:p,blockPos),blockPos+1)

                        if D(diagInd,diagInd) < 1e-12

                            D(diagInd,diagInd) = 1e-12;

                        end

                    end

                    for blockInd = blockPos'

                        block = D(blockInd:blockInd+1,blockInd:blockInd+1);

                        block_a = block(1);

                        block_b = block(2);

                        block_d = block(4);

                        lambda = (block_a+block_d)/2 - sqrt(4*block_b^2 + (block_a - block_d)^2)/2;

                        D(blockInd:blockInd+1,blockInd:blockInd+1) = block+eye(2)*(lambda+1e-12);

                    end

                    d(perm) = -L' \ (D \ (L \ g(perm)));

                else

                    % Take Newton step if Hessian is pd,

                    % otherwise take a step with negative curvature

                    [R,posDef] = chol(H);

                    if posDef == 0

                        d = -R\(R'\g);

                    else

                        if debug

                            fprintf('Taking Direction of Negative Curvature\n');

                        end

                        [V,D] = eig(H);

                        u = V(:,1);

                        d = -sign(u'*g)*u;

                    end

                end

            else

                % Solve with Conjugate Gradient

                cgMaxIter = p;

                cgForce = min(0.5,sqrt(norm(g)))*norm(g);

                % Select Preconditioner

                if cgSolve == 1

                    % No preconditioner

                    precondFunc = [];

                    precondArgs = [];

                elseif cgSolve == 2

                    % Diagonal preconditioner

                    precDiag = diag(H);

                    precDiag(precDiag < 1e-12) = 1e-12 - min(precDiag);

                    precondFunc = @precondDiag;

                    precondArgs = {precDiag.^-1};

                elseif cgSolve == 3

                    % L-BFGS preconditioner

                    if i == 1

                        old_dirs = zeros(length(g),0);

                        old_stps = zeros(length(g),0);

                        Hdiag = 1;

                    else

                        [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);

                    end

                    g_old = g;

                    if useMex

                        precondFunc = @lbfgsC;

                    else

                        precondFunc = @lbfgs;

                    end

                    precondArgs = {old_dirs,old_stps,Hdiag};

                elseif cgSolve > 0

                    % Symmetric Successive Overelaxation Preconditioner

                    omega = cgSolve;

                    D = diag(H);

                    D(D < 1e-12) = 1e-12 - min(D);

                    precDiag = (omega/(2-omega))*D.^-1;

                    precTriu = diag(D/omega) + triu(H,1);

                    precondFunc = @precondTriuDiag;

                    precondArgs = {precTriu,precDiag.^-1};

                else

                    % Incomplete Cholesky Preconditioner

                    opts.droptol = -cgSolve;

                    opts.rdiag = 1;

                    R = cholinc(sparse(H),opts);

                    if min(diag(R)) < 1e-12

                        R = cholinc(sparse(H + eye*(1e-12 - min(diag(R)))),opts);

                    end

                    precondFunc = @precondTriu;

                    precondArgs = {R};

                end

                % Run cg with the appropriate preconditioner

                if isempty(HvFunc)

                    % No user-supplied Hessian-vector function

                    [d,cgIter,cgRes] = conjGrad(H,-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs);

                else

                    % Use user-supplied Hessian-vector function

                    [d,cgIter,cgRes] = conjGrad(H,-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs,HvFunc,{x,varargin{:}});

                end

                if debug

                    fprintf('CG stopped after %d iterations w/ residual %.5e\n',cgIter,cgRes);

                    %funEvals = funEvals + cgIter;

                end

            end

        case TENSOR % Tensor Method

            if numDiff

                % Compute 3rd-order Tensor Numerically

                [junk1 junk2 junk3 T] = autoTensor(x,useComplex,funObj,varargin{:});

            else

                % Use user-supplied 3rd-derivative Tensor

                [junk1 junk2 junk3 T] = feval(funObj, x, varargin{:});

            end

            options_sub.Method = 'newton';

            options_sub.Display = 'none';

            options_sub.TolX = tolX;

            options_sub.TolFun = tolFun;

            d = minFunc(@taylorModel,zeros(p,1),options_sub,f,g,H,T);

            if any(abs(d) > 1e5) || all(abs(d) < 1e-5) || g'*d > -tolX

                if debug

                    fprintf('Using 2nd-Order Step\n');

                end

                [V,D] = eig((H+H')/2);

                D = diag(D);

                D = max(abs(D),max(max(abs(D)),1)*1e-12);

                d = -V*((V'*g)./D);

            else

                if debug

                    fprintf('Using 3rd-Order Step\n');

                end

            end

    end

    if ~isLegal(d)

        fprintf('Step direction is illegal!\n');

        pause;

        return

    end

    % ****************** COMPUTE STEP LENGTH ************************

    % Directional Derivative

    gtd = g'*d;

    % Check that progress can be made along direction

    if gtd > -tolX

        exitflag=2;

        msg = 'Directional Derivative below TolX';

        break;

    end

    % Select Initial Guess

    if i == 1

        if method < NEWTON0

            t = min(1,1/sum(abs(g)));

        else

            t = 1;

        end

    else

        if LS_init == 0

            % Newton step

            t = 1;

        elseif LS_init == 1

            % Close to previous step length

            t = t*min(2,(gtd_old)/(gtd));

        elseif LS_init == 2

            % Quadratic Initialization based on {f,g} and previous f

            t = min(1,2*(f-f_old)/(gtd));

        elseif LS_init == 3

            % Double previous step length

            t = min(1,t*2);

        elseif LS_init == 4

            % Scaled step length if possible

            if isempty(HvFunc)

                % No user-supplied Hessian-vector function,

                % use automatic differentiation

                dHd = d'*autoHv(d,x,g,0,funObj,varargin{:});

            else

                % Use user-supplid Hessian-vector function

                dHd = d'*HvFunc(d,x,varargin{:});

            end

            funEvals = funEvals + 1;

            if dHd > 0

                t = -gtd/(dHd);

            else

                t = min(1,2*(f-f_old)/(gtd));

            end

        end

        if t <= 0

            t = 1;

        end

    end

    f_old = f;

    gtd_old = gtd;

    % Compute reference fr if using non-monotone objective

    if Fref == 1

        fr = f;

    else

        if i == 1

            old_fvals = repmat(-inf,[Fref 1]);

        end

        if i <= Fref

            old_fvals(i) = f;

        else

            old_fvals = [old_fvals(2:end);f];

        end

        fr = max(old_fvals);

    end

    computeHessian = 0;

    if method >= NEWTON

        if HessianIter == 1

            computeHessian = 1;

        elseif i > 1 && mod(i-1,HessianIter) == 0

            computeHessian = 1;

        end

    end

    % Line Search

    f_old = f;

    if LS < 3 % Use Armijo Bactracking

        % Perform Backtracking line search

        if computeHessian

            [t,x,f,g,LSfunEvals,H] = ArmijoBacktrack(x,t,d,f,fr,g,gtd,c1,LS,tolX,debug,doPlot,LS_saveHessianComp,funObj,varargin{:});

        else

            [t,x,f,g,LSfunEvals] = ArmijoBacktrack(x,t,d,f,fr,g,gtd,c1,LS,tolX,debug,doPlot,1,funObj,varargin{:});

        end

        funEvals = funEvals + LSfunEvals;

    elseif LS < 6

        % Find Point satisfying Wolfe

        if computeHessian

            [t,f,g,LSfunEvals,H] = WolfeLineSearch(x,t,d,f,g,gtd,c1,c2,LS,25,tolX,debug,doPlot,LS_saveHessianComp,funObj,varargin{:});

        else

            [t,f,g,LSfunEvals] = WolfeLineSearch(x,t,d,f,g,gtd,c1,c2,LS,25,tolX,debug,doPlot,1,funObj,varargin{:});

        end

        funEvals = funEvals + LSfunEvals;

        x = x + t*d;

    else

        % Use Matlab optim toolbox line search

        [t,f_new,fPrime_new,g_new,LSexitFlag,LSiter]=...

            lineSearch({'fungrad',[],funObj},x,p,1,p,d,f,gtd,t,c1,c2,-inf,maxFunEvals-funEvals,...

            tolX,[],[],[],varargin{:});

        funEvals = funEvals + LSiter;

        if isempty(t)

            exitflag = -2;

            msg = 'Matlab LineSearch failed';

            break;

        end

        if method >= NEWTON

            [f_new,g_new,H] = funObj(x + t*d,varargin{:});

            funEvals = funEvals + 1;

        end

        x = x + t*d;

        f = f_new;

        g = g_new;

    end

    % Output iteration information

    if verboseI

        fprintf('%10d %10d %15.5e %15.5e %15.5e\n',i,funEvals*funEvalMultiplier,t,f,sum(abs(g)));

    end

    if logfile

        fid = fopen(logfile, 'a');

        if (fid > 0)

            fprintf(fid, '-- %10d %10d %15.5e %15.5e %15.5e\n',i,funEvals*funEvalMultiplier,t,f,sum(abs(g)));

            fclose(fid);

        end

    end

    % Output Function

    if ~isempty(outputFcn)

        callOutput(outputFcn,x,'iter',i,funEvals,f,t,gtd,g,d,sum(abs(g)),varargin{:});

    end

    % Update Trace

    trace.fval(end+1,1) = f;

    trace.funcCount(end+1,1) = funEvals;

    % Check Optimality Condition

    if sum(abs(g)) <= tolFun

        exitflag=1;

        msg = 'Optimality Condition below TolFun';

        break;

    end

    % ******************* Check for lack of progress *******************

    if sum(abs(t*d)) <= tolX

        exitflag=2;

        msg = 'Step Size below TolX';

        break;

    end

    if abs(f-f_old) < tolX

        exitflag=2;

        msg = 'Function Value changing by less than TolX';

        break;

    end

    % ******** Check for going over iteration/evaluation limit *******************

    if funEvals*funEvalMultiplier > maxFunEvals

        exitflag = 0;

        msg = 'Exceeded Maximum Number of Function Evaluations';

        break;

    end

    if i == maxIter

        exitflag = 0;

        msg='Exceeded Maximum Number of Iterations';

        break;

    end

end

if verbose

    fprintf('%s\n',msg);

end

if nargout > 3

    output = struct('iterations',i,'funcCount',funEvals*funEvalMultiplier,...

        'algorithm',method,'firstorderopt',sum(abs(g)),'message',msg,'trace',trace);

end

% Output Function

if ~isempty(outputFcn)

    callOutput(outputFcn,x,'done',i,funEvals,f,t,gtd,g,d,sum(abs(g)),varargin{:});