function [xp, fvalp, funp, fcost, iarm, levmar_hist, nfail, nunew] ... = lev_mar_exec(f, x, fun, jac, xc, gc, ... iteration_data, nuold) % LEV_MAR_EXEC % v1.01 % ilev counter incremented as it should have been in the first place % 1/10/2014 % % Example of the executive_function option. % %function [xp, fvalp, funp, fcost, iarm, levmar_hist, nfail, nunew] ... % = lev_mar_exec(f, x, fun, jac, xc, gc, ... % iteration_data, nuold) % % I/O follows the instructions in the manual/book. % % Get the options you need from iteration_data. % options=iteration_data.options; h=iteration_data.h; complete_history=options.complete_history; core_data=iteration_data.core_data; parallel=options.parallel; least_squares=options.least_squares; if least_squares == 0 disp('lev_mar_exec solves nonlinear least squares problems.') disp('Please set the least_squares option to 1') disp('and make sure you have a least squares problem to solve.') end % nfail=1; nu = nuold; % % Now it's a plain vanilla Levenberg-Marquardt step. % iarm=0; ilev=0; ilevmax=5; mf=length(fun); funmat=[]; failvec=[]; xarray=[]; fcost=0; ared=-1; xp=x; funp=fun; n=length(x); fvalp=fun'*fun/2; % % Gt the epsilon-inactive indices. % epsb=1.d-6; ilist=max(( x > epsb) , (x < 1 - epsb)); pr=diag(ilist); rjac=jac*pr; grad=jac'*fun; fval=fun'*fun/2; sbinding=-(eye(n) - pr)*grad; % while ilev < ilevmax & ared < 0 lmm=(nu*eye(n) + jac'*jac); % % Compute Levenberg-Marquardt step as per Section 3.9.4. % big_matrix=[jac*pr ; sqrt(nu)*pr]; [mr,nr]=size(pr); big_rhs=[fun; zeros(mr,1)]; lm_direction=-big_matrix\big_rhs; xt= x + lm_direction; xt = min(max(xt,0),1); step = xt-x; % % Evaluate function at trial point. % [funt, ifail, icount]=feval(f,xt,h,core_data); % % Get the data organized for the complete_history update. % funmat=[funmat, funt]; failvec=[failvec, ifail]; xarray=[xarray, xt]; ilev = ilev + 1; % % If the function returns a value, then % adjust the Levenberg parameter using the trial step. % if ifail == 0 % fvalt=.5*funt'*funt; fcost=fcost+icount; ared=fval-fvalt; pred = -step'*grad - step'*lmm*step/2; if ared < 0 nu = nu*2; iarm=iarm+1; else xp=xt; funp=funt; fvalp=fvalt; nfail=0; if ared/pred > .75 nu=nu/2; end if ared/pred < .25 nu=2*nu; end break; end else nu = nu*2; iarm=iarm+1; end end % % Clean up after yourself. % levmar_hist=build_history(xarray,funmat,failvec); % nunew=nu;