From 5425245ec138ac3868b2327803cfaf2e85d7ed80 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Mon, 12 May 2014 17:15:49 +0200 Subject: [PATCH] adding lmmcp --- matlab/dynare_solve.m | 9 + matlab/lmmcp/catstruct.m | 148 +++++ matlab/lmmcp/dyn_lmmcp.m | 77 +++ matlab/lmmcp/dyn_lmmcp_func.m | 56 ++ matlab/lmmcp/get_complementarity_conditions.m | 53 ++ matlab/lmmcp/lmmcp.m | 615 ++++++++++++++++++ 6 files changed, 958 insertions(+) create mode 100644 matlab/lmmcp/catstruct.m create mode 100644 matlab/lmmcp/dyn_lmmcp.m create mode 100644 matlab/lmmcp/dyn_lmmcp_func.m create mode 100644 matlab/lmmcp/get_complementarity_conditions.m create mode 100644 matlab/lmmcp/lmmcp.m diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m index 3033f1143..1b68c688b 100644 --- a/matlab/dynare_solve.m +++ b/matlab/dynare_solve.m @@ -164,6 +164,15 @@ elseif options_.solve_algo == 3 else [x,info] = csolve(func,x,[],1e-6,500,varargin{:}); end +elseif options_.solve_algo == 10 + olmmcp = options_.lmmcp; + olmmcp.Display = 'iter'; + [x,fval,exitflag] = lmmcp(func,x,olmmcp.lb,olmmcp.ub,olmmcp,varargin{:}); + if exitflag == 1 + info = 0; + else + info = 1; + end else error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9]') end diff --git a/matlab/lmmcp/catstruct.m b/matlab/lmmcp/catstruct.m new file mode 100644 index 000000000..8248366e4 --- /dev/null +++ b/matlab/lmmcp/catstruct.m @@ -0,0 +1,148 @@ +function A = catstruct(varargin) +% CATSTRUCT Concatenate or merge structures with different fieldnames +% X = CATSTRUCT(S1,S2,S3,...) merges the structures S1, S2, S3 ... +% into one new structure X. X contains all fields present in the various +% structures. An example: +% +% A.name = 'Me' ; +% B.income = 99999 ; +% X = catstruct(A,B) +% % -> X.name = 'Me' ; +% % X.income = 99999 ; +% +% If a fieldname is not unique among structures (i.e., a fieldname is +% present in more than one structure), only the value from the last +% structure with this field is used. In this case, the fields are +% alphabetically sorted. A warning is issued as well. An axample: +% +% S1.name = 'Me' ; +% S2.age = 20 ; S3.age = 30 ; S4.age = 40 ; +% S5.honest = false ; +% Y = catstruct(S1,S2,S3,S4,S5) % use value from S4 +% +% The inputs can be array of structures. All structures should have the +% same size. An example: +% +% C(1).bb = 1 ; C(2).bb = 2 ; +% D(1).aa = 3 ; D(2).aa = 4 ; +% CD = catstruct(C,D) % CD is a 1x2 structure array with fields bb and aa +% +% The last input can be the string 'sorted'. In this case, +% CATSTRUCT(S1,S2, ..., 'sorted') will sort the fieldnames alphabetically. +% To sort the fieldnames of a structure A, you could use +% CATSTRUCT(A,'sorted') but I recommend ORDERFIELDS for doing that. +% +% When there is nothing to concatenate, the result will be an empty +% struct (0x0 struct array with no fields). +% +% NOTE: To concatenate similar arrays of structs, you can use simple +% concatenation: +% A = dir('*.mat') ; B = dir('*.m') ; C = [A ; B] ; +% +% See also CAT, STRUCT, FIELDNAMES, STRUCT2CELL, ORDERFIELDS + +% for Matlab R13 and up +% version 3.0 (mar 2013) +% (c) Jos van der Geest +% email: jos@jasen.nl + +% (C) 2013 Christophe Gouel + +% History +% Created in 2005 +% Revisions +% 2.0 (sep 2007) removed bug when dealing with fields containing cell +% arrays (Thanks to Rene Willemink) +% 2.1 (sep 2008) added warning and error identifiers +% 2.2 (oct 2008) fixed error when dealing with empty structs (Thanks to +% Lars Barring) +% 3.0 (mar 2013) fixed problem when the inputs were array of structures +% (thanks to Tor Inge Birkenes for pointing this out). +% Rephrased the help section as well. + +error(nargchk(1,Inf,nargin)) ; +N = nargin ; + +if ~isstruct(varargin{end}), + if isequal(varargin{end},'sorted'), + sorted = 1 ; + N = N-1 ; + error(nargchk(1,Inf,N)) ; + else + error('catstruct:InvalidArgument','Last argument should be a structure, or the string "sorted".') ; + end +else + sorted = 0 ; +end + +sz0 = [] ; % used to check that all inputs have the same size + +% used to check for a few trivial cases +NonEmptyInputs = false(N,1) ; +NonEmptyInputsN = 0 ; + +% used to collect the fieldnames and the inputs +FN = cell(N,1) ; +VAL = cell(N,1) ; + +% parse the inputs +for ii=1:N, + X = varargin{ii} ; + if ~isstruct(X), + error('catstruct:InvalidArgument',['Argument #' num2str(ii) ' is not a structure.']) ; + end + + if ~isempty(X), + % empty structs are ignored + if ii > 1 && ~isempty(sz0) + if ~isequal(size(X), sz0) + error('catstruct:UnequalSizes','All structures should have the same size.') ; + end + else + sz0 = size(X) ; + end + NonEmptyInputsN = NonEmptyInputsN + 1 ; + NonEmptyInputs(ii) = true ; + FN{ii} = fieldnames(X) ; + VAL{ii} = struct2cell(X) ; + end +end + +if NonEmptyInputsN == 0 + % all structures were empty + A = struct([]) ; +elseif NonEmptyInputsN == 1, + % there was only one non-empty structure + A = varargin{NonEmptyInputs} ; + if sorted, + A = orderfields(A) ; + end +else + % there is actually something to concatenate + FN = cat(1,FN{:}) ; + VAL = cat(1,VAL{:}) ; + FN = squeeze(FN) ; + VAL = squeeze(VAL) ; + MatlabVersion = version; + if str2double(MatlabVersion(end-5:end-2))<2013 % Equivalent to, but faster than if verLessThan('matlab','8.1') + [UFN,ind] = unique(FN) ; + else + [UFN,ind] = unique(FN,'legacy') ; + end + + if numel(UFN) ~= numel(FN), + warning('catstruct:DuplicatesFound','Fieldnames are not unique between structures.') ; + sorted = 1 ; + end + + if sorted, + VAL = VAL(ind,:) ; + FN = FN(ind,:) ; + end + + A = cell2struct(VAL, FN); + A = reshape(A, sz0) ; % reshape into original format +end + + + diff --git a/matlab/lmmcp/dyn_lmmcp.m b/matlab/lmmcp/dyn_lmmcp.m new file mode 100644 index 000000000..ce3db1c7b --- /dev/null +++ b/matlab/lmmcp/dyn_lmmcp.m @@ -0,0 +1,77 @@ +function [endo_simul,info] = dyn_lmmcp(M,options,oo) + +% Copyright (C) 2014 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +[lb,ub,eq_index] = get_complementarity_conditions(M); + +lead_lag_incidence = M.lead_lag_incidence; + +ny = M.endo_nbr; + +max_lag = M.maximum_endo_lag; + +nyp = nnz(lead_lag_incidence(1,:)) ; +iyp = find(lead_lag_incidence(1,:)>0) ; +ny0 = nnz(lead_lag_incidence(2,:)) ; +iy0 = find(lead_lag_incidence(2,:)>0) ; +nyf = nnz(lead_lag_incidence(3,:)) ; +iyf = find(lead_lag_incidence(3,:)>0) ; + +nd = nyp+ny0+nyf; +nrc = nyf+1 ; +isp = [1:nyp] ; +is = [nyp+1:ny+nyp] ; +isf = iyf+nyp ; +isf1 = [nyp+ny+1:nyf+nyp+ny+1] ; +stop = 0 ; +iz = [1:ny+nyp+nyf]; + +periods = options.periods; +steady_state = oo.steady_state; +params = M.params; +endo_simul = oo.endo_simul; +exo_simul = oo.exo_simul; +i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)'); +i_cols_A1 = find(lead_lag_incidence(2:3,:)'); +i_cols_T = nonzeros(lead_lag_incidence(1:2,:)'); +i_cols_j = 1:nd; +i_upd = ny+(1:periods*ny); + +x = endo_simul(:); + +model_dynamic = str2func([M.fname,'_dynamic']); +z = x(find(lead_lag_incidence')); +[res,A] = model_dynamic(z, exo_simul, params, steady_state,2); +nnzA = nnz(A); + +LB = repmat(lb,periods,1); +UB = repmat(ub,periods,1); + +Y0 = endo_simul(:,1); +YT = endo_simul(:,end); +x = endo_simul(:,2:end-1); +x = x(:); + +func_handle = @(x) dyn_lmmcp_func(x,model_dynamic, Y0, YT, exo_simul, ... + params, steady_state, periods, ny, ... + lead_lag_incidence, i_cols_A1, i_cols_1, ... + i_cols_T, i_cols_j,nnzA,eq_index); + +[x, info] = lmmcp(func_handle,x,LB,UB); + +endo_simul = [Y0 reshape(x,ny,periods) YT]; \ No newline at end of file diff --git a/matlab/lmmcp/dyn_lmmcp_func.m b/matlab/lmmcp/dyn_lmmcp_func.m new file mode 100644 index 000000000..936f98306 --- /dev/null +++ b/matlab/lmmcp/dyn_lmmcp_func.m @@ -0,0 +1,56 @@ +function [F,A] = dyn_lmmcp_func(x, model_dynamic, Y0, YT, exo_simul, params, ... + steady_state, periods, ny, lead_lag_incidence, ... + i_cols_A1, i_cols_1, i_cols_T, i_cols_j, ... + nnzA,eq_index) + +% Copyright (C) 2014 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + + Y = [Y0; x; YT]; + + F = zeros(periods*ny,1); + if nargout == 2 + A = sparse([],[],[],periods*ny,periods*ny,periods*nnzA); + end + + i_rows = 1:ny; + i_cols = find(lead_lag_incidence'); + i_cols_A = i_cols; + + for it = 2:(periods+1) + + [res,jacobian] = model_dynamic(Y(i_cols),exo_simul, params, ... + steady_state,it); + F(i_rows) = res(eq_index); + + if nargout == 2 + if it == 2 + A(i_rows,i_cols_A1) = jacobian(eq_index,i_cols_1); + elseif it == periods+1 + A(i_rows,i_cols_A(i_cols_T)) = jacobian(eq_index,i_cols_T); + else + A(i_rows,i_cols_A) = jacobian(eq_index,i_cols_j); + end + end + + i_rows = i_rows + ny; + i_cols = i_cols + ny; + if nargout == 2 && it > 2 + i_cols_A = i_cols_A + ny; + end + end + diff --git a/matlab/lmmcp/get_complementarity_conditions.m b/matlab/lmmcp/get_complementarity_conditions.m new file mode 100644 index 000000000..a58f8e337 --- /dev/null +++ b/matlab/lmmcp/get_complementarity_conditions.m @@ -0,0 +1,53 @@ +function [lb,ub,eq_index] = get_complementarity_condition(M) + +% Copyright (C) 2014 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + + +etags = M.equations_tags; +ub = inf(M.endo_nbr,1); +lb = -ub; +eq_index = (1:M.endo_nbr)'; +for i=1:size(etags,1) + if strcmp(etags{i,2},'mcp') + [b,m] = strsplit(etags{i,3},{' ','<','>'}); + if length(m) == 1 + if any(m{1} == '<') + k = find(strcmp(b{1},cellstr(M.endo_names))); + if isempty(k) + error(sprintf(['Complementarity condition %s: variable %s is ' ... + 'not recognized',etags{i,3},b{1}])) + end + ub(k) = str2num(b{2}); + eq_index(etags{i,1}) = k; + eq_index(k) = etags{i,1}; + elseif any(m{1} == '>') + k = find(strcmp(b{1},cellstr(M.endo_names))); + if isempty(k) + error(sprintf(['Complementarity condition %s: variable %s is ' ... + 'not recognized',etags{i},b{1}])) + end + lb(k) = str2num(b{2}); + eq_index(etags{i,1}) = k; + eq_index(k) = etags{i,1}; + end + else + error(sprintf('Complementarity condition %s can''t be parsed',etags{i,3})) + end + end +end + diff --git a/matlab/lmmcp/lmmcp.m b/matlab/lmmcp/lmmcp.m new file mode 100644 index 000000000..1122b83ed --- /dev/null +++ b/matlab/lmmcp/lmmcp.m @@ -0,0 +1,615 @@ +function [x,FVAL,EXITFLAG,OUTPUT,JACOB] = lmmcp(FUN,x,lb,ub,options,varargin) +% LMMCP solves a mixed complementarity problem. +% +% LMMCP uses a semismooth least squares formulation. The method applies a +% Levenberg-Marquardt/Gauss-Newton algorithm to a least-squares formulation. +% +% X = LMMCP(FUN,X0) tries to solve the system of nonlinear equations F(X)=0 and +% starts at the vector X0. FUN accepts a vector X and return a vector F of equation +% values F evaluated at X and, as second output if required, a matrix J, the +% Jacobian evaluated at X. +% +% X = LMMCP(FUN,X0,LB,UB) solves the mixed complementarity problem of the form: +% LB =X => F(X)>0, +% LB<=X<=UB => F(X)=0, +% X =UB => F(X)<0. +% +% X = LMMCP(FUN,X0,LB,UB,OPTIONS) solves the MCP problem using the options +% defined in the structure OPTIONS. Main fields are +% Display : control the display of iterations, 'none' (default), +% 'iter-detailed' or 'final-detailed' +% Switch from phase I to phase II +% preprocess : activate preprocessor for phase I (default = 1) +% presteps : number of iterations in phase I (default = 20) +% Termination parameters +% epsilon2 : termination value of the merit function (default = 1E-16) +% MaxIter : maximum number of iterations (default = 500) +% tmin : safeguard stepsize (default = 1E-12) +% Stepsize parameters +% m : number of previous function values to use in the nonmonotone +% line search rule (default = 10) +% kwatch : maximum number of steps (default = 20 and should not be +% smaller than m) +% watchdog : activate the watchdog strategy (default = 1) +% Ther are other minor parameters. Please see the code for their default values +% and interpretation. +% +% [X,FVAL] = LMMCP(FUN,X0,...) returns the value of the equations FUN at X. +% +% [X,FVAL,EXITFLAG] = LMMCP(FUN,X0,...) returns EXITFLAG that describes the exit +% conditions. Possible values are +% 1 : LMMCP converged to a root +% 0 : Too many iterations +% -1 : +% +% [X,FVAL,EXITFLAG,OUTPUT] = LMMCP(FUN,X0,...) returns the structure OUTPUT that +% contains the number of iterations (OUTPUT.iterations), the value of the merit +% function (OUTPUT.Psix), and the norm of the derivative of the merit function +% (OUTPUT.normDPsix). +% +% [X,FVAL,EXITFLAG,OUTPUT,JACOB] = LMMCP(FUN,X0,...) returns JACOB the Jacobian +% of FUN evaluated at X. +% +% More details of the main program may be found in the following paper: +% +% Christian Kanzow and Stefania Petra: On a semismooth least squares formulation of +% complementarity problems with gap reduction. Optimization Methods and Software +% 19, 2004, pp. 507-525. +% +% In addition, the current implementation uses a preprocessor which is the +% projected Levenberg-Marquardt step from the following preprint: +% +% Christian Kanzow and Stefania Petra: Projected filter trust region methods for a +% semismooth least squares formulation of mixed complementarity +% problems. Optimization Methods and Software +% 22, 2007, pp. 713-735. +% +% A user's guide is also available: +% +% Christian Kanzow and Stefania Petra (2005). +% LMMCP --- A Levenberg-Marquardt-type MATLAB Solver for Mixed Complementarity Problems. +% University of Wuerzburg. +% http://www.mathematik.uni-wuerzburg.de/~kanzow/software/UserGuide.pdf +% +% This is a modification by Christophe Gouel of the original files, which can be +% downloaded from: +% http://www.mathematik.uni-wuerzburg.de/~kanzow/software/LMMCP.zip +% +% copyright: Christian Kanzow and Stefania Petra +% Institute of Applied Mathematics and Statistics +% University of Wuerzburg +% Am Hubland +% 97074 Wuerzburg +% GERMANY +% +% e-mail: kanzow@mathematik.uni-wuerzburg.de +% petra@mathematik.uni-wuerzburg.de + +% ------------Initialization---------------- +defaultopt = struct(... + 'beta', 0.55,... + 'Big', 1e10,... + 'delta', 5,... + 'deltamin', 1,... + 'Display', 'none',... + 'epsilon1', 1e-6,... + 'epsilon2', 1e-16,... + 'eta', 0.95,... + 'kwatch', 20,... + 'lambda1', 0.1,... + 'm', 10,... + 'MaxIter', 500,... + 'null', 1e-10,... + 'preprocess', 1,... + 'presteps', 20,... + 'sigma', 1e-4,... + 'sigma1', 0.5,... + 'sigma2', 2,... + 'tmin', 1e-12,... + 'watchdog', 1); + +if nargin < 4 + ub = inf(size(x)); + if nargin < 3 + lb = -inf(size(x)); + end +end + +if nargin < 5 || isempty(options) || ~isstruct(options) + options = defaultopt; +else + warning('off','catstruct:DuplicatesFound') + options = catstruct(defaultopt,options); +end + + +switch options.Display + case {'off','none'} + verbosity = 0; + case {'iter','iter-detailed'} + verbosity = 2; + case {'final','final-detailed'} + verbosity = 1; + otherwise + verbosity = 0; +end + +% parameter settings +eps1 = options.epsilon1; +eps2 = options.epsilon2; +null = options.null; +Big = options.Big; + +% maximal number of iterations +kmax = options.MaxIter; + +% choice of lambda +lambda1 = options.lambda1; +lambda2 = 1-lambda1; + +% steplength parameters +beta = options.beta; +sigma = options.sigma; +tmin = options.tmin; + +% parameters watchdog and nonmonotone line search; redefined later +m = options.m; +kwatch = options.kwatch; +watchdog = options.watchdog; % 1=watchdog strategy active, otherwise not + +% parameters for preprocessor +preprocess = options.preprocess; % 1=preprocessor used, otherwise not +presteps = options.presteps; % maximum number of preprocessing steps + +% trust-region parameters for preprocessor +delta = options.delta; +deltamin = options.deltamin; +sigma1 = options.sigma1; +sigma2 = options.sigma2; +eta = options.eta; + +% initializations +k = 0; +k_main = 0; + +% compute a feasible starting point by projection +x = max(lb,min(x,ub)); + +n = length(x); +OUTPUT.Dim = n; + +% definition of index sets I_l, I_u and I_lu +Indexset = zeros(n,1); +I_l = lb>-Big & ub>Big; +I_u = lb<-Big & ub-Big & ub 1 + fprintf(' k Psi(x) || DPsi(x) || stepsize\n'); + disp('====================================================================') + disp('********************* Output at starting point *********************') + fprintf('%4.0f %24.5e %24.5e\n',k,Psix,normDPsix); +end + +% +% Preprocessor using local method +% + +if preprocess==1 + + if verbosity > 1 + disp('************************** Preprocessor ****************************') + end + + normpLM=1; + while (k < presteps) && (Psix > eps2) && (normpLM>null) + k=k+1; + + % choice of Levenberg-Marquardt parameter, note that we do not use + % the condition estimator for large-scale problems, although this + % may cause numerical problems in some examples + + i=0; + mu=0; + if n<100000 + i=1; + mu=1e-16; + if condest(DPhix'*DPhix)>1e25 + mu=1e-6/(k+1); + end + end + if i==1 + pLM= [ DPhix ; sqrt(mu)*speye(n)]\[-Phix;sparse(n,1)]; + else + pLM=-DPhix\Phix; + end + normpLM=norm(pLM); + + % compute the projected Levenberg-Marquard step onto box Xk + lbnew=max(min(lb-x,0),-delta); + ubnew=min(max(ub-x,0),delta); + d=max(lbnew,min(pLM,ubnew)); + xnew=x+d; + + % function evaluations etc. + [Fxnew,DFxnew] = feval(FUN,xnew,varargin{:}); + Phixnew=Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset); + Psixnew=0.5*(Phixnew'*Phixnew); + normPhixnew=norm(Phixnew); + + % update of delta + if normPhixnew<=eta*normPhix + delta=max(deltamin,sigma2*delta); + elseif normPhixnew>5*eta*normPhix + delta=max(deltamin,sigma1*delta); + end + + % update + x=xnew; + Fx=Fxnew; + DFx=DFxnew; + Phix=Phixnew; + Psix=Psixnew; + normPhix=normPhixnew; + DPhix=DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset); + DPsix=DPhix'*Phix; + normDPsix=norm(DPsix,inf); + + % output at each iteration + t=1; + if verbosity > 1 + fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t); + end + end +end + +% terminate program or redefine current iterate as original initial point +if preprocess==1 && Psix 0 + fprintf('Psix = %1.4e\nnormDPsix = %1.4e\n',Psix,normDPsix); + disp('Approximate solution found.') + end + EXITFLAG = 1; + FVAL = Fx; + OUTPUT.iterations = k; + OUTPUT.Psix = Psix; + OUTPUT.normDPsix = normDPsix; + JACOB = DFx; + return +elseif preprocess==1 && Psix>=eps2 + x=x0; + Phix=Phix0; + Psix=Psix0; + DPhix=DPhix0; + DPsix=DPsix0; + if verbosity > 1 + disp('******************** Restart with initial point ********************') + fprintf('%4.0f %24.5e %24.5e\n',k_main,Psix0,normDPsix0); + end +end + +% +% Main algorithm +% + +if verbosity > 1 + disp('************************** Main program ****************************') +end + +while (k < kmax) && (Psix > eps2) + + % choice of Levenberg-Marquardt parameter, note that we do not use + % the condition estimator for large-scale problems, although this + % may cause numerical problems in some examples + + i=0; + if n<100 + i=1; + mu=1e-16; + if condest(DPhix'*DPhix)>1e25 + mu=1e-1/(k+1); + end + end + + % compute a Levenberg-Marquard direction + + if i==1 + d= [ DPhix ; sqrt(mu)*eye(n)]\[-Phix;zeros(n,1)]; + else + d=-DPhix\Phix; + end + + % computation of steplength t using the nonmonotone Armijo-rule + % starting with the 6-th iteration + + % computation of steplength t using the monotone Armijo-rule if + % d is a 'good' descent direction or k<=5 + + t = 1; + xnew = x+d; + Fxnew = feval(FUN,xnew,varargin{:}); + Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset); + Psixnew = 0.5*(Phixnew'*Phixnew); + const = sigma*DPsix'*d; + + while (Psixnew > MaxPsi + const*t) && (t > tmin) + t = t*beta; + xnew = x+t*d; + Fxnew = feval(FUN,xnew,varargin{:}); + Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset); + Psixnew = 0.5*(Phixnew'*Phixnew); + end + + % updatings + x=xnew; + Fx=Fxnew; + Phix=Phixnew; + Psix=Psixnew; + [~,DFx]=feval(FUN,x,varargin{:}); + DPhix=DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset); + DPsix=DPhix'*Phix; + normDPsix=norm(DPsix); + k=k+1; + k_main=k_main+1; + + if k_main<=5 + aux(mod(k_main,m)+1)=Psix; + MaxPsi=Psix; + else + aux(mod(k_main,m)+1) = Psix; + MaxPsi = max(aux); + end; + + % updatings for the watchdog strategy + if watchdog ==1 + if Psixkwatch + x=xbest; + Phix=Phibest; + Psix=Psibest; + DPhix=DPhibest; + DPsix=DPsibest; + normDPsix=normDPsibest; + MaxPsi=Psix; + end; + end; + + if verbosity > 1 + % output at each iteration + fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t); + end +end; + +% final output +if Psix<=eps2 + EXITFLAG = 1; + if verbosity > 0, disp('Approximate solution found.'); end +elseif k>=kmax + EXITFLAG = 0; + if verbosity > 0, disp('Maximum iteration number reached.'); end +elseif normDPsix<=eps1 + EXITFLAG = -1; % Provisoire + if verbosity > 0, disp('Approximate stationary point found.'); end +else + EXITFLAG = -1; % Provisoire + if verbosity > 0, disp('No solution found.'); end +end + +FVAL = Fx; +OUTPUT.iterations = k; +OUTPUT.Psix = Psix; +OUTPUT.normDPsix = normDPsix; +JACOB = DFx; + +% Subfunctions + +function y = Phi(x,Fx,lb,ub,lambda1,lambda2,n,Indexset) + +y = zeros(2*n,1); +phi_u = sqrt((ub-x).^2+Fx.^2)-ub+x+Fx; +LZ = false(n,1); % logical zero + +I0 = Indexset==0; +y(I0) = -lambda1*Fx(I0); +y([LZ; I0]) = -lambda2*Fx(I0); + +I1 = Indexset==1; +y(I1) = lambda1*(-x(I1)+lb(I1)-Fx(I1)+sqrt((x(I1)-lb(I1)).^2+Fx(I1).^2)); +y([LZ; I1]) = lambda2*max(0,x(I1)-lb(I1)).*max(0,Fx(I1)); + +I2 = Indexset==2; +y(I2) = -lambda1*phi_u(I2); +y([LZ; I2]) = lambda2*max(0,ub(I2)-x(I2)).*max(0,-Fx(I2)); + +I3 = Indexset==3; +y(I3) = lambda1*(sqrt((x(I3)-lb(I3)).^2+phi_u(I3).^2)-x(I3)+lb(I3)-phi_u(I3)); +y([LZ; I3]) = lambda2*(max(0,x(I3)-lb(I3)).*max(0,Fx(I3))+max(0,ub(I3)-x(I3)).*max(0,-Fx(I3))); + + +function H = DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset) +% DPHI evaluates an element of the C-subdifferential of operator Phi + +null = 1e-8; +beta_l = zeros(n,1); +beta_u = zeros(n,1); +alpha_l = zeros(n,1); +alpha_u = zeros(n,1); + + +z = zeros(n,1); +H2 = sparse(n,n); + +I = abs(x-lb)<=null & abs(Fx)<=null; +beta_l(I) = 1; +z(I) = 1; + +I = abs(ub-x)<=null & abs(Fx)<=null; +beta_u(I) = 1; +z(I) = 1; + +I = x-lb>=-null & Fx>=-null; +alpha_l(I) = 1; + +I = ub-x>=-null & Fx<=null; +alpha_u(I) = 1; + +Da = zeros(n,1); +Db = zeros(n,1); + +I = 1:n; + +I0 = Indexset==0; +Da(I0) = 0; +Db(I0) = -1; +H2(I0,:) = -DFx(I0,:); + +I1 = Indexset==1; +denom1 = zeros(n,1); +denom2 = zeros(n,1); +if any(I1) + denom1(I1) = max(null,sqrt((x(I1)-lb(I1)).^2+Fx(I1).^2)); + denom2(I1) = max(null,sqrt(z(I1).^2+(DFx(I1,:)*z).^2)); +end + +I1b = Indexset==1 & beta_l==0; +Da(I1b) = (x(I1b)-lb(I1b))./denom1(I1b)-1; +Db(I1b) = Fx(I1b)./denom1(I1b)-1; +I1b = Indexset==1 & beta_l~=0; +if any(I1b) + Da(I1b) = z(I1b)./denom2(I1b)-1; + Db(I1b) = (DFx(I1b,:)*z)./denom2(I1b)-1; +end + +I1a = I(Indexset==1 & alpha_l==1); +if any(I1a) + H2(I1a,:) = repmat(x(I1a)-lb(I1a),1,n).*DFx(I1a,:)+sparse(1:length(I1a),I1a, ... + Fx(I1a),length(I1a),n,length(I1a)); +end + +I2 = Indexset==2; +denom1 = zeros(n,1); +denom2 = zeros(n,1); +if any(I2) + denom1(I2) = max(null,sqrt((ub(I2)-x(I2)).^2+Fx(I2).^2)); + denom2(I2) = max(null,sqrt(z(I2).^2+(DFx(I2,:)*z).^2)); +end + +I2b = Indexset==2 & beta_u==0; +Da(I2b) = (ub(I2b)-x(I2b))./denom1(I2b)-1; +Db(I2b) = -Fx(I2b)./denom1(I2b)-1; +I2b = Indexset==2 & beta_u~=0; +if any(I2b) + Da(I2b) = -z(I2b)./denom2(I2b)-1; + Db(I2b) = -(DFx(I2b,:)*z)./denom2(I2b)-1; +end + +I2a = I(Indexset==2 & alpha_u==1); +if any(I2a) + H2(I2a,:) = repmat(x(I2a)-ub(I2a),1,n).*DFx(I2a,:)+sparse(1:length(I2a),I2a, ... + Fx(I2a),length(I2a),n,length(I2a)); +end + +I3 = Indexset==3; +phi = zeros(n,1); +ai = zeros(n,1); +bi = zeros(n,1); +ci = zeros(n,1); +di = zeros(n,1); +denom1 = zeros(n,1); +denom2 = zeros(n,1); +denom3 = zeros(n,1); +denom4 = zeros(n,1); +if any(I3) + phi(I3) = -ub(I3)+x(I3)+Fx(I3)+sqrt((ub(I3)-x(I3)).^2+Fx(I3).^2); + denom1(I3) = max(null,sqrt((x(I3)-lb(I3)).^2+phi(I3).^2)); + denom2(I3) = max(null,sqrt(z(I3).^2+(DFx(I3,:)*z).^2)); + denom3(I3) = max(null,sqrt((ub(I3)-x(I3)).^2+Fx(I3).^2)); + denom4(I3) = max(null,sqrt(z(I3).^2)); +end + +I3bu = Indexset==3 & beta_u==0; +ci(I3bu) = (x(I3bu)-ub(I3bu))./denom3(I3bu)+1; +di(I3bu) = Fx(I3bu)./denom3(I3bu)+1; +I3bu = Indexset==3 & beta_u~=0; +if any(I3bu) + ci(I3bu) = 1+z(I3bu)./denom2(I3bu); + di(I3bu) = 1+(DFx(I3bu,:)*z)./denom2(I3bu); +end + +I3bl = Indexset==3 & beta_l==0; +ai(I3bl) = (x(I3bl)-lb(I3bl))./denom1(I3bl)-1; +bi(I3bl) = phi(I3bl)./denom1(I3bl)-1; +I3bl = Indexset==3 & beta_l~=0; +if any(I3bl) + ai(I3bl) = z(I3bl)./denom4(I3bl)-1; + bi(I3bl) = (ci(I3bl).*z(I3bl)+(di(I3bl,ones(1,n)).*DFx(I3bl,:))*z)./denom4(I3bl)-1; +end + +Da(I3) = ai(I3)+bi(I3).*ci(I3); +Db(I3) = bi(I3).*di(I3); + +I3a = I(Indexset==3 & alpha_l==1 & alpha_u==1); +if any(I3a) + H2(I3a,:) = repmat(-lb(I3a)-ub(I3a)+2*x(I3a),1,n).*DFx(I3a,:)+... + 2*sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a)); +end +I3a = I(Indexset==3 & alpha_l==1 & alpha_u~=1); +if any(I3a) + H2(I3a,:) = repmat(x(I3a)-lb(I3a),1,n).*DFx(I3a,:)+... + sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a)); +end +I3a = I(Indexset==3 & alpha_l~=1 & alpha_u==1); +if any(I3a) + H2(I3a,:) = repmat(x(I3a)-ub(I3a),1,n).*DFx(I3a,:)+... + sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a)); +end + +%H1 = sparse(1:n,1:n,Da,n,n,n)+Db(:,ones(n,1)).*DFx; +H1 = bsxfun(@times,Db,DFx); +H1 = spdiags(diag(H1)+Da,0,H1); +H = [lambda1*H1; lambda2*H2];