AAfun.m0000644001023000102300000000011607045615720013041 0ustar hwolkowihwolkowicz%% function function [y] = dnRk(x,R,k); y = (delsq(numgrid(R,k))) * x; Afun.m0000644001023000102300000000011307040701057012727 0ustar hwolkowihwolkowiczfunction [v] = Afun(u) global A counter v = A*u; counter = counter + 1; Dfun.m0000644001023000102300000000011507041077545012745 0ustar hwolkowihwolkowiczfunction [v] = Dfun(u) global A D counter v = D*u; counter = counter + 1; arruntrs.m0000644001023000102300000000024607046615627013742 0ustar hwolkowihwolkowiczglobal A D counter %par=struct('matA','Afun','linear',b,'radius',r,'gaptol','1e-9'); par.gaptol=1e-9; 'in run' keyboard [xopt,fopt,lopt,par] = artrust(matA,b,r,par); artrust.m0000644001023000102300000010106007046616450013555 0ustar hwolkowihwolkowiczfunction [xopt,fopt,lopt,par] = ... artrust(matA,a,s,par) %
``` DATE=Feb 4, 2000
%NEWTRUST Finds an approximate solution to the trust region subproblem.
%
%  ARTRUST solves: min      f(x) = x'Ax - 2a'x
%                 such that        s^2-x'x >=0    (Lagr. mult. lopt <=0)
%%  which is the trust region subproblem. This can be used within
%%     a trust region algorithm for min f(x).
%
% The algorithm returns an approximation xopt of the true optimum x*;
%   and a corresponding value fopt=f(xopt) so that:
%         norm(xopt) <= (1+normtol)s
%   and two constants mulow,muup so that
%   mulow <= f(x*) <= muup,   mulow <= fopt <= muup,  AND
%  |f(x*) - fopt| <=  |muup-mulow| <= |2*dgaptol*muup| <= |2*dgaptol*fopt|
%  So that f(x*) >= (1-2*dgaptol)*fopt
%
%  [xopt,fopt] = artrust(A,a,s) returns the approximations xopt and
%    fopt for the data A,a,s
%
%  [xopt,fopt,lopt,par] = artrust(A,a,s,par) returns the approximations
%    fopt and xopt as well as an approximation lopt to the optimal
%    Lagrange multiplier. The structure par contains parameters on
%    input and output (described below).
%
% %%%%%%%%%%%%%%%%%%%%%% INPUT VARIABLES:
% A   The first input argument is either a real square
%     symmetric matrix (which can be full or sparse) or a
%     string containing the name of an M-file which applies a linear
%     operator to the columns of a given matrix.
%     For example, newtrust('fft',...) is much faster than EIGS(F,...)
%     where F is the explicit FFT matrix. This function is passed
%     to a sparse eigenvalue routine.
%
% a   is an n-vector which is the linear term of the quadratic.
% s   is the radius of the trust region.
%
% par is a structure with several optional inputs.
%     par.zerotol
%     par.mulow, par.muup
%%    ??????
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% optional files needed: Afun.m,Dfun.m, and this file artrust.m  ???
%%%        option file conj_grad.m, see commented out part below  ???
%%
% BEFORE USING NEWTRUST, THE FOLLOWING INSTRUCTION MUST BE ENTERED
%           > global A D counter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Trust Region Subroutine Algorithm  Trust Region Subroutine Algorithm  Algorithm with Links to Documentation  Charles Fortin and Henry Wolkowicz  Dept. Comb. & Opt., Univ. of Waterloo
%This file is available with URL http://orion.math.uwaterloo.ca/~hwolkowi/henry/software/trustreg.d/newtrust.html
% In Documentation:  Complete Documentation File, Variable Definitions,  In (this file) Algorithm:  Tolerances,  Initialization, While Loop, End of While Loop,
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%
%
global A D counter

if ~ischar(matA)
[n,n] = size(A);
if any(size(A) ~= n)
error('Matrix A must be a square matrix or a string.')
return
end
end
n=length(a);
%%%%%%%
%if s < zerotol,   %????too small???? what is too small here???
%   xopt=0;
%   fopt=0;
%   lopt=0;
%end                  %?????????????
%%%%%%%
if n==1,
if A > 0 & abs(xopt) <= s,
xopt=a/A;
lopt=0;
else
xopt=sign(a)*s;
lopt=A+b/xopt;
end
fopt= A*xopt^2-2*a*xopt;
return
end
%%%%%%%%%%

%
%
%%%%%%%%%%%%%%%%%%%%%% TOLERANCES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tolerances - these can and should be changed by the user.
%
if nargin < 4,
par.gaptol=1e-9;   % default tolerance
dgaptol=par.gaptol;    % duality gap tolerance
else
dgaptol=par.gaptol;     % duality gap tolerance
end
intol= dgaptol;
normtol  =intol;         % allowable error in feasibility for stopping
zerotol=intol/log(n^100); % tolerance for zero for y(1) for hard case
% since y is norm 1, tolerance decreases with dim.
kzerotol=zerotol;       % tolerance for zero for k(t)
hzerotol=1e-4;      % hard case detection  ???
iterbnd=30;             % upper bound on iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%END TOLERANCES%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%INITIALIZATION%%%%%%%%%%%%%%
t0=cputime;             % saved to get total cputime usage
counter = 0;            % counts matrix mults
a=full(a);              %  set to full to avoid conflicts
na=norm(a);             % vector norm
if na    <   zerotol,
disp('WARNING: linear part -a- is approximately 0');
end

%%%First eigenvalue calculation
OPTIONS.tol=1e-13;   % convergence tolerance
OPTIONS.issym = 1; % because we know the matrices A and D are symmetric
OPTIONS.isreal = 1; % A is real
OPTIONS.disp = 0; % no display of the output for eigs.m
v=full(sum(A))';v=-v/norm(v);  % good/cheap normalized estimate ???
%i.e. largest eig of -A using v=-Ae, e=ones
OPTIONS.v0 =v ; % initial eigenvector estimate for sparse A
%%%a and t are passed to Afun ????
tarn10=cputime; % keep for the first eig calc time
if ischar(matA)
[va,lambda,flageigs]=areigs(matA,n,1,'SA',OPTIONS); % computes smallest
else
[va,lambda,flageigs]=areigs(A,1,'SA',OPTIONS); % computes smallest
end
tarn1f=cputime-tarn10;   % save initial eig calc cputime for statistics
disp(['elapsed time for first eig calc:  ',num2str(tarn1f)]);
if (flageigs>0),
disp(['condition est. of eigs:  ',num2str(condest(A))]);
end
%           %eigenvalue (lambda) and corresponding eigenvector (va) for A
avnorm=norm(A*va-lambda*va);
%if avnorm > zerotol*(abs(lambda)+1),   %  ???zerotol???check???
%      disp(['USING eig, output from eigs norm error ',num2str(avnorm)]);
%      [V,Lam]=eig(full(A));
%      lambda=min(diag(Lam));
%      indx=find(lambda==diag(Lam) );va=V(:,indx);va=va/norm(va);
%      avnorm=norm(A*va-lambda*va);
%      disp(['USING eig new norm error is  ',num2str(avnorm)]);
%end

lambdamin = lambda; % save the min eigenvalue of A
vaa=va'*a;
% ensure va'*a is not negative
if (vaa < 0)
va = -va;vaa=-vaa;
end
y=[1;va];y=y/norm(y);  % initial normalized estimate eigenv. D(t^*)
flag = 0; %flag will change value if the almost hard case (case 2) occurs

alow=lambda;aup=lambda;
% lower estimate on mu*
if lambda<0
mulow = alow*s^2-2*na*s;  % lower estimate for optimal value
else
% mulow = -2*na*s;  %  replaced by min on boundary or interior
% since interior opt value is a'invAa-2a'invAa=-a'invAa
mulow = min(alow*s^2-2*na*s,-na*s);
end
% upper estimate on mu*
if lambda > 0
if (( (vaa)/(lambda )  <=  s)),   %  optimum may be on boundary or interior
testmu = s^2*lambda - 2*s*na;
muup  = -(vaa)^2/lambda;
if muup < testmu,   %  apply conj grad - opt interior
a=sparse(a);
% apply conjugate gradient method to  find the optimum.
xopt=A\a;     % use matlab function for now or uncomment conj_grad

fopt  = -a'*xopt;   % optimal value since A xopt = a
lopt = 0;
topt = 99999;
muup=99999; mulow=99999;iter=0;
gap=[topt topt
lopt lopt
topt topt
];
t0f=cputime-t0;
flag = 2; % To know answer is the unconstrained minimum
disp(['elapsed time: ', num2str(t0f)]);
disp(['Problem solved using conjugate gradients  ']);
relopt=-norm(  (A-lopt*eye(size(A,1)))*xopt-a)/fopt;
%rel opt cond. norm grad of Lagrangian / value of Lagrang. at opt.
format long e
disp(['relative opt. condition:  ', num2str(relopt)])
return;    %OPTIMALITY - PROBLEM IS SOLVED
end
bestxf = ((vaa)/lambda)*va; % vaa is >=0 now
else
muup = s^2*lambda - 2*s*(vaa);  % <=0 by above if
bestxf = s*va;
end
else
[muup,indice] = min([0,lambda*s^2-2*s*(vaa)]);
if indice == 1
bestxf = zeros(size(a,1),1);
else
bestxf = s*va;
end
end
bestfxf=muup; % current best opt value corresp to bestxf
% bestxf is feasible (possibly interior)
mubnds=[mulow muup];    % keep history of bounds - latest at top
% This is only for debugging help.
gap=[mulow muup];       % duality gap - gets updated when a bound changes

if (lambda > 0)
low = min(alow - na/s,-na*s);   %lower bound on t
else
low = alow - na/s;
end
high = aup + na*s;      % upper bound on t
newlow=low;newhigh=high;   %new points found by vertical cuts
vertflag=0;             % to prevent 2 vertical cuts in a row

ss1=s^2+1;
D = [ 0 -a'
-a   A];           % augmented matrix
ystar = 1 / sqrt( s*s +1);      % this is value we are after for y0(t)
iter=0;         % iteration count
b6 = '      ';              % define 6 blanks

% ptgood, ptbad variables: used for saving interpolation data on the
%  good (easy) and bad (hard) side of maximum of k(t).
%  gy,by saves last y vector on good and bad side
% The latest data is put into row 1.
%   Only two rows are needed. The rest is saved for debugging help.
%    t, y(1), lambda,  valid-data-indicator ; is data saved.
ptgood=[];
gy=[];by=[];    % The latest vectors y from the good and bad side resp.
side=0;        % indicator for consecutive t found
%  -i for i bad sides and +i for i  good sides
tslopes=[];    % Save slopes of kt on the [bad good] sides - for
%              debugging and for case where kt is almost piecewise linear.
kts=[];      % saving kt bounds and t value - for debugging
t=99999;indexg=0;indexb=0;  % initialize in case of 0 iterations
tind=[];     % follow trail of t changes - for debugging
%
disp('iter ----- t  ----  rel dual gap -- rel. primal opt cond');
%txt=[int2str(iter),b6,num2str(t),b6,num2str(dgap),b6,num2str(dbestgap)];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%NOTES on ALGORITHM%%%%%%%%%%%%%%
%   Following are the functions of interest:
%    h(lambda)  =  s^2lambda -a'(A-lambda I)^{-1} a;      maximize - concave
%    hp(lambda)  = s^2 -a'(A-lambda I)^{-2} a;solve = 0 - concave, nonincreasing
%    phi(lambda)  = 1/s-1/||(A-lambda I)^{-1} a||;solve = 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    k(t)  =  (s^2+1)lambda -t         maximize - concave
%    k'(t) = (s^2+1)y0^2(t) -1         solve = 0  - same sign as f(t)
%    psi(t)  = -1/y0(t) + sqrt(s^2+1)    solve = 0 - concave, nonincreasing
%    y0(t) = first component of y normalized   - nonincreasing, 1to0
%    lambda=lambda(t)    concave function
%
% The algorithm is based on maximizing k(t).
% We interpolate on psi(t) whenever possible.
%
% The general iteration interpolates psi(t) whenever
%   possible. Duality is used to calculate bounds, i.e. the function k(t) is
%   used when problems with the hard case occur.

% Initialize the stopping criteria
t1=abs( 1-y(1)^2*(1+s^2))>( (2*normtol+normtol^2) *y(1)^2*s^2 );
dbestgap= (bestfxf-mulow)/(abs(bestfxf)+1);
t2= dbestgap   > 2*dgaptol;
dgap=  abs(muup-mulow)/(abs(muup)+1) ;
t3=  (muup-mulow) >= -2*dgaptol*muup;
%t3= dgap >  dgaptol; % Duality gap criteria    - old way
t4=iter <=iterbnd; % iteration bound
t5=abs(high-low)/(abs(high)+1) > zerotol; %interval of uncertainty for t
%%%%%%%%%%%%%%%%%%%END INITIALIZATION%%%%%%%%%%%%%
%%
%
%%%%%%%%%%%%%%%%%%%%%%%WHILE LOOP%%%%%%%%%%%
%      PART 1. Find New Point t,
%              Vertical Cut,
%              Triangle Interpolation,
%              Good Side Interpol.,
%      PART 2. Updating Begins.,
%             Lambda has wrong sign,
%             Good Side Update,
%       PART 3.Tolerance Updates, End of While Loop, Start of Output,
%
%% check feasibility and iteration count and duality gap
while ( (t1&t2) | t3) & (t4 & t5),  % start of while%
iter=iter+1; % one more iteration is done
indexg=size(ptgood,1); % number of points on the good side
%
%%%%%%%%%  FIND a new point t.
newlow=max(newlow,low);newhigh=min(newhigh,high);
if newlow > low,
t=newlow;
tind=[tind -2];    % indicator for t value - for debugging
elseif newhigh < high,
t=newhigh;
tind=[tind -1];
else
t=(max(newlow,low)+min(newhigh,high))/2;      % midpoint default
tind=[tind 1];
end

if indexb > 0 & indexg > 0,         %
%  The last point on right and left that has eigenvalue data
% always do triangle interpol and vertical cut if possible
kth=ss1*ptgood(1,3)-ptgood(1,1);
thslope=ss1*ptgood(1,2)^2-1;

%
%
%%%%%%%%%%%%%%%%%%%%%%%TRIANGLE INTERPOLATION%%%%%%%%%%%

% Always find interpolant using triangle-with associated upper bound.
tslopes=[tlslope thslope
tslopes];         % for debugging purposes - possible use
% when both sides are almost linear
(thslope*(-ptgood(1,1))+kth)     ];
% rhs of interpolation
ktt=[ 1  (-tlslope)
1  (-thslope)];
kttans=ktt\mrhs; % find intersection of the two tangent lines at
% the last bad and good side points. The point of
% intersection is (kttans(2),kttans(1)). Hence, kttans(1)
% is an upper bound on mu^* and kttans(2) is an approximation
% for t^*
ktup=kttans(1);
%%%%%%%?????testing the next if ????!!!
%%%???change!!!!in case near hard case - this is always too big
%%%???a step away from good side?????
if ptgood(1,2) >zerotol,   % avoided hard case here - try
tt=kttans(2);   % This is alternate interpol point using triangle
else
tt=.9*newhigh+.1*tt;
end

if (ktup <=muup) & (ktup >= mulow),
muup=ktup;           % new upper bound on optimal value
end
if max(newlow,low)<=tt & tt <=min(newhigh,high);
t=tt;
tind=[tind 2];
end
%  %%%%%%%%%%%%%%%%%%%%End of triangle for new interpolant

%
%
%%%%%%%%%%%%%%%%%%%%%%%VERTICAL CUT%%%%%%%%%%%
%
if ptgood(1,2) >zerotol,   % avoided hard case here - try
%%% smaller values for (was) hzerotol!!!
if vertflag==1,                 %
vertflag=0;                 %
else                 %
vertflag=1;          % prevents 2 vertical cuts in a row
%% errors can occur if eigenvalue not calculated properly
%% or in hard case - here and below!!!
if kth > (1+kzerotol)*ktl;  % opt cond for k(t) - check???
tempnewlow=low + (kth-ktl)/tlslope;
if newlow <= tempnewlow & tempnewlow <= newhigh,       %
newlow=tempnewlow;
t=newlow;
tind=[tind -2];
else
disp('WARNING newlow not in interval ?????')
%%????this is not explained satisf. why??? roundoff
%%%that important???
newlow=low;            %
end            %
elseif ktl > (1+kzerotol)*kth;  %opt cond for k(t)-check?
tempnewhigh=high + (ktl-kth)/thslope;
if newlow <= tempnewhigh & tempnewhigh <= newhigh,  %
newhigh=tempnewhigh;
t=newhigh;           %
tind=[tind -1];
else           %
disp('WARNING newhigh not in interval ?????')
end                           %
else                             %
disp('two values of k up to tolerance found');
% arrange proper output still????RETRN!!!
end                 %
end                  %end vertical flag check
end          %  end of hard case check above vert flag -  temp
%%%%%%%%%%%%%%%%%%%%%%%end vertical cut%%%%%%%%%%%

end     % end of if indexb > 0 and indexg > 0

%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% INVERSE INTERPOLATION - good side
%????how can these two lines not be executed???even though the if is executed???
if side>0,   % last iterate was on the GOOD SIDE

if  sum(ptgood(:,4)) >=2, %two t-iterates with data good side
tint=ptgood(1:2,1);
yint=ptgood(1:2,2);
f1=-1/yint(1)   + 1/ystar;
f2=-1/yint(2)   + 1/ystar;
slope=(tint(1)-tint(2))/(f1-f2);
tt=tint(1) - f1* slope;    % linear inverse interpolation
disp('before NaN')
keyboard
if newlow <=tt&tt<=newhigh,     % roundoff error problem check
%% We now decide whether to chose the new point
t=tt;
tind=[tind 3];     % indicator for t value - for debugging
end
if side>4&(min(ptgood(1:2,2)) < (10^(side-1)*zerotol) ), %
%too many iterates on good side and it is near bad case
t=.95*newlow+.05*newhigh;
%????use weights from k' if indexb>0?????
tind=[tind 4];     % indicator for t value - for debugging
end               %
elseif sum(ptgood(:,4)) >=1 & indexb>0 >= 1,
if ptgood(1,2)>zerotol & tind(length(tind)) >= 0,
% not near hard case  and not vertical cut
% interpolate with one point from each side if not hard case
tintg=ptgood(1,1);
psig=-1/ptgood(1,2)+sqrt(ss1);
tt=(tintg*psib-tintb*psig)/(psib-psig);
if newlow <= tt & tt <= newhigh,
t=tt;
tind=[tind 31];
end
else
%use spline on k(t) - to do still
tind=[tind 32];
end
end   % end of if sum(ptgood) >=2 else >=1

%%%%%%%%%%%%%%%%%%%%%%%%%
%
elseif side<0,   % last iterate was on the BAD SIDE
%INVERSE INTERPOLATION
% t(psi) inverse function of psi(t), then we linearly interpolate
% value t(0), which is t^*, using two points from the bad side.
f1=-1/yint(1)   + 1/ystar;
f2=-1/yint(2)   + 1/ystar;
% If f1=f2, then t value repeated by e.g. t=low in triangle.
% This indicates that the optimal t is very close to t=low.
if f2>f1,   % correct monotonicity occurred
slope=(tint(1)-tint(2))/(f1-f2);
tt=tint(1) - f1* slope;    % linear inverse interpolation
if newlow<=tt&tt<(newhigh-zerotol),
% Use the new point - it is not too large
% This should go to good side
t=tt;
tind=[tind 5];     % indicator for t value - for debugging
elseif side<-4,   % need a point on good side desperately
% was following two lines - need to speed this up using a
%   CONCAVE spline  for k(t) in hard case!!!!
t=.95*high + .05*low;
tind=[tind 7];     % indicator for t value - for debugging
end
end                   %
elseif  indexg>0 & ptgood(1,2)>hzerotol & tind(length(tind)) >= 0,
% not near hard case and not vertical cut
% this interpolation makes no sense if near hard case
% interpolate with one point from each side if not hard case
tintg=ptgood(1,1);
psig=-1/ptgood(1,2)+sqrt(ss1);
tt=(tintg*psib-tintb*psig)/(psib-psig);
if newlow <= tt & tt <= newhigh,
t=tt;
tind=[tind 31];
if side<-3,   % need a point on good side desperately
% was following two lines - need to speed this up using a
%   CONCAVE spline  for k(t) in hard case!!!!
t=.95*high + .05*low;
tind=[tind 7];     % indicator for t value - for debugging
end
else
%In case near hard case above !!!!
%  %use spline on k(t) - to do still
tind=[tind 32];
end
end    %  end for sum(ptbad) >=2 and >=1
end    % end for side check
disp(['side: ', num2str(side)]);
disp([sprintf('trail:'), sprintf('%d,%d,',tind) ])

%

%  A new point t has been found. Start updates.
%
% Recall t is the most recent update for t^*

D( 1, 1) = t;

tarn20=cputime;

OPTIONS.tol=1e-13;   % convergence tolerance
OPTIONS.v0 = y;       %  default
if tind(length(tind))==3,  % finished interpol on good side
OPTIONS.v0 = gy;
elseif tind(length(tind))==31,  % finished interpol using both sides
OPTIONS.v0 = by;
elseif tind(length(tind))==7 & indexg > 0, %  very close to high
OPTIONS.v0 = gy;
elseif tind(length(tind))==-2, %  close to low
OPTIONS.v0 = by;
elseif tind(length(tind))==-1, %  close to high
OPTIONS.v0 = gy;
elseif tind(length(tind))==4 & indexb > 0, %  very close to newlow
OPTIONS.v0 = by;
elseif tind(length(tind))==5 & indexg > 0, %  go to good side
OPTIONS.v0 = gy;
elseif (t-low) < ( (high-t) - 100*zerotol ) & indexb > 0,
OPTIONS.v0 = by;
elseif (t-low) > (  (high-t) + 100*zerotol )& indexg > 0,
OPTIONS.v0 = gy;
end
if ischar(matA)
[va,lambda,flageigs]=areigs('Dfun',n+1,1,'SA',OPTIONS); % computes smallest
else
[va,lambda,flageigs]=areigs(D,1,'SA',OPTIONS); % computes smallest
end
%[y,lambda,flageigs] = eigs('Dfun',n+1,1,'SR',OPTIONS);
tarn2f=cputime-tarn20;  % save initial Lanczos call cputime for statistics
disp(['tind(last) and elapsed time for second eig calc:  ', ...
num2str(tind(length(tind))),' ,  ',num2str(tarn2f)]);
if (flageigs>0),
disp('flag WARNING from eigs EIGENVALUE CALCULATION BAD');
disp(['condition est. of eigs:  ',num2str(condest(D))]);
end
%                   % y and lambda are the new estimates
%                   % for the eigenvector of D(t^*) and for lambda^*
%dynorm=norm(D*y-lambda*y);
%if dynorm > zerotol*(abs(lambda)+1),   %  ???zerotol???check???
%      disp(['USING eig, output from eigs norm error ',num2str(dynorm)]);
%      [V,Lam]=eig(full(D));
%      lambda=min(diag(Lam));
%      indx=find(lambda==diag(Lam) );y=V(:,indx);y=y/norm(y);
%      dynorm=norm(D*y-lambda*y);
%      disp(['USING eig new norm error is  ',num2str(dynorm)]);
%end

if y(1)<0,  % ensure the sign of the first component is positive
y=-y;
end

%
%  multiplier sign check - to update lower bounds
if lambda > 0,   % A is known to be pos def - multiplier has wrong sign

if y(1) >= ystar,
% means y(2:n+1)/y(1) is the optimum solution for all points outside
% trust region of radius (1-y(1)^2)/y(1)^2 or on its boundary. If this
% radius is less than s^2, as is case here since this is equivalent to
% y(1) >= ystar, then the optimum we are looking for is inside
% the trust region and A is positive semidefinite.
% We will then apply a conjugate
% gradient method to find the optimum.
% Hence, y(2:n+1)/y(1) is feasible so optimal with lambda=0.

a=sparse(a);
%     % apply conjugate gradient method to  find the optimum.
xopt=A\a;     % use the matlab function for now or uncomment conj_grad
gap  = -a'*xopt;   % optimal value since A xopt = a
fopt=gap;
lopt = 0;
topt = 99999;
gap=[topt topt
lopt lopt
topt topt
];
t0f=cputime-t0;
flag = 2; % To know answer is the unconstrained minimum
disp(['elapsed time: ', num2str(t0f)]);
disp(['Problem solved using conjugate gradients  ']);
relopt=-norm(  (A-lopt*eye(size(A,1)))*xopt-a)/fopt;
%rel opt cond. norm grad of Lagrangian / value of Lagrang. at opt.
format long e
disp(['relative opt. condition:  ', num2str(relopt)])
return;    %OPTIMALITY - PROBLEM IS SOLVED
end
else     % the sign of the multiplier is correct
%   so apply duality theory
kt=ss1*lambda-t;         % value of k(t), dual functional. This gives
% a lower bound on mu^*
kts=[t lambda kt
kts];
if (kt>mulow)&(kt<=muup),
mulow=max(mulow,kt);           % new lower bound on optimal value
end
%
end        % end of if lambda > 0 loop
%%  end multiplier sign check

%
%
if y(1) >= ystar,        % bad side
%%%%???? when y(1) very small and ystar very small, then error can
%%% occur about knowing which side - roundoff problem can occur
%%%   here???? analysis needed.
if t==newhigh;      %% roundoff error problem here ????why???
%%%newhigh was on wrong side
newhigh=high;
end
low = t;           % new lower bound on t^*
ptbad];       % Keep points for first interpolation.
by=y;
%%
%%
%%%%%%%%%%%%%
%% Begin bound for PSDP from bad side - approx. negative curvature
%  direction and move to the boundary.
if indexg>0,
if  min(ptgood(:,3)) <= 0,   % multiplier has correct sign on good side
normx2=(1-y(1)^2)/(y(1)^2); % square of norm of current soln.
xf=y(2:(n+1))/y(1);   % current soln. inside trust region

% In the hard case (case 2), we use MoSo formula page 558 - however
% with negative eigenvalue and eigenvector, to move to take a step to
% the boundary. In the easy case and hard case (case 1), we take a step
% to the boundary by taking the solution between the line joining the
% xf and the x-point - gy(2:n+1)/gy(1,1) -  outside the trust region
% associate with the last good side t-point.

if (gy(1) < zerotol)
% This is the almost hard case (case 2)
z = gy(2:n+1);
normz=norm(z); % normz should be almost 1

else
% This is the easy case or almost hard case (case 1)
z = gy(2:n+1)/gy(1,1)-by(2:n+1)/by(1,1);
normz=norm(z);

end

xyprod=(z'*xf)/normz;  % normalized inner product
ndiff=s^2-normx2;
sxy=sign(xyprod);
if sxy==0,
sxy=1;
end
%%%??????check whether this is really the correct direction that
%%% keeps the optimality conditions - seems perhaps not???
%% or perhaps because of roundoff??? z should be in null space
%%% of gradient of Lagrangian ????? or because it is on good side???
tau=ndiff/(xyprod+sxy*sqrt(xyprod^2+ndiff) );
xnew=xf+(tau/normz)*z;  % xnew is the point obtained by moving to the
% boundary
fxnew=xnew'*A*xnew-2*xnew'*a; %should improve at point xf
% get new upper bound using the feasible xnew - need to allow almost
%  best xnew as the bestxf
if fxnew=mulow),
muup=fxnew;
end

end      % end if <= 0 line mult. check
end        %end of if indexg line
%%%%%%%%%% end of PSDP bound
%
%
%%%%%%%%%%%%%Good side update
%%%%%%%%%%%%%%%%%%%
else                    %  current point is on good side
if t==newlow;      % roundoff error problem correction - but
newlow=low; why did it occur????
end
high = t;       % we have new upper bound on t
ptgood=[t y(1) lambda 1
ptgood   ] ;
side=max(1,side+1);  % consecutive good sides
gy=y;
if  ptgood(1,3) <= 0,   % current multiplier has correct sign
if indexb > 0,
%% begin type 1 bound for PSDP from the good side
normx2=(1-by(1)^2)/(by(1)^2);  % use last bad side point
xf=by(2:(n+1))/by(1);

% In the hard case (case 2), we use MoSo formula page 558 - however
% with negative eigenvalue and eigenvector, to move to take a step to
% the boundary. In the easy case and hard case (case 1), we take a step
% to the boundary by taking the solution between the line joining the
% xf and the x-point - gy(2:n+1)/gy(1,1) -  outside the trust region
% associate with the last good side t-point.

if (gy(1) < zerotol)
% This is the almost hard case (case 2)
z = gy(2:n+1);
normz=norm(z); %normz
else
% This is the easy case or the almost hard case (case1)
z = gy(2:n+1)/gy(1,1)-by(2:n+1)/by(1,1);
normz=norm(z);
end

xyprod=(z'*xf)/normz;  % normalized inner product
ndiff=s^2-normx2;
sxy=sign(xyprod);
if sxy==0,
sxy=1;
end
tau=ndiff/(xyprod+sxy*sqrt(xyprod^2+ndiff) );
xnew=xf+(tau/normz)*z;  %
fxnew=xnew'*A*xnew-2*xnew'*a;
% get new upper bound using the feasible xnew - need to allow almost
%  best xnew as the bestxf
if fxnew=mulow),
muup=fxnew;
end

else      % else of indexb>0 line
% Need upper bound as well in case triangle was not used.
%
%% begin type 2 bound for PSDP
xf=y(2:(n+1));    % norm is sqrt of 1-y_0^2
xnew=(s/sqrt(1-y(1)^2))*xf;
% xnew is the point obtained by moving to the  boundary
fxnew=xnew'*A*xnew-2*xnew'*a;
% get new upper bound using the feasible xnew
if (fxnew=mulow),
bestxf=xnew;
bestfxf=fxnew;
muup=fxnew;
end
%% reached end type 2 bound for PSDP
end     % ind of if indexb > 0 line
end      % end if <= 0 line - multiplier sign check
end        % end of if bad side
%%%%%%%%%%%%
%% save new bounds and gap
if mulow > mubnds(1,1) | muup < mubnds(1,2),
mubnds= [ mulow muup
mubnds ];
gap=[mulow muup];    % duality gap - gets updated when a bound changes
end

%
% update tolerances
t1=abs( 1-y(1)^2*(1+s^2))>( (2*normtol+normtol^2) *y(1)^2*s^2 );
% use y(1)=ystar at optimum
%t1=abs( 1-y(1)^2*(1+s^2))>(normtol*y(1)^2*s^2 ); % old Feas criteria
dbestgap= (bestfxf-mulow)/(abs(bestfxf)+1);
t2= dbestgap   > 2*dgaptol; %??????why 2???change: rel dual opt cond
% if t1 feasibility fails, but bestfxf works, then feasibility is fine
if muup==0,
muupt=muup+1;
else
muupt=muup;
end
reldgap=  abs(muup-mulow)/(abs(muupt)) ;
t3=  (muup-mulow) >= -2*dgaptol*muup;
t4=iter <=iterbnd; % iteration bound criteria
t5=abs(high-low)/(abs(high)+1) > zerotol; % gap for t^* criteria
txt=[int2str(iter),b6,num2str(t),b6,num2str(reldgap),b6,num2str(dbestgap)];
disp( txt);
%     end of while loop

end         % end of while for gap and iteration check
%
%%%%%%%%%%%%%%%%%%%%%%%%%END GENERAL ITERATION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%OUTPUT DATA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fopt=(muup+mulow)/2;
if ~t2,           % bestxf is almost feasible
xopt=bestxf;
fopt=bestfxf;
else
xopt = y(2:n+1)/y(1);
end
topt=t;lopt=lambda;
if indexg ==0,
ptgood=[99999 99999 99999 0];
end
if indexb ==0,
end
t0f=cputime-t0;
gap=[gap                    % lower and upper bound for mu^*
ptbad(1,3) ptgood(1,3) % lower and upper bound for lambda^*
low high               % lower and upper bound for t^*
tarn1f t0f ];          % tarn1f is elapsed time for first eigs call
% t0f is total computation time

disp(['elapsed time: ', num2str(t0f)]); % t0f is total computation time
disp(['number of matrix-vector multiplications:  ', num2str(counter)])
disp(['number of iterations:  ', num2str(iter)])
relopt=-norm(  (A-lopt*eye(size(A,1)))*xopt-a)/fopt;
%rel opt cond. norm of grad of Lagrangian / value of Lagrang. at opt.
format long e
disp(['relative opt. condition:  ', num2str(relopt)])
disp(['interval of optimal value: [ ',num2str(mulow),' to ',...
num2str(muup),' ]']);

disp('final results for  ');
disp('iter ----- t  ----  rel dual gap -- rel. dual opt cond');
txt=[int2str(iter),b6,num2str(t),b6,num2str(reldgap),b6,num2str(dbestgap)];
disp( txt);

if (abs(lambdamin - lopt)/(abs(lambdamin)+1) < 1e-2)
flag = 1;
end

% Top of algorithm page. This directory on orion
%
%%%%%%%%%%%%%%%%%%%%%%%%%WARNING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Before running 'conj_grad', 'A' has to be defined and set to be
% global. We also need to define the global variable counter to take
% into account the number of matrix vector multiplications when this
% function is used with 'newtrust.m'. The following instruction need
% be done
%
% >> clear A counter
% >> global A D counter
% >> A = [.....];
% >> counter = 0;
%
% If 'newtrust' is the function we are interested in and 'conj_grad'
% might be called, then before using 'newtrust', do the following
%
% >> clear A D counter
% >> global A D counter
% >> A = [...]
% (Note: counter now doesn't need to be initialized as it will be
%  in 'newtrust')
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function solves the problem
%
% min (1/2)*x^tAx - a^tx
%
% and A is assumed positive definite
%
%
% INPUT:    a     = a vector
%           x0    = a initial estimate for x_opt
%
% OUTPUT:   x_opt = the global minimum
%
% STOPPING CRITERIA: norm(A*x-a) < zerotol:= max(10^(-15),10^(-5)*norm(a));
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global A D counter

% Initialize variables

%zerotol = 1e-2;
zerotol = max(10^(-15),10^(-5)*norm(a));
%zerotol = norm(a)*min(0.1,norm(a)^(0.1));

x = x0; % x is the current iterate

J = A*x0-a; % d is the gradient of the current iterate

nJ = norm(J); % norm of the gradient of the current iterate

J_1 = J; % d will be the gradient of the previous iterate

nJ_1 = norm(J_1); % norm of the gradient of previous iterate

d_1 = 0; % previous step direction

counter = counter +2;  % counter counts the number of matrix vector
% multiplications

while (nJ > zerotol)

d = J + (nJ/nJ_1)^2*d_1; % current step direction

d_1 = d;

r = (J'*d)/(2*(d'*A*d)); % step lenght

x_1 = x; % previous iterate is x_1

J_1 = J;

nJ_1 = nJ;

x = x - r*d;

J = A*x-a;

nJ = norm(J);

counter = counter + 2;

end

x_opt = x;

dnRK.m0000644001023000102300000000011607045615435012710 0ustar  hwolkowihwolkowicz%% function
function [y] = dnRk(x,R,k);
y = (delsq(numgrid(R,k))) * x;
gould2.m0000644001023000102300000001310207045454633013246 0ustar  hwolkowihwolkowiczfunction [x_min,f_min,fcteval,compteur,iteration,temps,comptons,drapeau] = trmgould(x0,R,prob)

% Call: [x_min,f_min,fcteval,compteur,iteration,temps,comptons,drapeau] = trmgould(x0,R,prob)
% This is the trust region method found in Gould et al. 
% This algorithm is used to minimize an unconstrained function. The function's
% objective should be available using the file 'objective.m' and the function's
% function uses the subroutine 'newtrust.m' to solve the trust region
% subproblem
%
% INPUT VARIABLES:
%
% x0 ... initial estimate of a minimizer of the function
% R  ... initial radius of the trust region
% prob ... number of the problem to be solve (used to access the files
%
% OUTPUT VARIABLES:
%
% x_min   ... an approximate of a minimizer of the function
% f_min   ... objective at x_min. This approximates the minimum of the function
% fcteval ... number of fonctions evaluations
% compteur .. number of matrix-vector multiplications
% iteration . number of iterations done by the trust region method
% temps   ... time taken to find the approximate solution
% comptons .. relative distance of the optimal Lagrange multipliers of the
%             trust region subproblems to the minimum eigenvalue of the
%             Hessians
% drapeau ... vector of integers (0,1,2) used to distinguish the different
%             cases for the trust region subproblemes

% GLOBAL VARIABLES
global A D counter % global variables used by the newtrust4b.m subroutine

global A counter x

b = ' '; % define blanks
n = size(x0,1);

temps = cputime;

% Constant used to update the trust region radius (i.e. evaluate the
% performance of the trust region step
n1 = 0.01;
n2 = 0.95;
gam1 = 0.5;
gam2 = 2;
s = R;

x = x0; % x is the current approximation to a minimizer of the function

[fopt,g,A] = objgradhess(x,prob);; % fopt is the objective evaluated at x (i.e.
% the current approximation to the minimum
% of the function, g is the gradient and A
% is the Hessian.

% initiaze variables
fcteval = 1;
iteration = 1;
compteur = 0;

disp(['iteration',b,b,b,b,'time elapsed so far',b,b,b,b,'f(x)',b,b,b,b,b,b,b,b,'norm of the gradient']);
disp([b,b,b,b,num2str(iteration),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(cputime-temps),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(fopt),b,b,b,b,b,b,b,b,b,b,b,num2str(norm(g)),b,b,b,b]);

if (norm(g)/(1+abs(fopt)) < 0.00001) % relative error on the norm of the

% Initial approximate is good enough
x_min = x;
f_min = fopt;
return; % END OF ALGORITHM
end

while (norm(g)/(abs(fopt)+1) > 0.00001)

epsilon = 10^(-6);
dgaptol = min(epsilon,10^(-4)*norm(g)); % relative duality gap. This is used
% so that initially we don't ask a lot
% of precision for the trust region
% subproblems, but as we get close to
% the solution (gradient is small), we
% do.
par.gaptol=dgaptol;
% SOLVE THE TRUST REGION SUBPROBLEM
[d,qopt,lopt,lambdamin,iter,drapeau(iteration)] = newtrust(-g,s,par);

% d is the step direction
% qopt is the predicted value of the objective at x+d (predicted by the
% lopt is the Lagranga multiplier of the optimal solution to the trust
% region subproblem.
% lambdamin is the minimum eigenvalue of the current matrix A (the Hessian at
% x).
% drapeau(iteration) is an integer that tells us which case occurred in the
% trust region subproblem (hard case, case 1 and 2 and easy case).

comptons(iteration) = (lambdamin-lopt)/(abs(lambdamin) + 1); % relative
% distance of the optimal Lagrange multiplier to the
% minimum eigenvalue of A (Hessian at x). This is just
% used for the purpose of analysing the iterates.;

compteur = compteur + counter; % counter is the number of matrix-vector
% multiplication done in the last trust region
% subproblem. So update compteur, the total
% number of matrix-vector multiplications

f_predicted = fopt + 1/2*qopt;
f_real = objective(x+d,prob);

r = (fopt - f_real)/(fopt - f_predicted); % r is a measure of our improvement
% r=1 is good, r=0 is bad
if (r > n1)
% improvement of the objective function is good enough so take a step
x = x+d;
fcteval = fcteval + 1;
end

if (r > n2)
% improvement of the objective function is great, increase the trust region
s = gam2*s;
end

if (r < n1)
% improvement of the objective function is not satisfactory. Reduce the
s = gam1*s;
end

disp('new values');
r
s

iteration = iteration + 1; % one more iteration has been done

% display info on the iteration
disp(['iteration',b,b,b,b,'time elapsed so far',b,b,b,b,'f(x)',b,b,b,b,b,b,b,b,'norm of the gradient']);
disp([b,b,b,b,num2str(iteration),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(cputime-temps),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(fopt),b,b,b,b,b,b,b,b,b,b,b,num2str(norm(g)),b,b,b,b]);

end

% return the approximation to the optimal solutions
f_min = fopt;
x_min = x;

temps = cputime - temps; % temps = time to solve the problem

hprob1.m0000644001023000102300000000304507046131630013240 0ustar  hwolkowihwolkowiczclear
format compact
format long e
global A D counter mults SHIFTA ;
A=sparse(A);
D=sparse(D);
dgaptol=1e-9;

N = 32 ;
n = N*N ;
e = ones(n,1) ;
d = e ;
for i = N:N:n
d(i) = 0 ;
end

A = spdiags([ -e -d 2*e ], [-N -1 0 ], n , n) ;
A = (A+A')  - 5*speye(n) ;
b = rand(n,1) ;

r   = 100 ;
tmults = 0 ;
tcounter = 0 ;
avecount = 0 ;
titers = 0 ;
aveiters = 0 ;
viters = 0 ;
itermax = 0 ;
for it = 1:20
b = rand(n,1) ;
A=sparse(A);
D=sparse(D);
[xopt,fopt,lopt,lambdamin,muup,mulow,iter,flag] = ...
newtrust(b,r,dgaptol);
relfeasxopt(it)=abs(1-norm(xopt)/r);
reloptcond(it)=-norm(  (A-lopt* eye(size(A,1)) )*xopt-b )/fopt;
absoptcond(it)=-norm(  (A-lopt* eye(size(A,1)) )*xopt-b );
mults
xopt = xopt*(r/norm(xopt)) ;
t = b - A*xopt ;
mu = t'*xopt/(r^2) ;
err(it) = norm(t-xopt*mu);
tmults = tmults + mults ;
tcounter=tcounter+counter;
titers=titers+iter;
viters(it)=iter;
if iter > itermax,
save xprb1 A r b dgaptol viters
itermax=iter;
end
end
avemults = tmults/20;
avecount=tcounter/20;
!rm prb1.diary
diary prb1.diary
disp('average of matrix multiplications')
[avecount]
disp('average of iterations')
[sum(viters)/length(viters)]
disp('individual iterations per problem')
[viters]
disp('relative feasibility for each problem')
[relfeasxopt]
disp('absol. optimality error for each problem')
[absoptcond]
disp('relative optimality error for each problem')
[reloptcond]
diary off
hprob2.m0000644001023000102300000000515707045315477013263 0ustar  hwolkowihwolkowiczclear
format compact
format long e
global A D counter mults SHIFTA ;
A=sparse(A);
D=sparse(D);
n = 1000;
%tmults = zeros(2,1) ;
tcounter = zeros(2,1) ;        % added by henry - tmults commented out
avecount = zeros(2,1) ;        % added by henry - tmults commented out
titers = zeros(2,1) ;        % added by henry - tmults commented out
aveiters = zeros(2,1) ;        % added by henry - tmults commented out
relfeasxopt=0;
reloptcond=0;
reldualoptcond=0;
ij=0;
itermax=0;
for r = [ 100 10 ]    %  added by henry for outside loop
for i = 1:20
d = rand(n,1) - .5 ;
u = rand(n,1) - .5 ;
u = u*(sqrt(2)/norm(u)) ;
b = rand(n,1) - .5 ;
b = b/norm(b) ;
y = d.*u ;
c = u'*y ;
A = diag(d) - y*u' - u*y' + (u*c)*u' ;
d = diag(A) ;
B = abs(A - spdiags(d,,n,n)) ;
SHIFTA = max(d+sum(B)') ;

OPTIONS.issym = 1; % know the matrix is symmetric
OPTIONS.disp = 0; % no display of the output for eig

[v,d] = eigs('mult',n,1,OPTIONS) ;
d(1)+SHIFTA

j = 0 ;

%   for r = [ 100 10 ]
r
%       [lopt, xopt, topt, fopt, iter, gap,flag_conj,tarn1f,t0f] = ...
%           newtrust4(b,r) ;
%       updated newtrust4b code with updated output line, added by henry
A=sparse(A);
D=sparse(D);
[xopt,fopt,lopt,lambdamin,muup,mulow,iter,flag] = ...
newtrust(b,r,dgaptol);
xAopt=  (A-lopt* eye(size(A,1))  ) \b;
relfeasxopt(ij,i)=abs(r-norm(xopt))/r;
reloptcond(ij,i)=-norm(  (A-lopt* eye(size(A,1)) )*xopt-b )/fopt;
philopt=lopt*r^2-b'*xAopt;    %dual function value
reldualoptcond(ij,i)=(r^2-norm(xAopt)^2)/philopt;
mults
xopt = xopt*(r/norm(xopt)) ;
t = b - A*xopt ;
mu = t'*xopt/(r^2) ;
err = norm(t-xopt*mu);
j = j + 1 ;
%       tmults(j) = tmults(j) + mults ;
viters(ij,i)=iter;         % total record of iterations
if iter > itermax,
save xprb2 A b r dgaptol viters
itermax=iter;
end
end
%avemults = tmults/20
end
!rm prb2.diary
diary prb2.diary
disp('average of matrix multiplications')
[avecount]
disp('individual iterations per problem')
[viters]
disp('relative feasibility for each problem')
[relfeasxopt]
disp('relative optimality error for each problem')
[reloptcond]
diary off
hprob3.m0000644001023000102300000000405707045315622013252 0ustar  hwolkowihwolkowiczclear
format compact
format long e

global A D counter mults SHIFTA ;
A=sparse(A);
D=sparse(D);
tcounter = zeros ;        % added by henry - tmults commented out
avecount = zeros ;        % added by henry - tmults commented out
titers = zeros ;        % added by henry - tmults commented out
aveiters = zeros ;        % added by henry - tmults commented out
itermax=0;
N = 16 ;
n = N*N ;
e = ones(n,1) ;
d = e ;
for i = N:N:n
d(i) = 0 ;
end

A = spdiags([ -e -d 2*e ], [-N -1 0 ], n , n) ;
A = (A+A')  - 5*speye(n) ;
tmults = 0 ;
d = diag(A) ;
B = abs(A - spdiags(d,,n,n)) ;
SHIFTA = max(d+sum(B)') ;
OPTIONS.issym = 1; % know the matrix is symmetric
OPTIONS.disp = 0; % no display of the output for eig

[v,d] = eigs('mult',n,1,OPTIONS) ;
d(1)+SHIFTA
v = v/norm(v) ;

for i = 1:20
b = rand(n,1) ;
b = b - v*(v'*b) ;
r   = 100 ;
%   [lopt, xopt, topt, fopt, iter, gap,flag_conj,tarn1f,t0f] = ...
%       newtrust4(b,r) ;
A=sparse(A);
D=sparse(D);
[xopt,fopt,lopt,lambdamin,muup,mulow,iter,flag] = ...
newtrust(b,r,dgaptol);
relfeasxopt(i)=abs(r-norm(xopt))/r;
reloptcond(i)=-norm(  (A-lopt* eye(size(A,1)) )*xopt-b )/fopt;
iter
mults
tmults = tmults + mults ;
xopt = xopt*(r/norm(xopt)) ;
t = b - A*xopt ;
mu = t'*xopt/(r^2) ;
err = norm(t-xopt*mu)
viters(i)=iter;         % total record of iterations
if iter > itermax,
save xprb3 A b r dgaptol viters
itermax=iter;
end

end
avemults = tmults/20
!rm prb3.diary
diary prb3.diary
disp('average of matrix multiplications')
[avecount]
disp('individual iterations per problem')
[viters]
disp('relative feasibility for each problem')
[relfeasxopt]
disp('relative optimality error for each problem')
[reloptcond]
diary off
initialize.m0000644001023000102300000000213607040701060014200 0ustar  hwolkowihwolkowiczfunction [x0] = initialize(n,prob)

% This function initializes x0 for the trust region method depending on the
% problem considered
%
% n is the problem size
% prob is the problem number (refer to objective.m)
%

if prob == 1

x0 = -1*ones(n,1);

end

if prob == 2

x0 = 1/(n+1)*[1:n]';

end

if prob == 3

x0 = 1/n*[1:n]';

end

if prob == 4

div = n/4;

for k=1:div
x0(4*k-3) = 3;
x0(4*k-2) = -1;
x0(4*k-1) = 0;
x0(4*k)   = 1;
end

x0 = x0';

end

if prob == 5

x0 = 1/(n+1)*[1:n]';

end

if prob == 6

x0 = ones(n,1);

end

if prob == 7

x0 = zeros(n,1);

end

if prob == 8

for i=1:n
x0(i) = (i*(i-n-1))/(n+1)^2;
end

x0 = x0';

end

if prob == 9

x0 = -1*ones(n,1);

end

if prob == 10

x0 = zeros(n,1);

end

if prob == 11

x0 = 0.5*ones(n,1);

end

if prob == 12

x0 = ones(n,1);

end

if prob == 13

x0 = ones(n,1);

end

if prob == 14

x0 = ones(n,1);

end

if prob == 15

x0 = zeros(n,1);
x0(n) = 3/2;

end

if prob == 16

x0 = zeros(n,1);

end

if prob == 17

for i=1:n/2
x0(2*i-1) = 2*i;
x0(2*i) = 2*i;
end

x0 = x0';

endmult.m0000644001023000102300000000011707043057314013026 0ustar  hwolkowihwolkowiczfunction [P] = mult(B)

global A D counter mults SHIFTA ;
P = A*B - SHIFTA*B ;
newtrust4.m0000644001023000102300000005027607035530752014042 0ustar  hwolkowihwolkowiczfunction [lopt, xopt, topt, fopt, iter, gap,flag_conj,tarn1f,t0f] = newtrust4(a,s)
% Call:  [lopt, xopt, topt, fopt, iter, gap,flag_conj,tarn1f,t0f] = newtrust4(a,s)
% All elimination of old variables from newtrust.m have been made. A new step
% to the boundary has been included too. Now conjuguate gradient has been
% include in the update phase of the algorithm. 'newtrust.m' uses the

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%WARNING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                          %
% BEFORE WORKING WITH NEWTRUST4 IN MATLAB, THE FOLLOWING INSTRUCTION HAS   %
% TO BE DONE:                                                              %
%                                                                          %
%                            >> global A D counter                         %
%                                                                          %
% THIS TELLS MATLAB THAT BOTH A,D AND COUNTER ARE CONSIDERED GLOBAL IN THE %
% PROGRAM                                                                  %
%                                                                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%WARNING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   Call:  [lopt, xopt, topt, fopt, iter, gap] = newtrust4(a,s)
%
% solves: min      f(x) = x'Ax - 2a'x
%      such that        x'x <= s^2.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%      This code is available by anonymous ftp to orion.uwaterloo.ca
%        and over WWW   URL:   http://orion.uwaterloo.ca/~hwolkowi
%
%      More information on this algorithm can be found in the document
%      'newtrust_guide.ps'. In particular, the equivalent of this matlab
%      file is rewritten and more comments are included. A 'theory' section
%      also explains the different techniques used in the algorithm.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   The user can change the stopping tolerances, for the optimal value,
%      for feasibility, and the iteration count, below.
%   The user can also change the formulae for the various initial
%      bounds, below.
%
%  INPUT:
%         A     ... symmetric matrix - assumed to be in sparse format
%         a     ... column vector - it is changed to full format below
%         s>0   ...   radius of trust region
%
% OUTPUT:
%         lopt ... dual opt multiplier,
%         xopt ... opt. primal solution
%         topt ... right t-value to generate opt. solution;
%                     t=99999 means that lopt=0 used or 0 iterations with no optimal t,
%                       i.e. an unconstrained optimum was found.
%   WARNING: if lopt=0 is found, then finding the unconstrained
%            optimum is not implemented yet. This should be done using
%            conjugate gradients and is left to the user.
%         fopt ... value of objective function at a feasible point,
%              ...   value is < gap(2)+dgaptol if assigned a value;
%              ...   takes value 99999 if not assigned a value.
%
%        iter  ... number of iterations - success occurred if < iterbnd
%        gap   ... the optimal value is in the interval [gap(1,1),gap(1,2)]
%              ... the optimal lambda is in the interval [gap(2,1),gap(2,2)]
%              ... the optimal t value is in the interval [gap(3,1),gap(3,2)]
%              ... the elapsed time for first eig calc is in gap(4,1)
%              ... the elapsed time for the algorithm is in gap(4,2)
%       flag_conj. a boolean variable to know if conjuguated gradient has
%                  been used.

%
% The following statements are put into subprograms since they should
%   be done efficiently, i.e. subprograms used are:
%     Now using the matlab function "eigs" to compute the smallest
%     eigenvalue of a matrix.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%GLOBAL VARIABLES%%%%%%%%%%%%%%%%%%%%%%%
%
global A D counter
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%TOLERANCES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tolerances - these can and should be changed by the user.
%
%    The tolerance for the duality gap is for significant figures,
%     i.e. it is a relative duality gap.
%     The significant figures of agreement are computed for
%       the duality gap as follows:
%      signif = max( -log_10  |primal-dual|/(|primal|+1)  , 0  )
%
%    The tolerances for ***LANCZOS** should be significantly higher
%         than those used here, e.g. 4 here and 7 for Lanczos
%
% Stopping criteria:  (i)   feasibility,      ||x|| <= (1+normtol)s
%                     (ii)  duality gap,      |gap(2)-gap(1)|/(|gap(2)+1)
%                                                               <= dgaptol
%                     (iii) iteration count,  iter > iterbnd
%                     (iv) interval of uncertainty for t is < zerotol
%      Stop when:  ( (i) and (ii) )  or (iii) or (iv)
%
dgaptol  =1e-4;         % duality gap for stopping criteria
normtol  =1e-4;         % allowable error in feasibility for stopping
iterbnd=30;            % upper bound on iterations
n = size( A,1);    % problem size
zerotol=(1e-3)/n;           % tolerance for zero for y(1) for hard case
% since y is norm 1, tolerance decreases with dim.
%%%%%%%%%%%%%%%%%%%%%%%%%%%END TOLERANCES%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%INITIALIZATION%%%%%%%%%%%%%%
t0=cputime;             % saved to get total cputime usage

% Initialize "counter". "counter" gives the number of matrix-vector
% multiplications
counter = 0;

% Initial estimates  -

a=full(a);     %  set to full to avoid conflicts
na=norm(a);      % vector norm

tarn10=cputime; % We want to keep for statistical purpose the first Arnoldi
% cputime

OPTIONS.issym = 1; % know the matrix is symmetric
OPTIONS.disp = 0; % no display of the output for Arnoldi
OPTIONS.tol = 1e-1;
keyboard;
[va,lambda] = eigs('Afun',n,1,'SR',OPTIONS); % computes smallest
%eigenvalue and corresponding eigenvector for A

% lambda is the initial estimate for the optimal lambda of h(lambda).

flag_conj = 0; % flag_conj is 0 the unconstrained minimum is not the answer
% and flag_conj = 1 if it is (conjuguated gradient is then
% used to find it).

tarn1f=cputime-tarn10;   % save initial Arnoldi cputime for statistics
disp(['elapsed time for first Arnoldi',num2str(tarn1f)]);

y=[1;va];     % initial estimate for the eigenvector of D(t^*)
alow=lambda;aup=lambda;

% lower estimate on mu*
mulow = alow*s^2-2*na*s;  % lower estimate for optimal value

bestxf=(s/norm(va))*va; % current best feasible solution

% upper estimate on mu*
if lambda > 0
if (~(abs(va'*a)/(lambda)>s))
muup = -abs(va'*a)/lambda;
else
muup = s^2*lambda - 2*s*abs(va'*a);
end
else
muup = min(0,lambda*s^2-2*s*abs(va'*a));
end

% upper estimate for optimal value-use feasible eigenvector
bestfxf=muup; % current best opt value corresp to bestxf
mubnds=[mulow muup];    % keep history of bounds - latest at top
% This is only for debugging help.
if (lambda > 0)
low = min(alow - na/s,-na*s);
else
low = alow - na/s;  % low: always describes lower bound on t-gets updated
end
high = aup + na*s;  % high: ... upper bound on t-gets updated
%    End of initial estimates
%
%   More initialization of variables

gap=[mulow muup];         % duality gap - gets updated when a bound changes
ss1=s^2+1;
D = [ 0 -a'
-a   A];           % augmented matrix

ystar = 1 / sqrt( s*s +1);      % this is value we are after for y0(t)
iter=0;         % iteration count
b = '      ';              % define blanks
%
% ptgood, ptbad variables: used for saving interpolation data on the
%  good (easy) and bad (hard) side of maximum of k(t).
% The latest data is put into row 1.
%   Only two rows are needed. The rest is saved for debugging help.
%    t y(1) lambda  is data saved.
ptgood=[];

gy=[];by=[];    % The latest vectors y from the good and bad side resp.
side=0;        % indicator for consecutive t found
%  -i for i bad sides and +i for i  good sides
tslopes=[];    % Save slopes of kt on the [bad good] sides - for
%               debugging and for case where kt is almost piecewise linear.
kts=[];      % saving kt bounds and t value - for debugging
bestfxf=muup;         % for first while test below
t=99999;indexg=0;indexb=0;  % initialize in case of 0 iterations
tind=[];     % follow trail of t changes - for debugging
%
disp( 'iteration  - -   t  - -  relative duality gap - - feasibility gap');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%NOTES on ALGORITHM%%%%%%%%%%%%%%
%   Following are the functions of interest:
%    h(lambda)  =  s^2lambda -a'(A-lambda I)^{-1} a;      maximize - concave
%    hp(lambda)  = s^2 -a'(A-lambda I)^{-2} a;solve = 0 - concave, nonincreasing
%    phi(lambda)  = 1/s-1/||(A-lambda I)^{-1} a||;solve = 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    k(t)  =  (s^2+1)lambda -t         maximize - concave
%    k'(t) = (s^2+1)y0^2(t) -1         solve = 0  - same sign as f(t)
%    psi(t)  = -1/y0(t) + sqrt(s^2+1)    solve = 0 - concave, nonincreasing
%    y0(t) = first component of y normalized   - concave, nonincreasing, 1to0
%    lambda=lambda(t)    concave function
%
% The algorithm is based on maximizing k(t).
% We interpolate on psi(t) whenever possible.
%
% The general iteration interpolates psi(t) whenever
%   possible. Duality is used to calculate bounds, i.e. the function k(t) is
%   used when problems with the hard case occur.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%END INITIALIZATION%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%START GENERAL ITERATION%%%%%%%%%%%%%%
%

% Update the tolerances.
t1=1;
t2= ( (bestfxf-muup)/(abs(bestfxf)+1) )   > 2*dgaptol;
t3= ( (muup-mulow)/(abs(muup)+1) )   > dgaptol;
t4=iter <=iterbnd;
t5=abs(high-low)/(abs(high)+1) > zerotol;
%
%keyboard;

%  check feasibility and iteration count and duality gap
while ( (t1&t2) | t3) & (t4 & t5),

iter=iter+1;
indexg=size(ptgood,1);

%  Find a new point t.
t=(low+high)/2;      % midpoint default
tind=[tind 1];    % indicator for t value - for debugging
if indexb > 0 & indexg>0,
% Always find an interpolant using triangle - with an associated upper bound.
thslope=ss1*ptgood(1,2)^2-1;
tslopes=[tlslope thslope
tslopes];         % for debugging purposes - possible use
% when both sides are almost linear
kth=ss1*ptgood(1,3)-ptgood(1,1);     % value at t upper
(thslope*(-ptgood(1,1))+kth)     ];   % rhs of interpolation
ktt=[ 1  (-tlslope)
1  (-thslope)];
kttans=ktt\mrhs;
ktup=kttans(1);
tt=kttans(2);   % This is alternate interpol point using triangle
if low <=tt&tt<=high,  % roundoff error problem check
%  low = tt can occur in hard case and indicates low is almost optimal.
%  An error in divide by zero can occur since iterate is repeated in
%     interpolation.
t=tt;
tind=[tind 2];     % indicator for t value - for debugging
end
if (ktup <=muup) & (ktup >= mulow),
muup=ktup;           % new upper bound on optimal value
end
end     % if indexb > 0
%  End of triangle for new interpolant
%%%%
if side>0,   % last iterate was on the good side

if  indexg>=2,
if (lambda<0) & (min(ptgood(1:2,2)) > (10^(side-1)*zerotol) ),
% two easy case points - prevent too many good side iterates
%   interpolate on good side
tint=ptgood(1:2,1);
yint=ptgood(1:2,2);
f1=-1/yint(1)   + 1/ystar;
f2=-1/yint(2)   + 1/ystar;
%  Slope on the good side should always guarantee f1-f2 not =0.
slope=(tint(1)-tint(2))/(f1-f2);
tt=tint(1) - f1* slope;    % linear inverse interpolation
if low <=tt&tt<=high,     % roundoff error problem check
t=tt;
tind=[tind 3];     % indicator for t value - for debugging
end
elseif side>4,    % We need a point on the bad side - too many
%iterates on good side
t=.95*low+.05*high;
tind=[tind 4];     % indicator for t value - for debugging
end   % end of lambda < 0
end   % end of if indexg>=2

elseif side<0,   % last iterate was on the bad side
if indexb >=2,   % at least two bad side points-
%interpolate from bad side if we have a good side point already.
f1=-1/yint(1)   + 1/ystar;
f2=-1/yint(2)   + 1/ystar;
% If f1=f2, then a t value has been repeated by e.g. t=low in triangle.
% This indicates that the optimal t is very close to t=low.
if f2>f1,   % correct monotonicity occurred
slope=(tint(1)-tint(2))/(f1-f2);
tt=tint(1) - f1* slope;    % linear inverse interpolation
if low<=tt&tt<(high-zerotol),
% Use the new point - it is not too large
t=tt;
tind=[tind 5];     % indicator for t value - for debugging
elseif indexg==0&side>-4,   % need a point on good side
t=.95*high+.05*low;
tind=[tind 6];     % indicator for t value - for debugging
elseif indexg==0&side<-4,   % need a point on good side desperately
t=high;
tind=[tind 7];     % indicator for t value - for debugging
end
elseif side>-4,
t=.95*low + .05*high;  % try to stay on bad side but not too long
tind=[tind 8];     % indicator for t value - for debugging
end
end    %  end for indexb>=2
end    % end for side check
disp(['side: ', num2str(side)]);
disp(['trail for t: ', num2str(tind)]);

%  A new point t has been found. Start updates.
%
D( 1, 1) = t;

tarn20=cputime;

OPTIONS.v0 = y;
[y,lambda] = eigs('Dfun',n+1,1,'SR',OPTIONS);

tarn2f=cputime-tarn20;   % save initial Lanczos call cputime for statistics
disp(['elapsed time for second Arnoldi',num2str(tarn2f)]);

if y(1)<0,  % ensure the sign of the first component
y=-y;
end

%  multiplier sign check - to update lower bounds
if lambda > 0,   % A is known to be pos def - multiplier has wrong sign
%
if y(1) >= ystar,   % means x/y0 is feasible so optimal with lambda=0
%        Finding unconstrained optimum not implemented yet.!!!
a=sparse(a);

gap  = -a'*xopt;   % optimal value since A xopt = a
fopt=gap;
lopt = 0;
topt = 99999;
gap=[topt topt
lopt lopt
topt topt
];
t0f=cputime-t0;
flag_conj = 1; % To know answer is the unconstrained minimum
disp(['elapsed time: ', num2str(t0f)]);
return;    %OPTIMALITY*******************
end
else     % the sign of the multiplier is correct
%   so apply duality theory
kt=ss1*lambda-t;           % value of k(t), dual functional
kts=[t lambda kt
kts];
if (kt>mulow)&(kt<=muup),
mulow=max(mulow,kt);           % new lower bound on optimal value
end
%
end        % end of if lambda > 0 loop
%%  end multiplier sign check

%
if y(1) >= ystar,        % bad side
low = t;           % new lower bound
ptbad];       % Keep points for first interpolation.
by=y;
%%
%%
%%%%%%%%%%%%%
%% Begin bound for PSDP from bad side - approx. negative curvature
%  direction and move to the boundary.
if indexg>0,
if  ptgood(1,3) <= 0,   % multiplier has correct sign on good side
normx2=(1-y(1)^2)/(y(1)^2); % square of norm of current soln.
% Use MoSo formula page 558 - however with negative eigenvalue and eigenvector
xf=y(2:(n+1))/y(1);   % current soln. inside trust region

% We have modified here the original step to the boundary. A different step is
% taken if the easy case or the hard case (case 1) occurs. Otherwise, the ori-
% ginal step is taken.

if (gy(1) < zerotol)
% This is the almost hard case (case 2)
z = gy(2:n+1);
normz=norm(z);
%disp('Hard case primal step');

else
% This is the easy case or almost hard case (case 1)
z = gy(2:n+1)/gy(1,1)-by(2:n+1)/by(1,1);
normz=norm(z);
%disp('Easy case primal step');
end

xyprod=(z'*xf)/normz;  % normalized inner product
ndiff=s^2-normx2;
sxy=sign(xyprod);
if sxy==0,
sxy=1;
end
tau=ndiff/(xyprod+sxy*sqrt(xyprod^2+ndiff) );
xnew=xf+(tau/normz)*z;  %normalized direction
fxnew=xnew'*A*xnew-2*xnew'*a;
% get new upper bound using the feasible xnew - need to allow almost
%  best xnew as the bestxf
if fxnew=mulow),
muup=fxnew;
end

end      % end if <= 0 line mult. check
end        %end of if indexg line
%%%%%%%%%% end of PSDP bound
%
%%%%%%%%%%%%%%%%%%%
else                    %  current point is on good side
high = t;       % we have new upper bound on t
ptgood=[t y(1) lambda
ptgood   ] ;
gy=y;
if  ptgood(1,3) <= 0,   % current multiplier has correct sign
if indexb > 0,
%% begin type 1 bound for PSDP from the good side
normx2=(1-by(1)^2)/(by(1)^2);  % use last bad side point
% Use MoSo formula page 558 - however with negative eigenvalue and eigenvector
xf=by(2:(n+1))/by(1);

% We have modified here the original step to the boundary. A different step is
% taken if the easy case or the hard case (case 1) occurs. Otherwise, the ori-
% ginal step is taken.

if (gy(1) < zerotol)
% This is the almost hard case (case 2)
z = gy(2:n+1);
%disp('hard case primal step');

normz=norm(z);   % This is norm 1 only if real hard case!!!
else
% This is the easy case or the almost hard case (case1)
z = gy(2:n+1)/gy(1,1)-by(2:n+1)/by(1,1);
normz=norm(z);
%disp('easy case primal step');
end

xyprod=(z'*xf)/normz;  % normalized inner product
ndiff=s^2-normx2;
sxy=sign(xyprod);
if sxy==0,
sxy=1;
end
tau=ndiff/(xyprod+sxy*sqrt(xyprod^2+ndiff) );
xnew=xf+(tau/normz)*z;  %normalized direction
fxnew=xnew'*A*xnew-2*xnew'*a;
% get new upper bound using the feasible xnew - need to allow almost
%  best xnew as the bestxf
if fxnew=mulow),
muup=fxnew;
end

else      % else of indexb>0 line
% Need upper bound as well in case triangle was not used.
%
%% begin type 2 bound for PSDP
xf=y(2:(n+1));    % norm is sqrt of 1-y_0^2
xnew=(s/sqrt(1-y(1)^2))*xf;    % feasible xnew
fxnew=xnew'*A*xnew-2*xnew'*a;
% get new upper bound using the feasible xnew
if (fxnew=mulow),
bestxf=xnew;
bestfxf=fxnew;
muup=fxnew;
end
%% reached end type 2 bound for PSDP
end     % ind of if indexb > 0 line
end      % end if <= 0 line - multiplier sign check
end        % end of if bad side
%%%%%%%%%%%%
%% save new bounds and gap
if mulow > mubnds(1,1) | muup < mubnds(1,2),
mubnds= [ mulow muup
mubnds ];
gap=[mulow muup];    % duality gap - gets updated when a bound changes
end
%
% update tolerances
t1=abs( 1-y(1)^2*(1+s^2)) > (normtol*y(1)^2*s^2 );
dbestgap= (bestfxf-muup)/(abs(bestfxf)+1);
t2= dbestgap   > 2*dgaptol;
dgap=  abs(muup-mulow)/(abs(muup)+1) ;
t3= dgap >  dgaptol;
t4=iter <=iterbnd;
t5=abs(high-low)/(abs(high)+1) > zerotol;
txt=[int2str(iter),b,b,b,num2str(t),b,b,num2str(dgap),b,b,num2str(dbestgap)];
disp( txt);

end         % end of while for gap and iteration check
%
%%%%%%%%%%%%%%%%%%%%%%%%%END GENERAL ITERATION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%OUTPUT DATA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fopt=(muup+mulow)/2;
if ~t2,           % bestxf is almost feasible
xopt=bestxf;
fopt=bestfxf;
else
xopt = y(2:n+1)/y(1);
end
topt=t;lopt=lambda;
if indexg ==0,
ptgood=[99999 99999 99999];
end
if indexb ==0,
end
t0f=cputime-t0;
gap=[gap
low high
tarn1f tarn10
];               % tarnf is elapsed time for first eigs call
disp(['elapsed time: ', num2str(t0f)]); % t0f is total computation time
disp(['number of matrix-vector multiplications:', num2str(counter)])
newtrust4b.m0000644001023000102300000005776411350446405014211 0ustar  hwolkowihwolkowicz%%% muup must always correspond to bestxf??? bestfxf???
function [xopt,fopt,lopt,lambdamin,flag] = newtrust4b(a,s,dgaptol)
%%%%%%
% Call:
%[xopt,fopt,lopt,lambdamin,flag] = newtrust4b(a,s,dgaptol)
%%%%%%
%
% solves: min      f(x) = x'Ax - 2a'x
%      such that        x'x <= s^2.
%
%%%%%%%
%% files needed: Afun.m,Dfun.m,conj_grad.m, and this file newtrust4b.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%WARNING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                          %
% BEFORE WORKING WITH NEWTRUST4 IN MATLAB, THE FOLLOWING INSTRUCTION HAS   %
% TO BE DONE:                                                              %
%                                                                          %
%                            >> global A D counter                         %
%                                                                          %
% THIS TELLS MATLAB THAT BOTH A,D AND COUNTER ARE CONSIDERED GLOBAL IN THE %
% PROGRAM                                                                  %
%                                                                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%WARNING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%      This code is available by anonymous ftp to orion.uwaterloo.ca
%        and over WWW   URL:   http://orion.uwaterloo.ca/~hwolkowi
%
%     More information on this algorithm can be found in the document
%     'newtrust_guide.ps'. In particular, the equivalent of this matlab
%     file is rewritten and more comments are included. A 'theory' section
%      also explains the different techniques used in the algorithm.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   The user can change the stopping tolerances, for the optimal value,
%      for feasibility, and the iteration count, below.
%   The user can also change the formulae for the various initial
%      bounds, below.
%
%  INPUT:
%         A     ... symmetric matrix - assumed to be in sparse format
%                            global variable due to Lanczos calls
%         a     ... column vector - it is changed to full format below
%         s>0   ...   radius of trust region
%        dgaptol... duality gap tolerance - changes throughtout the trust
%                   algorithm.
%
% OUTPUT:
%         lopt ... dual opt multiplier vector,
%         xopt ... opt. primal solution vector,
%         fopt ... value of objective function at a feasible point,
%              ...   value is < gap(1,1)+dgaptol if assigned a value;
%              ...   takes value 99999 if not assigned a value.
%       flag   ... a boolean variable to know if solution is in the interior
%                  (flag ==2) and if the hard case (case 2) (or near hard case
%                  (case2)) occurs (flag ==1) and otherwise (flag == 0).
%    lambdamin ... minimum eigenvalue of A
%
% OTHER VARIABLES OF INTEREST:
%
%        iter  ... number of iterations - success occurred if < iterbnd
%        gap   ... gap, optimal value is in interval [gap(1,1),gap(1,2)]
%
%             on output, vector gap is:
%
%             gap=[mulow muup        % lower and upper bounds on mu^* (the optimum
%                                    % of the quadratic objective )
%             ptbad(1,3) ptgood(1,3) % lower and upper bounds on lambda^* (multiplier)
%             low high               % lower and upper bounds on t^*
%             tarn1f t0f]
%       tarn1f ... = cputime-tarn10;  save initial eigenvalue cputime for statistics
%          t0f ... = cputime - t0;  total computation time
%         topt ... right t-value to generate opt. solution;
%                     t=99999 means that lopt=0 used or
%                     0 iterations with no optimal t,
%                       i.e. an unconstrained optimum was found.
%                     t is parameter in eigenvalue problem.
%
%    ptgood, ptbad variables: ... used for saving interpolation data
%                     bestxf  ... current best approximation to the x-optimum
%                     bestfxf  ... current best approximation to the optimum of f(x)
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%GLOBAL VARIABLES%%%%%%%%%%%%%%%%%%%%%%%
%
global A D counter
sparse A D;
if not(issparse(A))| not(issparse(D)),
disp(['A or D not in sparse format in newtrust4b']);
keyboard
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%TOLERANCES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tolerances - these can and should be changed by the user.
%
%    The tolerance for the duality gap is for significant figures,
%     i.e. it is a relative duality gap.
%     The significant figures of agreement are computed for
%       the duality gap as follows:
%      signif = max( -log_10  |primal-dual|/(|primal|+1)  , 0  )
%
%    The tolerances for ***LANCZOS** should be significantly higher
%         than those used here, e.g. 4 here and 7 for Lanczos
%
% Stopping criteria:  (i)   feasibility,      ||x|| <= (1+normtol)s
%            ???      (ii)  ????????,      dbestgap= (bestfxf-muup)/(|bestfxf|+1)
%            ???                                                <= 2*dgaptol
%                     (iii)  duality gap,      |gap(2)-gap(1)|/(|gap(2)|+1)
%                                                               <= dgaptol
%                     (iv) iteration count,  iter > iterbnd
%                     (v) interval of uncertainty for t is < zerotol
%      Stop when:  ( ((i) and (ii)) or (iii))  and ((iv) and (v))
%
%dgaptol                % duality gap for stopping criteria
normtol  =1e-9;         % allowable error in feasibility for stopping
iterbnd=30;             % upper bound on iterations
n = size( A,1);         % problem size
zerotol=(1e-4)/n;       % tolerance for zero for y(1) for hard case
% since y is norm 1, tolerance decreases with dim.
%%%%%%%%%%%%%%%%%%%%%%%%%%%END TOLERANCES%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%INITIALIZATION%%%%%%%%%%%%%%
t0=cputime;             % saved to get total cputime usage

% Initialize "counter". "counter" gives the number of matrix-vector
% multiplications
counter = 0;

% Initial estimates  -

a=full(a);     %  set to full to avoid conflicts
na=norm(a);      % vector norm

tarn10=cputime; % We want to keep for statistical purpose the first eig calc
% cputime

OPTIONS.issym = 1; % because we know the matrices A and D are symmetric
OPTIONS.disp = 0; % no display of the output for eig.m
%v=min(diag(A));
%d=diag(diag(A)+2*(abs(v)+1)*ones(n,1));
%v=d\a;
v=full(sum(A))';       % good/cheap estimate for sparse A ???
OPTIONS.v0 =v ; % initial eigenvector estimate

%%%%% adding this can accelerate the computations of the first eigenvalue
%lambda = min(eig(A));
%R = A;
%for i=1:n
%	R(i,i) = R(i,i) -lambda;
%end
%R = sparse(R);

%[va, lambda2] = eigs(R,1,0,OPTIONS);

%%%% you also need to put the next line as a comment

[va,lambda] = eigs('Afun',n,1,'SA',OPTIONS); % computes smallest
%eigenvalue (lambda) and corresponding eigenvector (va) for A

lambdamin = lambda; % save the min eigenvalue of A

flag = 0; %flag will change value if the almost hard case (case 2) occurs or if
%the solution is in the interior.

tarn1f=cputime-tarn10;   % save initial eig calc cputime for statistics
disp(['elapsed time for first eig calc',num2str(tarn1f)]);

y=[1;va];     % initial estimate for the eigenvector of D(t^*)
alow=lambda;aup=lambda;

%two lines replaced by following 6 lines jan 10/00
%% lower estimate on mu*
%mulow = alow*s^2-2*na*s;  % lower estimate for optimal value

% lower estimate on mu*
if lambda<0
mulow = alow*s^2-2*na*s;  % lower estimate for optimal value
else
mulow = -2*na*s;
end

bestxf=(s/norm(va))*va; % current best feasible solution

% upper estimate on mu* -optimal value- uses feasible eigenvector

% ensure va'*a is not negative
if (va'*a < 0)
va = -va;
end

% upper estimate on mu*
if lambda > 0
if (~((va'*a)/(lambda)>s))
muup = -(va'*a)/lambda;
bestxf = -((va'*a)/lambda)*va;
else
muup = s^2*lambda - 2*s*(va'*a);
bestxf = s*va;
end
else
[muup,indice] = min([0,lambda*s^2-2*s*(va'*a)]);
if indice == 1
bestxf = zeros(size(a,1),1);
else
bestxf = s*va;
end
end

bestfxf=muup; % current best opt value corresp to bestxf
mubnds=[mulow muup];    % keep history of bounds - latest at top
% This is only for debugging help.
if (lambda > 0)
low = min(alow - na/s,-na*s);
else
low = alow - na/s;  % low: always describes lower bound on t-gets updated
end

high = aup + na*s;  % high: ... upper bound on t-gets updated

%    End of initial estimates
%
%   More initialization of variables

gap=[mulow muup];         % duality gap - gets updated when a bound changes
ss1=s^2+1;
D = [ 0 -a'
-a   A];           % augmented matrix

ystar = 1 / sqrt( s*s +1);      % this is value we are after for y0(t)
iter=0;         % iteration count
b = '      ';              % define blanks
%
% ptgood, ptbad variables: used for saving interpolation data on the
%  good (easy) and bad (hard) side of maximum of k(t).
% The latest data is put into row 1.
%   Only two rows are needed. The rest is saved for debugging help.
%    t y(1) lambda  is data saved.
ptgood=[];

gy=[];by=[];    % The latest vectors y from the good and bad side resp.
side=0;        % indicator for consecutive t found
%  -i for i bad sides and +i for i  good sides
tslopes=[];    % Save slopes of kt on the [bad good] sides - for
%               debugging and for case where kt is almost piecewise linear.
kts=[];      % saving kt bounds and t value - for debugging
bestfxf=muup;         % for first while test below
t=99999;indexg=0;indexb=0;  % initialize in case of 0 iterations
tind=[];     % follow trail of t changes - for debugging
%
disp( 'iteration  - -   t  - -  relative duality gap - - feasibility gap');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%NOTES on ALGORITHM%%%%%%%%%%%%%%
%   Following are the functions of interest:
%    h(lambda)  =  s^2lambda -a'(A-lambda I)^{-1} a;      maximize - concave
%    hp(lambda)  = s^2 -a'(A-lambda I)^{-2} a;solve = 0 - concave, nonincreasing
%    phi(lambda)  = 1/s-1/||(A-lambda I)^{-1} a||;solve = 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    k(t)  =  (s^2+1)lambda -t         maximize - concave
%    k'(t) = (s^2+1)y0^2(t) -1         solve = 0  - same sign as f(t)
%    psi(t)  = -1/y0(t) + sqrt(s^2+1)    solve = 0 - concave, nonincreasing
%    y0(t) = first component of y normalized   - nonincreasing, 1to0
%    lambda=lambda(t)    concave function
%
% The algorithm is based on maximizing k(t).
% We interpolate on psi(t) whenever possible.
%
% The general iteration interpolates psi(t) whenever
%   possible. Duality is used to calculate bounds, i.e. the function k(t) is
%   used when problems with the hard case occur.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%END INITIALIZATION%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%START GENERAL ITERATION%%%%%%%%%%%%%%
%

% Update the tolerances.
t1=1;    % feasibility tolerance using y(1)
t2= ( (bestfxf-muup)/(abs(bestfxf)+1) )   > 2*dgaptol; %%???why???
t3= ( (muup-mulow)/(abs(muup)+1) )   > dgaptol; % duality gap tolerance
t4=iter <=iterbnd; % iteration bound
t5=abs(high-low)/(abs(high)+1) > zerotol; %interval of uncertainty for t

%  check feasibility and iteration count and duality gap
while ( (t1&t2) | t3) & (t4 & t5),

iter=iter+1; % one more iteration is done
indexg=size(ptgood,1); % number of points on the good side

%  Find a new point t.
t=(low+high)/2;      % midpoint default
tind=[tind 1];    % indicator for t value - for debugging

if indexb > 0 & indexg > 0,
% TRIANGLE INTERPOLATION

% Always find an interpolant using triangle - with an associated upper bound.
thslope=ss1*ptgood(1,2)^2-1;% slope of k(t) at the last good side point
tslopes=[tlslope thslope
tslopes];         % for debugging purposes - possible use
% when both sides are almost linear
kth=ss1*ptgood(1,3)-ptgood(1,1); % value of k(t)at last good side point
(thslope*(-ptgood(1,1))+kth)     ];   % rhs of interpolation
ktt=[ 1  (-tlslope)
1  (-thslope)];
kttans=ktt\mrhs; % This finds the intersection of the two tangent lines at
% the last bad and good side points. The point of
% intersection is (kttans(2),kttans(1)). Hence, kttans(1)
% is an upper bound on mu^* and kttans(2) is an approximation
% for t^*
ktup=kttans(1);
tt=kttans(2);   % This is alternate interpol point using triangle
if low <=tt&tt<=high,  % roundoff error problem check
%  low = tt can occur in hard case and indicates low is almost optimal.
%  An error in divide by zero can occur since iterate is repeated in
%     interpolation.
t=tt; % new update for the approximation of t^*
tind=[tind 2];     % indicator for t value - for debugging
end
if (ktup <=muup) & (ktup >= mulow),
muup=ktup;           % new upper bound on optimal value
end
end     % if indexb > 0
%  End of triangle for new interpolant

if side>0,   % last iterate was on the good side

if  indexg>=2, % At least two t-iterates are on the good side
if (lambda<0) & (min(ptgood(1:2,2)) > (10^(side-1)*zerotol) ),
% two easy case points - prevent too many good side iterates
% INVERSE INTERPOLATION - interpolate on good side
% Let t(psi) be the inverse function of psi(t), then we linearly interpolate
% the value t(0), which is t^*, using the two points from the good side.
tint=ptgood(1:2,1);
yint=ptgood(1:2,2);
f1=-1/yint(1)   + 1/ystar;
f2=-1/yint(2)   + 1/ystar;
%  Slope on the good side should always guarantee f1-f2 not =0.
slope=(tint(1)-tint(2))/(f1-f2);
tt=tint(1) - f1* slope;    % linear inverse interpolation
if low <=tt&tt<=high,     % roundoff error problem check
t=tt;
tind=[tind 3];     % indicator for t value - for debugging
end
elseif side>4,    % We need a point on the bad side - too many
%iterates on good side
t=.95*low+.05*high;
tind=[tind 4];     % indicator for t value - for debugging
end   % end of lambda < 0
end   % end of if indexg>=2

elseif side<0,   % last iterate was on the bad side
if indexb >=2,   % at least two bad side points-
%INVERSE INTERPOLATION
% Let t(psi) be the inverse function of psi(t), then we linearly interpolate
% the value t(0), which is t^*, using the two points from the bad side.

%?????  %interpolate from bad side if we have a good side point already.

f1=-1/yint(1)   + 1/ystar;
f2=-1/yint(2)   + 1/ystar;
% If f1=f2, then a t value has been repeated by e.g. t=low in triangle.
% This indicates that the optimal t is very close to t=low.
if f2>f1,   % correct monotonicity occurred
slope=(tint(1)-tint(2))/(f1-f2);
tt=tint(1) - f1* slope;    % linear inverse interpolation
if low<=tt&tt<(high-zerotol),
% Use the new point - it is not too large
t=tt;
tind=[tind 5];     % indicator for t value - for debugging
elseif indexg==0&side>-4,   % need a point on good side
t=.95*high+.05*low;
tind=[tind 6];     % indicator for t value - for debugging
elseif indexg==0&side<-4,   % need a point on good side desperately
t=high;
tind=[tind 7];     % indicator for t value - for debugging
end
elseif side>-4,
t=.95*low + .05*high;  % try to stay on bad side but not too long
tind=[tind 8];     % indicator for t value - for debugging
end
end    %  end for indexb>=2
end    % end for side check
disp(['side: ', num2str(side)]);
disp(['trail for t: ', num2str(tind)]);

%  A new point t has been found. Start updates.
%
% Recall t is the most recent update for t^*

D( 1, 1) = t;

tarn20=cputime;

OPTIONS.v0 = y; % y is the most recent estimate for the eigenvector of D(t^*)
[y,lambda] = eigs('Dfun',n+1,1,'SA',OPTIONS); % y and lambda are the new estimates
% for the eigenvector of D(t^*) and for lambda^*

tarn2f=cputime-tarn20;   % save initial Lanczos call cputime for statistics
disp(['elapsed time for second eig calc',num2str(tarn2f)]);

if y(1)<0,  % ensure the sign of the first component is positive
y=-y;
end

%  multiplier sign check - to update lower bounds
if lambda > 0,   % A is known to be pos def - multiplier has wrong sign

if y(1) >= ystar,
% This means that y(2:n+1)/y(1) is the optimum solution for all points outside
% the trust region of radius (1-y(1)^2)/y(1)^2 or on its boundary. If this
% radius is less than s^2, as it is the case here since this is equivalent to
% y(1) >= ystar, then this means that the optimum we are looking for is inside
% the trust region and A is positive semidefinite. We will then apply a conjugate
% gradient method to find the optimum. Hence, y(2:n+1)/y(1) is feasible so
% optimal with lambda=0.

a=sparse(a);

% find the optimum.

gap  = -a'*xopt;   % optimal value since A xopt = a
fopt=gap;
lopt = 0;
topt = 99999;
gap=[topt topt
lopt lopt
topt topt
];
t0f=cputime-t0;
flag = 2; % To know answer is the unconstrained minimum
disp(['elapsed time: ', num2str(t0f)]);
return;    %OPTIMALITY - PROBLEM IS SOLVED
end
else     % the sign of the multiplier is correct
%   so apply duality theory
kt=ss1*lambda-t;           % value of k(t), dual functional. This gives
% a lower bound on mu^*
kts=[t lambda kt
kts];
if (kt>mulow)&(kt<=muup),
mulow=max(mulow,kt);           % new lower bound on optimal value
end
%
end        % end of if lambda > 0 loop
%%  end multiplier sign check

%
if y(1) >= ystar,        % bad side
low = t;           % new lower bound on t^*
ptbad];       % Keep points for first interpolation.
by=y;
%%
%%
%%%%%%%%%%%%%
%% Begin bound for PSDP from bad side - approx. negative curvature
%  direction and move to the boundary.
if indexg>0,
if  ptgood(1,3) <= 0,   % multiplier has correct sign on good side
normx2=(1-y(1)^2)/(y(1)^2); % square of norm of current soln.
xf=y(2:(n+1))/y(1);   % current soln. inside trust region

% In the hard case (case 2), we use MoSo formula page 558 - however
% with negative eigenvalue and eigenvector, to move to take a step to
% the boundary. In the easy case and hard case (case 1), we take a step
% to the boundary by taking the solution between the line joining the
% xf and the x-point - gy(2:n+1)/gy(1,1) -  outside the trust region
% associate with the last good side t-point.

if (gy(1) < zerotol)
% This is the almost hard case (case 2)
z = gy(2:n+1);
normz=norm(z); % normz should be almost 1

else
% This is the easy case or almost hard case (case 1)
z = gy(2:n+1)/gy(1,1)-by(2:n+1)/by(1,1);
normz=norm(z);

end

xyprod=(z'*xf)/normz;  % normalized inner product
ndiff=s^2-normx2;
sxy=sign(xyprod);
if sxy==0,
sxy=1;
end
tau=ndiff/(xyprod+sxy*sqrt(xyprod^2+ndiff) );
xnew=xf+(tau/normz)*z;  % xnew is the point obtained by moving to the
% boundary
fxnew=xnew'*A*xnew-2*xnew'*a;
% get new upper bound using the feasible xnew - need to allow almost
%  best xnew as the bestxf
if fxnew=mulow),
muup=fxnew;
end

end      % end if <= 0 line mult. check
end        %end of if indexg line
%%%%%%%%%% end of PSDP bound
%
%%%%%%%%%%%%%%%%%%%
else                    %  current point is on good side
high = t;       % we have new upper bound on t
ptgood=[t y(1) lambda
ptgood   ] ;
gy=y;
if  ptgood(1,3) <= 0,   % current multiplier has correct sign
if indexb > 0,
%% begin type 1 bound for PSDP from the good side
normx2=(1-by(1)^2)/(by(1)^2);  % use last bad side point
xf=by(2:(n+1))/by(1);

% In the hard case (case 2), we use MoSo formula page 558 - however
% with negative eigenvalue and eigenvector, to move to take a step to
% the boundary. In the easy case and hard case (case 1), we take a step
% to the boundary by taking the solution between the line joining the
% xf and the x-point - gy(2:n+1)/gy(1,1) -  outside the trust region
% associate with the last good side t-point.

if (gy(1) < zerotol)
% This is the almost hard case (case 2)
z = gy(2:n+1);
normz=norm(z); %normz
else
% This is the easy case or the almost hard case (case1)
z = gy(2:n+1)/gy(1,1)-by(2:n+1)/by(1,1);
normz=norm(z);
end

xyprod=(z'*xf)/normz;  % normalized inner product
ndiff=s^2-normx2;
sxy=sign(xyprod);
if sxy==0,
sxy=1;
end
tau=ndiff/(xyprod+sxy*sqrt(xyprod^2+ndiff) );
xnew=xf+(tau/normz)*z;  %
fxnew=xnew'*A*xnew-2*xnew'*a;
% get new upper bound using the feasible xnew - need to allow almost
%  best xnew as the bestxf
if fxnew=mulow),
muup=fxnew;
end

else      % else of indexb>0 line
% Need upper bound as well in case triangle was not used.
%
%% begin type 2 bound for PSDP
xf=y(2:(n+1));    % norm is sqrt of 1-y_0^2
xnew=(s/sqrt(1-y(1)^2))*xf; % xnew is the point obtained by moving to the
% boundary
fxnew=xnew'*A*xnew-2*xnew'*a;
% get new upper bound using the feasible xnew
if (fxnew=mulow),
bestxf=xnew;
bestfxf=fxnew;
muup=fxnew;
end
%% reached end type 2 bound for PSDP
end     % ind of if indexb > 0 line
end      % end if <= 0 line - multiplier sign check
end        % end of if bad side
%%%%%%%%%%%%
%% save new bounds and gap
if mulow > mubnds(1,1) | muup < mubnds(1,2),
mubnds= [ mulow muup
mubnds ];
gap=[mulow muup];    % duality gap - gets updated when a bound changes
end
%
% update tolerances
t1=abs( 1-y(1)^2*(1+s^2)) > (normtol*y(1)^2*s^2 ); % Feasibility criteria
dbestgap= (bestfxf-muup)/(abs(bestfxf)+1);
t2= dbestgap   > 2*dgaptol; %??????
dgap=  abs(muup-mulow)/(abs(muup)+1) ;
t3= dgap >  dgaptol; % Duality gap criteria
t4=iter <=iterbnd; % iteration bound criteria
t5=abs(high-low)/(abs(high)+1) > zerotol; % gap for t^* criteria
txt=[int2str(iter),b,b,b,num2str(t),b,b,num2str(dgap),b,b,num2str(dbestgap)];
disp( txt);

end         % end of while for gap and iteration check
%
%%%%%%%%%%%%%%%%%%%%%%%%%END GENERAL ITERATION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%OUTPUT DATA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fopt=(muup+mulow)/2;
if ~t2,           % bestxf is almost feasible
xopt=bestxf;
fopt=bestfxf;
else
xopt = y(2:n+1)/y(1);
end
topt=t;lopt=lambda;
if indexg ==0,
ptgood=[99999 99999 99999];
end
if indexb ==0,
end
t0f=cputime-t0;
gap=[gap                    % lower and upper bound for mu^*
ptbad(1,3) ptgood(1,3) % lower and upper bound for lambda^*
low high               % lower and upper bound for t^*
tarn1f t0f ];          % tarn1f is elapsed time for first eigs call
% t0f is total computation time

disp(['elapsed time: ', num2str(t0f)]); % t0f is total computation time
disp(['number of matrix-vector multiplications:', num2str(counter)])

if (abs(lambdamin - lopt)/(abs(lambdamin)+1) < 1e-2)
flag = 1;
end
ntrust.m0000644001023000102300000001257107046632145013417 0ustar  hwolkowihwolkowiczfunction [x,histout,costdata] = ntrust(x0,f,tol,maxit,resolution)
%
%
% function [x,histout,costdata] = ntrust(x0,f,tol,maxit,resolution)
%
%        trust region, Newton model, sparse algorithm
%
% Input: x0 = initial iterate
%        f = objective function,
%            the calling sequence for f should be
%            [fout,gout]=f(x) where fout=f(x) is a scalar
%              and gout = grad f(x) is a COLUMN vector
%        tol = termination criterion norm(grad) < tol
%        maxit = maximum iterations (optional) default = 100
%        resolution = estimated accuracy in functions/gradients (optional)
%                     default = 1.d-12
%                     The finite difference increment in the difference
%                     Hessian is set to sqrt(resolution).
%
%
% Output: x = solution
%         histout = iteration history
%             Each row of histout is
%         costdata = [num f, num grad, num hess]
%
% Requires: diffhess.m, dirdero.m
%
% set maxit, the resolution, and the difference increment
%
debug=0;
if nargin < 4
maxit = 100;
end
if nargin < 5
resolution = 1.d-12;
end
hdiff=sqrt(resolution);
%
maxit=100; itc=1; xc=x0; n=length(x0);
[fc,gc]=feval(f,xc);
numf=1; numg=1; numh=0;
ithist=zeros(maxit,4);
ithist(1,1)=norm(gc); ithist(1,2) = fc; ithist(1,4)=itc-1;
ithist(1,3)=0;
if debug == 1
ithist(itc,:)
end
%
% Iniitalize the TR radius, not a profound choice.

%
ijob=1;
jdata=[];
while(norm(gc) > tol & itc <= maxit)
if ijob == 1
hess=diffhess(xc,f,gc,hdiff);
numf=numf+n; numg=numg+n; numh=numh+1;
else
jdata=sdata;
end
itc=itc+1;
numf=numf+nf;
ijob=idid;
if idid == 2
sdata=jdata;
elseif idid == 3
elseif idid == 4
end
if idid==1
[fc,gc]=feval(f,xc);
numf=numf+1; numg=numg+1;
end
ithist(itc,1)=norm(gc); ithist(itc,2) = fc;
if debug == 1
ithist(itc,:)
end
end
x=xc; histout=ithist(1:itc,:);
costdata=[numf, numg, numh];
trfix(xc, f, hc, gc, fc, oldrad,ijob,jdata)
%
%
%     C. T. Kelley, Dec 15, 1997
%
%     This code comes with no guarantee or warranty of any kind.
%
%     function [xt, newrad, idid, sdata, nf]
%                = trfix(xc, f, hc, gc, oldrad, ijob, jdata)
%
%     Figure out what the new trust region radius and new point are
%
%     This code is called by ntrust.m
%     There is no reason for you to call this directly.
%
%     Input: xc = current point
%     f  = objective
%     hc = current Hessian
%     fc = current function value
%     ijob = what to do now: 1 = fresh start
%                            2 = TR radius reduction in progress
%                            3 = attempt TR radius expansion
%     jdata = Newton direction when ijob = 1 or 2, avoids recomputation
%     nf = number of function evaluations
%
%     Output: xp = new point
%     idid = result flag: 1 = step ok
%                         2 = TR radius reduced, step not ok
%                         3 = expansion attempted, save step and try
%                             to do better
%                         4 = expansion step failed, use the last good step
%     sdata = Newton direction to use in next call if idid > 1
%
%
%     Find the Cauchy point
%
%     bflag=1 means that the trial point is on the TR boundary and is not
%             the Newton point
%
nf=0;
bflag=0;
idid=1;
mu=gc'*hc*gc;
mu1=gc'*gc;
dsd=-gc;
if ijob == 1
dnewt=hc\dsd;
else
dnewt=jdata;
end
sdata=dnewt;
if mu > 0
sigma = mu1/mu;
end
cp = xc-sigma*gc;
else
%
%     If steepest descent direction is a direction of negative curvature
%     take a flying leap to the boundary.
%
bflag=1;
cp=xc-sigma*gc;
end
%
%     If CP is on the TR boundary, that's the trial point.
%     If it's not, compute the Newton point and complete the dogleg.
%
if bflag==1
xt=cp;
else
%
%     If we get to this point, CP is in the interior and the steepest
%     descent direction is a direction of positive curvature.
%
dsd=-gc; dnewt=hc\dsd;
xn=xc+dnewt;
mu2=dsd'*dnewt;
%
%     If the Newton direction goes uphill, revert to CP.
%
if mu2 <= 0
xt=cp;
%
%     If the Newton point is inside, take it.
%
xt=xn;
%
%    Newton is outside and CP is inside. Find the intersection of the
%    dog leg path with TR boundary.
%
else
d1=sigma*gc; d2=d1+dnewt;
xi=(-bco+sqrt((bco*bco) - 4*aco*cco))/(2*aco);
xt=cp + xi*(xn-cp);
bflag=1;
end
end
%
%
st=xt-xc; ft=feval(f,xt); ared=ft-fc; nf=nf+1;
pred=gc'*st + .5* (st'*hc*st);
if ared/pred < .25
xt=xc;
idid=2;
if ijob == 3 idid = 4; end
elseif ared/pred > .75 & bflag==1
idid=3;
end
xp=xt;
objective.m0000644001023000102300000001351407043101560014015 0ustar  hwolkowihwolkowiczfunction F = objective(x,prob)

% This function returns the objective value F(x) of F associated with problem
% number "prob"
%
% Here are the problems and their related number
%
% 1 = Broyden banded function BRYBND (31) in 
% 2 = Generalized Rosenbrock function GENROSE 
% 21 = GENROSE, but a scalar of 1000 is used instead in the function
% 22 = A function close to GENROSE (no reference)
% 23 = generalization of the banana function (no reference)
% 3 = Optimal Sensor Placement SENSORS 
% 4 = Extended Powell Singular Function (22) in 
% 5 = NONDIA in 
% 6 = TRIDIA in 
% 7 = EXTROS in 
% 8 = Discrete Boundary Value Function (28) in 
% 9 = Broyden Tridiagonal Function (30) in 
% 10= Watson Function (20) in 
% 11= Brown Almost-Linear Function (27) in 
%
% Note that for the last 2 problems, 'm' can be variable. It is fixed to 2*n
% in the code right now.
%
% 12= Linear Function-full rank (32) in 
% 13= Linear function-rank 1 (33) in 
% 14= Linear function-rank 1 with zero columns and rows (34) in 
% 15= A function I created so the hard case (case 2) would occur initially
% 16 & 17 = two other function where the hard case occurs
%
%
%  = J.J. More, B.S. Garbow and K.E. Hillstrom, Testing unconstrained
%       optimization software, 7(1981), no 1, pp. 17-41.
%
%  = S. Nash, Newton-type minimization via the Lanczos process, SIAM journal
%       on numerical analysis, 21, 1884, pp. 770-788.
%
%  = H. Zhang, X. Wang, Optimal Sensor placement, SIAM Review, vol. 35,
%       pp. 641, 1993.
%
%  = G. Liu and J. Nocedal, Test results of two limited memory methods for
%       large scale optimization, Report NAM 04, Department of Electrical
%       Engineering and Computer Science, Northwestern University, Evanston,
%       Illinois.

n = size(x,1); % x is assumed to be a column vector;

if prob == 1

if (n < 6)
disp('Error: n has to be greater or equal to 7');
F = 999999;
return;
end

f(1) = x(1)*(2+5*x(1)^2) + 1 - x(2)*(1+x(2));
f(2) = x(2)*(2+5*x(2)^2) + 1 - x(1)*(1+x(1))-x(3)*(1+x(3));
f(3) = x(3)*(2+5*x(3)^2) + 1 - sum(x(1:2).*(ones(2,1)+x(1:2)))-x(4)*(1+x(4));
f(4) = x(4)*(2+5*x(4)^2) + 1 - sum(x(1:3).*(ones(3,1)+x(1:3)))-x(5)*(1+x(5));
f(5) = x(5)*(2+5*x(5)^2) + 1 - sum(x(1:4).*(ones(4,1)+x(1:4)))-x(6)*(1+x(6));
f(n) = x(n)*(2+5*x(n)^2) + 1 - sum(x(n-5:n-1).*(ones(5,1)+x(n-5:n-1)));

for i=6:n-1
f(i) = x(i)*(2+5*x(i)^2) + 1 - sum(x(i-5:i-1).*(ones(5,1)+x(i-5:i-1)))-x(i+1)*(1+x(i+1));
end

F = sum(f.^2);

end

if prob == 2

for i=2:n
f(i) = 100*(x(i)-x(i-1)^2)^2 + (1-x(i))^2;
end

F = sum(f);

end

if prob == 21

const = 1000;

for i=2:n
f(i) = const*(x(i)-x(i-1)^2)^2 + (1-x(i))^2;
end

F = sum(f);

end

if prob == 22

for i=2:n
f(i) = 100*(x(i)-x(i-1))^2 + (1-x(i))^2;
end

F = sum(f);

end

if prob == 23

for i=2:n
f(i) = 100*(x(i)-x(i-1)^2)^2 + (1-x(i-1))^2;
end

F = sum(f);

end

if prob == 3

for i=1:n
for j=1:n
f(i,j) = sin(x(i))^2*sin(x(j))^2*sin(x(i)-x(j))^2;
end
end

F = -sum(sum(f));

end

if prob == 4

if ((n/4 ~= floor(n/4)) & (n/4 ~= ceil(n/4)))
disp('Error: n has to be a multiple of 4');
F = 999999;
return;
end

div = n/4;

for k=1:div
f(4*k-3) = x(4*k -3) + 10*x(4*k -2);
f(4*k-2) = sqrt(5)*(x(4*k -1)-x(4*k));
f(4*k-1) = (x(4*k-2)-2*x(4*k-1))^2;
f(4*k)   = sqrt(10)*(x(4*k-3) -x(4*k))^2;
end

F = sum(f.^2);

end

if prob == 5

for i=1:n
f(i) = 100*(x(1)-x(i)^2)^2+(1-x(i))^2;
end

F = sum(f);

end

if prob == 6

for i=2:n
f(i) = sqrt(i)*(2*x(i)-x(i-1));
end

F = sum(f.^2);

end

if prob == 7

if ((n/2 ~= floor(n/2)) & (n/2 ~= ceil(n/2)))
disp('Error: n has to be a multiple of 2');
F = 999999;
return;
end

div = n/2;

for k=1:div
f(k) = 100*(x(2*k)-x(2*k-1)^2)^2+(1-x(2*k-1))^2;
end

F = sum(f);

end

if prob == 8

if (n<3)
disp('Error: n has to be greater or equal to 3');
F = 999999;
return;
end

f(1) = 2*x(1) -x(2) + (x(1) + (n+2)/(n+1))^3/(2*(n+1)^2);
f(n) = 2*x(n) -x(n-1) + (x(n) + (2*n+1)/(n+1))^3/(2*(n+1)^2);

for i=2:n-1
f(i) = 2*x(i) -x(i-1) -x(i+1) +(x(i) + (n+i+1)/(n+1))^3/(2*(n+1)^2);
end

F = sum(f.^2);

end

if prob == 9

f(1) = (3-2*x(1))*x(1) -2*x(2) +1;
f(n) = (3-2*x(n))*x(n) -x(n-1) +1;

for i=2:n-1
f(i) = (3-2*x(i))*x(i) -x(i-1) -2*x(i+1) +1;
end

F = sum(f.^2);

end

if prob == 10

if (n~=31)
disp('Error: n has to be equal to 31');
F = 999999;
return;
end

index = [1:n-1]';

f(30) = x(1);
f(31) = x(2)-x(1)^2-1;
for i=1:29
f(i) = sum(index.*x(2:n).*((i/29*ones(n-1,1)).^(index-ones(n-1,1))) )-...
(sum(x.*((i/29*ones(n,1)).^([0:n-1]'))))^2 - 1;
end

F = sum(f.^2);

end

if prob == 11

for i=1:n-1
f(i) = x(i) + sum(x) - (n+1);
end
f(n) = prod(x) -1;

F = sum(f.^2);

end

if prob == 12

m = 2*n;

for i=1:n
f(i) = x(i) -2/m*sum(x) -1;
end

for i=n+1:m
f(i) = -2/m*sum(x) -1;
end

F = sum(f.^2);

end

if prob == 13

m=2*n;

for i=1:m
f(i) = i*sum(([1:n]').*x)-1;
end

F = sum(f.^2);

end

if prob == 14

m = 2*n;

f(1) = -1;
f(m) = -1;

for i=2:m-1
f(i) = (i-1)*sum(([2:n-1]').*x(2:n-1)) -1;
end

F = sum(f.^2);

end

if prob == 15

f(n) = (x(n) - 1)^2;

for i=1:n-1
f(i) = (x(i)^2-1)^2;
end

F = sum(f);

end

if prob == 16

if ((n/2 ~= floor(n/2)) & (n/2 ~= ceil(n/2)))
disp('Error: n has to be a multiple of 2');
F = 999999;
return;
end

F = 0;
a = 0.4238537991;

for i=1:n/2
F = F+(x(2*i-1)+a)^2*(x(2*i)+a)^2+(x(2*i-1)+a)^2+(x(2*i)+a)^2;
end

end

if prob == 17
if ((n/2 ~= floor(n/2)) & (n/2 ~= ceil(n/2)))
disp('Error: n has to be a multiple of 2');
F = 999999;
return;
end

F=0;

for i=1:n/2
F = F + (exp(x(2*i-1)+x(2*i))-1)^2;
end
end

% This function returns the objective value F(x) of F associated with problem
% number "prob", the gradient and the Hessien (returned in sparse format).
%
% Here are the problems and their related number
%
% 1 = Broyden banded function BRYBND (31) in 
% 2 = Generalized Rosenbrock function GENROSE 
% 21 = GENROSE, but a scalar of 1000 is used instead in the function
% 22 = a function close to GENROSE (no reference)
% 23 = generalization of the banana function (no reference)
% 3 = Optimal Sensor Placement SENSORS 
% 4 = Extended Powell Singular Function PWSING (22) in 
% 5 = NONDIA in 
% 6 = TRIDIA in 
% 7 = EXTROS in 
% 8 = Discrete Boundary Value Function (28) in 
% 9= Broyden Tridiagonal Function (30) in 
% 10= Watson Function (20) in 
% 11= Brown Almost-Linear Function (27) in 
%
% Note that for the last 2 problems, 'm' can be variable. It is fixed to 2*n
% in the code right now.
%
% 12= Linear Function-full rank (32) in 
% 13= Linear function-rank 1 (33) in 
% 14= Linear function-rank 1 with zero columns and rows (34) in 
% 15= A function I created so the hard case (case 2) would occur initially
%
%  = J.J. More, B.S. Garbow and K.E. Hillstrom, Testing unconstrained
%       optimization software, 7(1981), no 1, pp. 17-41.
%
%  = S. Nash, Newton-type minimization via the Lanczos process, SIAM journal
%       on numerical analysis, 21, 1884, pp. 770-788.
%
%  = H. Zhang, X. Wang, Optimal Sensor placement, SIAM Review, vol. 35,
%       pp. 641, 1993.
%
%  = G. Liu and J. Nocedal, Test results of two limited memory methods for
%       large scale optimization, Report NAM 04, Department of Electrical
%       Engineering and Computer Science, Northwestern University, Evanston,
%       Illinois.

n = size(x,1); % x is assumed to be a column vector;

if prob == 1

if (n < 6)
disp('Error: n has to be greater or equal to 7');
F = 999999;
return;
end

%  J(1,:) = [2,0,0,0,0,0];
%  J(2,:) = [1,3,0,0,0,0];
%  J(3,:) = [1,2,4,0,0,0];
%  J(4,:) = [1,2,3,5,0,0];
%  J(5,:) = [1,2,3,4,6,0];
%  J(n,:) = [n-5,n-4,n-3,n-2,n-1,0];

%  for i=6:n-1
%    J(i,:) = [i-5,i-4,i-3,i-2,i-1,i+1];
%  end

%  for i=1:n
%      somme = 0;
%   for j=1:6
%      if (J(i,j) ~= 0)
%          somme = somme + x(J(i,j))*(1+x(J(i,j)));
%
%      end
%   end
%   f(i) = x(i)*(2+5*x(i)^2) + 1 - somme;
%
%  end
%
%  F = sum(f.^2);

f(1) = x(1)*(2+5*x(1)^2) + 1 - x(2)*(1+x(2));
f(2) = x(2)*(2+5*x(2)^2) + 1 - x(1)*(1+x(1))-x(3)*(1+x(3));
f(3) = x(3)*(2+5*x(3)^2) + 1 - sum(x(1:2).*(ones(2,1)+x(1:2)))-x(4)*(1+x(4));
f(4) = x(4)*(2+5*x(4)^2) + 1 - sum(x(1:3).*(ones(3,1)+x(1:3)))-x(5)*(1+x(5));
f(5) = x(5)*(2+5*x(5)^2) + 1 - sum(x(1:4).*(ones(4,1)+x(1:4)))-x(6)*(1+x(6));
f(n) = x(n)*(2+5*x(n)^2) + 1 - sum(x(n-5:n-1).*(ones(5,1)+x(n-5:n-1)));

for i=6:n-1
f(i) = x(i)*(2+5*x(i)^2) + 1 - sum(x(i-5:i-1).*(ones(5,1)+x(i-5:i-1)))-x(i+1)*(1+x(i+1));
end

F = sum(f.^2);

%  for i=1:n
%    for k=1:6
%      if (J(i,k)~= 0)
%         gradf2(J(i,k),i) = 2*x(J(i,k)) + 1;
%      end
%    end
%  end

for i=6:n-1
end

%  H = cell(n,1);

%  for i=1:n
%     H{i} = zeros(n,n);
%     H{i}(i,i) = 16*x(i);
%     for k=1:6
%        if (J(i,k) ~= 0)
%          H{i}(J(i,k),J(i,k))=2;
%        end
%     end
%  end

%  Hessien = zeros(n,n);
%  for i=1:n
%  end

H1(1,1) = 16*x(1)*f(1)+sum(f(2:6));
H1(n-4,n-4) = 16*x(n-4)*f(n-4) + 2*sum([f(n-5),f(n-3:n)]);
H1(n-3,n-3) = 16*x(n-3)*f(n-3) + 2*sum([f(n-4),f(n-2:n)]);
H1(n-2,n-2) = 16*x(n-2)*f(n-2) + 2*sum([f(n-3),f(n-1:n)]);
H1(n-1,n-1) = 16*x(n-1)*f(n-1) + 2*sum([f(n-2),f(n)]);
H1(n,n)     = 16*x(n)*f(n) + 2*f(n-1);

for i=2:n-5
H1(i,i) = 16*x(i)*f(i) + 2*sum([f(i-1),f(i+1:i+5)]);
end

H1 = sparse(H1);

for i=2:n
end

H2 = sparse(H2);

Hessien = 2*(H1 + H2); Hessien = sparse(Hessien);

end

if prob == 2

const = 100;

for i=2:n
f(i) = const*(x(i)-x(i-1)^2)^2 + (1-x(i))^2;
end

F = sum(f);

G(1) = 4*const*x(1)^3-4*const*x(1)*x(2);
G(n) = 2*(const+1)*x(n) -2*const*x(n-1)^2 -2;

for i=2:n-1
G(i) = 2*(const+1)*x(i) - 2*const*x(i-1)^2-2+4*const*x(i)^3-4*const*x(i)*x(i+1);
end

G = G';

Hessien(1,1) = 12*const*x(1)^2-4*const*x(2);
Hessien(n,n) = 2*const+2;
Hessien(n,n-1) = -4*const*x(n-1);
Hessien(n-1,n) = -4*const*x(n-1);

Hessien = sparse(Hessien);

for i=2:n-1
Hessien(i,i-1) = -4*const*x(i-1);
Hessien(i-1,i) = -4*const*x(i-1);
Hessien(i,i) = 2*(1+const)+12*const*x(i)^2-4*const*x(i+1);
end

Hessien = sparse(Hessien);

end

if prob == 21

const = 1000;

for i=2:n
f(i) = const*(x(i)-x(i-1)^2)^2 + (1-x(i))^2;
end

F = sum(f);

G(1) = 4*const*x(1)^3-4*const*x(1)*x(2);
G(n) = 2*(const+1)*x(n) -2*const*x(n-1)^2 -2;

for i=2:n-1
G(i) = 2*(const+1)*x(i) - 2*const*x(i-1)^2-2+4*const*x(i)^3-4*const*x(i)*x(i+1);
end

G = G';

Hessien(1,1) = 12*const*x(1)^2-4*const*x(2);
Hessien(n,n) = 2*const+2;
Hessien(n,n-1) = -4*const*x(n-1);
Hessien(n-1,n) = -4*const*x(n-1);

Hessien = sparse(Hessien);

for i=2:n-1
Hessien(i,i-1) = -4*const*x(i-1);
Hessien(i-1,i) = -4*const*x(i-1);
Hessien(i,i) = 2*(1+const)+12*const*x(i)^2-4*const*x(i+1);
end

Hessien = sparse(Hessien);

end

if prob == 22

for i=2:n
f(i) = 100*(x(i)-x(i-1))^2 + (1-x(i))^2;
end

F = sum(f);

%  for i=2:n
%  end

%  G = full(G);

G(1) = -200*x(2) +200*x(1);
G(n) = 202*x(n) -200*x(n-1) -2;

for i=2:n-1
G(i) = 402*x(i) - 200*x(i-1) -200*x(i+1) -2;
end

G = G';

%  H = cell(n,1);

%  for i=2:n
%     H{i} = zeros(n,n);
%     H{i}(i,i) = 202;
%     H{i}(i-1,i-1) = 200;
%     H{i}(i,i-1) = -200;
%     H{i}(i-1,i) = -200;
%  end

%  Hessien = zeros(n,n);
%  for i=2:n
%     Hessien = Hessien + H{i};
%  end

%  Hessien = sparse(Hessien);

Hessien(1,1) = 200;
Hessien(n,n) = 202;
Hessien(n-1,n) = -200;
Hessien(n,n-1) = -200;

Hessien = sparse(Hessien);

for i=2:n-1
Hessien(i,i-1) = -200;
Hessien(i-1,i) = -200;
Hessien(i,i) = 402;
end

Hessien = sparse(Hessien);

end

if prob == 23

const = 100;

for i=2:n
f(i) = const*(x(i)-x(i-1)^2)^2 + (1-x(i-1))^2;
end

F = sum(f);

G(1) = 4*const*x(1)^3-4*const*x(1)*x(2)-2+2*x(1);
G(n) = 2*const*x(n) -2*const*x(n-1)^2;

if n>2
for i=2:n-1
G(i) = 2*const*x(i) - 2*const*x(i-1)^2-4*const*x(i)*x(i+1)+4*const*x(i)^3-2+2*x(i);
end
end

G = G';

Hessien(1,1) = 12*const*x(1)^2-4*const*x(2)+2;
Hessien(n,n) = 2*const;
Hessien(n,n-1) = -4*const*x(n-1);
Hessien(n-1,n) = -4*const*x(n-1);

Hessien = sparse(Hessien);

if n>2
for i=2:n-1
Hessien(i,i-1) = -4*const*x(i-1);
Hessien(i-1,i) = -4*const*x(i-1);
Hessien(i,i) = 2*const+12*const*x(i)^2-4*const*x(i+1)+2;
end
end

Hessien = sparse(Hessien);

end

if prob == 3

for i=1:n
for j=1:n
f(i,j) = sin(x(i))^2*sin(x(j))^2*sin(x(i)-x(j))^2;
end
end

F = -sum(sum(f));

for i=1:n
G(i) = -2*sum((sin(x)).^2.*(sin(2*x(i))*(sin(x(i)*ones(n,1)-x)).^2 +...
(sin(x(i)))^2*sin(2*(x(i)*ones(n,1)-x))));
end

G = G';

for i=1:n
Hessien(i,i) = -2*sum(sin(x).^2.*(2*cos(2*x(i))*sin(x(i)*ones(n,1)-x).^2 ...
+2*sin(2*x(i))*sin(2*(x(i)*ones(n,1)-x))+ ...
2*sin(x(i))^2*cos(2*(x(i)*ones(n,1)-x)))) +...
4*sin(x(i))^4;
for j=i+1:n
Hessien(i,j) = -2*(sin(2*x(j))*(sin(2*x(i))*(sin(x(i)-x(j)))^2 +...
(sin(x(i)))^2*sin(2*(x(i)-x(j)))) + ...
(sin(x(j)))^2*(-sin(2*x(i))*sin(2*(x(i)-x(j)))- ...
2*(sin(x(i)))^2*cos(2*(x(i)-x(j)))));
Hessien(j,i) = Hessien(i,j);
end
end

end

if prob == 4

if ((n/4 ~= floor(n/4)) & (n/4 ~= ceil(n/4)))
disp('Error: n has to be a multiple of 4');
F = 999999;
return;
end

div = n/4;

for k=1:div
f(4*k-3) = x(4*k -3) + 10*x(4*k -2);
f(4*k-2) = sqrt(5)*(x(4*k -1)-x(4*k));
f(4*k-1) = (x(4*k-2)-2*x(4*k-1))^2;
f(4*k)   = sqrt(10)*(x(4*k-3) -x(4*k))^2;
end

F = sum(f.^2);

for k=1:div
end

for k=1:div
H1(4*k-2,4*k-2) = 2*f(4*k-1);
H1(4*k-2,4*k-1) = -4*f(4*k-1);
H1(4*k-1,4*k-2) = -4*f(4*k-1);
H1(4*k-1,4*k-1) = 8*f(4*k-1);
H1(4*k-3,4*k-3) = 2*sqrt(10)*f(4*k);
H1(4*k,4*k)     = 2*sqrt(10)*f(4*k);
H1(4*k-3,4*k)   = -2*sqrt(10)*f(4*k);
H1(4*k,4*k-3)   = -2*sqrt(10)*f(4*k);
end

H1 = sparse(H1);

for i=2:n
end

H2 = sparse(H2);

Hessien = 2*(H1 + H2); Hessien = sparse(Hessien);

end

if prob == 5

for i=1:n
f(i) = 100*(x(1)-x(i)^2)^2+(1-x(i))^2;
end

F = sum(f);

G(1) = 200*sum(x(1)*ones(n-1,1) - x(2:n));
for i=2:n
G(i) = -400*(x(1)-x(i)^2)*x(i)-2*(1-x(i));
end

G = G';

Hessien(1,1) = 200*n;
for i=2:n
Hessien(i,1) = -400*x(i);
Hessien(1,i) = -400*x(i);
Hessien(i,i) = -400*(x(1)-x(i)^2)+800*x(i)^2+2;
end

Hessien = sparse(Hessien);

end

if prob == 6

for i=2:n
f(i) = sqrt(i)*(2*x(i)-x(i-1));
end

F = sum(f.^2);

G(1) = -4*(2*x(2)-x(1));
G(n) = 4*n*(2*x(n)-x(n-1));
for i=2:n-1
G(i) = 4*i*(2*x(i)-x(i-1)) - 2*(i+1)*(2*x(i+1)-x(i));
end

G = G';

Hessien(1,1) = 4;
Hessien(n,n) = 8*n;
Hessien(n-1,n) = -4*n;
Hessien(n,n-1) = -4*n;
for i=2:n-1
Hessien(i-1,i) = -4*i;
Hessien(i,i-1) = -4*i;
Hessien(i,i)   = 8*i+2*(i+1);
end

Hessien = sparse(Hessien);

end

if prob == 7

if ((n/2 ~= floor(n/2)) & (n/2 ~= ceil(n/2)))
disp('Error: n has to be a multiple of 2');
F = 999999;
return;
end

div = n/2;

for k=1:div
f(k) = 100*(x(2*k)-x(2*k-1)^2)^2+(1-x(2*k-1))^2;
end

F = sum(f);

for k=1:div
G(2*k-1) = -400*x(2*k-1)*(x(2*k)-x(2*k-1)^2)-2*(1-x(2*k-1));
G(2*k) = 200*(x(2*k)-x(2*k-1)^2);
end

G = G';

for k=1:div
Hessien(2*k-1,2*k-1) = 1200*x(2*k-1)^2-400*x(2*k) + 2;
Hessien(2*k-1,2*k) = -400*x(2*k-1);
Hessien(2*k,2*k-1) = -400*x(2*k-1);
Hessien(2*k,2*k)   = 200;
end

Hessien = sparse(Hessien);

end

if prob == 8

if (n<3)
disp('Error: n has to be greater or equal to 3');
F = 999999;
return;
end

f(1) = 2*x(1) -x(2) + (x(1) + (n+2)/(n+1))^3/(2*(n+1)^2);
f(n) = 2*x(n) -x(n-1) + (x(n) + (2*n+1)/(n+1))^3/(2*(n+1)^2);

for i=2:n-1
f(i) = 2*x(i) -x(i-1) -x(i+1) +(x(i) + (n+i+1)/(n+1))^3/(2*(n+1)^2);
end

F = sum(f.^2);

for i=2:n-1
end

for i=1:n
H1(i,i) = f(i)*(x(i)+(n+i+1)/(n+1));
end
H1 = 3/(n+1)^2*H1;

H1 = sparse(H1);

for i=2:n
end

H2 = sparse(H2);

Hessien = 2*(H1 + H2); Hessien = sparse(Hessien);

end

if prob == 9

f(1) = (3-2*x(1))*x(1) -2*x(2) +1;
f(n) = (3-2*x(n))*x(n) -x(n-1) +1;

for i=2:n-1
f(i) = (3-2*x(i))*x(i) -x(i-1) -2*x(i+1) +1;
end

F = sum(f.^2);

for i=2:n-1
end

for i=1:n
H1(i,i) = f(i);
end
H1 = -4*H1;

H1 = sparse(H1);

for i=2:n
end

H2 = sparse(H2);

Hessien = 2*(H1 + H2); Hessien = sparse(Hessien);

end

if prob == 10

if (n~=31)
disp('Error: n has to be equal to 31');
F = 999999;
return;
end

index = [1:n-1]';

f(30) = x(1);
f(31) = x(2)-x(1)^2-1;
for i=1:29
f(i) = sum(index.*x(2:n).*((i/29*ones(n-1,1)).^(index-ones(n-1,1))) )-...
(sum(x.*((i/29*ones(n,1)).^([0:n-1]'))))^2 - 1;
end

F = sum(f.^2);

%  for i=1:29
%    gradf(1,i) = -2*x(1) - 2*sum(x.*((i/29*ones(n,1)).^([0:n-1]'))) + 2*x(1);
%    for j=2:31
%      gradf(j,i) = (i-1)*i/29^(i-2) - ...
%                   2*x(i)*((i/29)^(i-1))^2 - ...
%                   2*sum((i/29)^(i-1)*x.*((i/29*ones(n,1)).^([0:n-1]'))) +...
%                   2*x(i)*((i/29)^2)^(i-1);
%    end
%  end
%

for i=1:29
g1 =0;
for j=2:n
g1 = g1+x(j)*(i/29)^(j-1);
end
for k=2:n
g2 = 0;
for h=1:n
g2 = g2+x(h)*(i/29)^(k-1)*(i/29)^(h-1);
end
g2 = g2 - x(k)*(i/29)^(k-1)*(i/29)^(k-1);
end

end

for i=1:n
H1(i,i) = -2*((1/29)^(i-1))^2;
for j=i+1:n
H1(i,j) = -2*(1/29)^(i+j-2);
H1(j,i) = H1(i,j);
end
end
H1 = f(1)*H1;

for k=2:29

for i=1:n
Htemp(i,i) = -2*((k/29)^(i-1))^2;
for j=i+1:n
Htemp(i,j) = -2*(k/29)^(i+j-2);
Htemp(j,i) = Htemp(i,j);
end
end

Htemp = f(k)*Htemp;
H1 = H1+Htemp;
end

H1(1,1) = H1(1,1)-2*f(31);

for i=2:n
end

Hessien = 2*(H1+H2);

end

if prob == 11

for i=1:n-1
f(i) = x(i) + sum(x) - (n+1);
end
f(n) = prod(x) -1;

F = sum(f.^2);

for i=1:n-1
end

for i=2:n-1
d = [1:i-1,i+1:n];
end

for i=1:n-1
for j=i+1:n
d = [1:i-1,i+1:j-1,j+1:n];
H1(i,j) = prod(x(d));
H1(j,i) = H1(i,j);
end
end

H1 = f(n)*H1;

for i=2:n
end

Hessien = 2*(H1+H2);

end

if prob == 12

m = 2*n;

for i=1:n
f(i) = x(i) -2/m*sum(x) -1;
end

for i=n+1:m
f(i) = -2/m*sum(x) -1;
end

F = sum(f.^2);

for i=1:n
end

for i=n+1:m
end

for i=2:m
end

end

if prob == 13

m=2*n;

for i=1:m
f(i) = i*sum(([1:n]').*x)-1;
end

F = sum(f.^2);

for i=1:m
end

for i=2:m
end

end

if prob == 14

m = 2*n;

f(1) = -1;
f(m) = -1;

for i=2:m-1
f(i) = (i-1)*sum(([2:n-1]').*x(2:n-1)) -1;
end

F = sum(f.^2);

for i=2:m-1
end

for i=2:m
end

end

if prob == 15

f(n) = (x(n) - 1)^2;

for i=1:n-1
f(i) = (x(i)^2-1)^2;
end

F = sum(f);

G(n) = 2*(x(n)-1);

for i=1:n-1
G(i) = 4*x(i)*(x(i)^2-1);
end

G = G';

Hessien(n,n) = 2;

for i=1:n-1
Hessien(i,i) = 12*x(i)^2-4;
end

Hessien = sparse(Hessien);

end

if prob == 16

if ((n/2 ~= floor(n/2)) & (n/2 ~= ceil(n/2)))
disp('Error: n has to be a multiple of 2');
F = 999999;
return;
end

F = 0;
a = 0.4238537991;

for i=1:n/2
F = F+(x(2*i-1)+a)^2*(x(2*i)+a)^2+(x(2*i-1)+a)^2+(x(2*i)+a)^2;
end

for i=1:n/2
G(2*i-1) = 2*(x(2*i-1)+a)*((x(2*i)+a)^2+1);
G(2*i  ) = 2*(x(2*i)+a)*((x(2*i-1)+a)^2+1);
end

G = G';

for i=1:n/2
Hessien(2*i-1,2*i-1) = 2*((x(2*i)+a)^2+1);
Hessien(2*i-1,2*i)   = 4*(x(2*i-1)+a)*(x(2*i)+a);
Hessien(2*i,2*i-1)   = 4*(x(2*i-1)+a)*(x(2*i)+a);
Hessien(2*i,2*i)     = 2*((x(2*i-1)+a)^2+1);
end

Hessien = sparse(Hessien);

end

if prob == 17
if ((n/2 ~= floor(n/2)) & (n/2 ~= ceil(n/2)))
disp('Error: n has to be a multiple of 2');
F = 999999;
return;
end

F=0;

for i=1:n/2
F = F + (exp(x(2*i-1)+x(2*i))-1)^2;
G(2*i-1) = 2*(exp(x(2*i-1)+x(2*i))-1)*exp(x(2*i-1)+x(2*i));
G(2*i) = G(2*i-1);
Hessien(2*i-1,2*i-1) = 2*exp(x(2*i-1)+x(2*i))*(2*exp(x(2*i-1)+x(2*i))-1);
Hessien(2*i,2*i) = Hessien(2*i-1,2*i-1);
Hessien(2*i,2*i-1) = Hessien(2*i-1,2*i-1);
Hessien(2*i-1,2*i) = Hessien(2*i-1,2*i-1);
end

G = G';
Hessien = sparse(Hessien);

endrun.m0000644001023000102300000000033107043371170012646 0ustar  hwolkowihwolkowiczclear global A
clear global D
clear global counter
clear all
pack
global A D counter
sparse A D
n=10,prob=23;
x0=80*ones(n,1);R=3;
[x_min,f_min,fcteval,compteur,iteration,temps,comptons,drapeau] = trmgould(x0,R,prob)
runtrs.m0000644001023000102300000000032111350447103013372 0ustar  hwolkowihwolkowiczclear all
profile on
global A D counter
n=100000;density=.0001;
A=sprandsym(n,density);
D=speye(n);
a=randn(n,1);
s=.1;
dgaptol=1e-12;
[xopt,fopt,lopt,lambdamin,flag] = newtrust4b(a,s,dgaptol);
profile report
sparsetrust.m0000644001023000102300000001345207046633705014461 0ustar  hwolkowihwolkowiczfunction [xc,histout,costdata] = sparsetrust(x0,f,parms,resolution)
%
%
% function [x,histout] = sparsetrust(x0,f,parms,resolution)
%
% RW Sparse-Lanzcos-Trust region algorithm
%
% Input: x0 = initial iterate
%        f = objective function,
%            the calling sequence for f should be
%            [fout,gout]=f(x) where fout=f(x) is a scalar
%              and gout = grad f(x) is a COLUMN vector
%        parms = [tol, eta, maxitnl, maxitl]
%        tol = termination criterion norm(grad) < tol
%        eta = forcing term in linear iteration (optional) default = .1
%        maxitnl = maximum nonlinear iterations (optional) default = 100
%        maxitl = maximum linear iterations (optional) default = 20
%        resolution = estimated accuracy in functions/gradients (optional)
%                     default = 1.d-12
%                     The finite difference increment in the difference
%                     Hessian is set to sqrt(resolution).
%
%
%
% Output: x = solution
%         histout = iteration history
%             Each row of histout is
%         costdata = [num f, num grad]
%
% Requires: dirdero.m
%
% trust region algorithm parameters
%
mu0=.25; mulow=.25; muhigh=.75;
%
% set maxit, the resolution, and the difference increment
%
debug=0;
tol=parms(1); np=length(parms); itc=1;
if np < 2; eta=.1; else eta = parms(2); end
if np < 3; maxit=100; else maxit = parms(3); end
if np < 4; maxitl=20; else maxitl = parms(4); end
if nargin < 5
resolution = 1.d-12;
end
numf=0; numg=0;
hdiff=sqrt(resolution);
t=100; itc=0; xc=x0; n=length(x0);
[fc,gc]=feval(f,xc); numf=1; numg=1;
%
paramstr=[eta,maxitl]; trcount=1;
while(norm(gc) > tol & itc <=maxit) & trcount < 30
itc=itc+1;
vd=size(dirs); itsl=vd(2); numf=numf+itsl; numg=numg+itsl;
xt=xc+s; ft=feval(f,xt); numf=numf+1; ared=ft-fc;
w= dirdero(xc, s, f, gc); numg=numg+1;
pred=gc'*s + .5*(w'*s); rat=ared/pred;
if rat > mu0
xc=xc+s; [fc,gc]=feval(f,xc); numf=numf+1;
else
for k=1:itsl dsums(k)=norm(sum(dirs(:,1:k),2)); end
trcount=1;
while rat <= mu0 & trcount < 30
xt=xc+s; ft=feval(f,xt); numf=numf+1; ared=ft-fc;
%
%      Only compute a new pred if ared < 0 and there's some hope
%      that rat > mu0
%
if ared < 0
w= dirdero(xc, s, f, gc); numg=numg+1; pred=gc'*s + .5*(w'*s);
end
rat=ared/pred;
itsl=kout; trcount=trcount+1;
if trcount > 30
%           histout=ithist(1:itc+1,:); costdata=[numf,numg];
disp(' stagnation in CGTR')
end
end
if trcount < 30
xc=xt; [fc,gc]=feval(f,xc); numf=numf+1; numg=numg+1;
end
end
end
histout=ithist(1:itc+1,:); costdata=[numf,numg];
%
% find the point of intersetion of the TR boundary and the PL path
%
st=dirs(:,1); inside=1;
if norm(st) > trrad | itsl == 1
kout=1;
else
for k=2:itsl
if norm(st+dirs(:,k)) > trrad  & inside == 1
kout=k;
alpha=(-bc + sqrt(bc*bc - 4*ac*cc))/(2*ac); st=st+alpha*p;
inside=0;
else
st=st+dirs(:,k);
end
end
end
%
%
%
function [x, directions]  = trcg(xc, f, gc, delta, params, pcv)
%
% Solve the trust region problem with preconditioned conjugate-gradient
%
% C. T. Kelley, January 13, 1997
%
% This code comes with no guarantee or warranty of any kind.
% function [x, directions]
%                    = trcg(xc, f, gc, delta)
%
%
%
% Input:        xc=current point
%               b=right hand side
%           f = function, the calling sequence is
%                gc has usually been computed
%                before the call to dirdero
%           params = two dimensional vector to control iteration
%                params(1) = relative residual reduction factor
%                params(2) = max number of iterations
%           pcv, a routine to apply the preconditioner
%                if omitted, the identity is used.
%                The format for pcv is
%                       function px = pcv(x).
%
% Output:   x = trial step
%           directions = array of search directions TR radius reduction
%
%

%
% initialization
%
n=length(xc); errtol = params(1); maxiters = params(2);
x=zeros(n,1); b=-gc; r=b - dirdero(xc, x, f, gc);
if nargin == 5
z=r;
else
z = feval(pcv, r);
end
rho=z'*r;
tst=norm(r);
terminate=errtol*norm(b);
it=1;
directions=zeros(n,1);
hatdel=delta*(1-1.d-6);
while((tst > terminate) & (it <= maxiters) & norm(x) <= hatdel)
%
%
%
if(it==1)
p = z;
else
beta=rho/rhoold;
p = z + beta*p;
%
% end if
%
end
w = dirdero(xc, p, f, gc);
alpha=p'*w;
%
% If alpha <=0 head to the TR boundary and return
%
ineg=0;
if(alpha <= 0)
ac=p'*p; bc=2*(x'*p); cc=x'*x - delta*delta;
alpha=(-bc + sqrt(bc*bc - 4*ac*cc))/(2*ac);
disp(' negative curvature')
else
alpha=rho/alpha;
if norm(x+alpha*p) > delta
ac=p'*p; bc=2*(x'*p); cc=x'*x - delta*delta;
alpha=(-bc + sqrt(bc*bc - 4*ac*cc))/(2*ac);
end
end
x=x+alpha*p;
directions(:,it)=alpha*p;
r = r - alpha*w;
tst=norm(r);
rhoold=rho;
if nargin < 6 z=r; else z = feval(pcv, r); end
rho=z'*r;
it=it+1;
%
% end while
%
end
startup.m0000644001023000102300000000010707044373556013560 0ustar  hwolkowihwolkowiczpath(path,'/fsys4/u/hwolkowicz/priv.d/software.d/matlab.d/sparfun.d');
trmgould.m0000644001023000102300000001306207043360300013675 0ustar  hwolkowihwolkowiczfunction [x_min,f_min,fcteval,compteur,iteration,temps,comptons,drapeau] = trmgould(x0,R,prob)

% Call: [x_min,f_min,fcteval,compteur,iteration,temps,comptons,drapeau] = trmgould(x0,R,prob)
% This is the trust region method found in Gould et al. 
% This algorithm is used to minimize an unconstrained function. The function's
% objective should be available using the file 'objective.m' and the function's
% function uses the subroutine 'newtrust.m' to solve the trust region
% subproblem
%
% INPUT VARIABLES:
%
% x0 ... initial estimate of a minimizer of the function
% R  ... initial radius of the trust region
% prob ... number of the problem to be solve (used to access the files
%
% OUTPUT VARIABLES:
%
% x_min   ... an approximate of a minimizer of the function
% f_min   ... objective at x_min. This approximates the minimum of the function
% fcteval ... number of fonctions evaluations
% compteur .. number of matrix-vector multiplications
% iteration . number of iterations done by the trust region method
% temps   ... time taken to find the approximate solution
% comptons .. relative distance of the optimal Lagrange multipliers of the
%             trust region subproblems to the minimum eigenvalue of the
%             Hessians
% drapeau ... vector of integers (0,1,2) used to distinguish the different
%             cases for the trust region subproblemes

% GLOBAL VARIABLES
global A D counter % global variables used by the newtrust4b.m subroutine

global A counter x

b = ' '; % define blanks
n = size(x0,1);

temps = cputime;

% Constant used to update the trust region radius (i.e. evaluate the
% performance of the trust region step
n1 = 0.01;
n2 = 0.95;
gam1 = 0.5;
gam2 = 2;
s = R;

x = x0; % x is the current approximation to a minimizer of the function

[fopt,g,A] = objgradhess(x,prob);; % fopt is the objective evaluated at x (i.e.
% the current approximation to the minimum
% of the function, g is the gradient and A
% is the Hessian.

% initiaze variables
fcteval = 1;
iteration = 1;
compteur = 0;

disp(['iteration',b,b,b,b,'time elapsed so far',b,b,b,b,'f(x)',b,b,b,b,b,b,b,b,'norm of the gradient']);
disp([b,b,b,b,num2str(iteration),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(cputime-temps),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(fopt),b,b,b,b,b,b,b,b,b,b,b,num2str(norm(g)),b,b,b,b]);

if (norm(g)/(1+abs(fopt)) < 0.00001) % relative error on the norm of the

% Initial approximate is good enough
x_min = x;
f_min = fopt;
return; % END OF ALGORITHM
end

while (norm(g)/(abs(fopt)+1) > 0.00001)

epsilon = 10^(-6);
dgaptol = min(epsilon,10^(-4)*norm(g)); % relative duality gap. This is used
% so that initially we don't ask a lot
% of precision for the trust region
% subproblems, but as we get close to
% the solution (gradient is small), we
% do.
% SOLVE THE TRUST REGION SUBPROBLEM
[d,qopt,lopt,lambdamin,iter,drapeau(iteration)] = newtrust(-g,s,dgaptol);

% d is the step direction
% qopt is the predicted value of the objective at x+d (predicted by the
% lopt is the Lagranga multiplier of the optimal solution to the trust
% region subproblem.
% lambdamin is the minimum eigenvalue of the current matrix A (the Hessian at
% x).
% drapeau(iteration) is an integer that tells us which case occurred in the
% trust region subproblem (hard case, case 1 and 2 and easy case).

comptons(iteration) = (lambdamin-lopt)/(abs(lambdamin) + 1); % relative
% distance of the optimal Lagrange multiplier to the
% minimum eigenvalue of A (Hessian at x). This is just
% used for the purpose of analysing the iterates.;

compteur = compteur + counter; % counter is the number of matrix-vector
% multiplication done in the last trust region
% subproblem. So update compteur, the total
% number of matrix-vector multiplications

f_predicted = fopt + 1/2*qopt;
f_real = objective(x+d,prob);

r = (fopt - f_real)/(fopt - f_predicted); % r is a measure of our improvement
% r=1 is good, r=0 is bad
if (r > n1)
% improvement of the objective function is good enough so take a step
x = x+d;
fcteval = fcteval + 1;
end

if (r > n2)
% improvement of the objective function is great, increase the trust region
s = gam2*s;
end

if (r < n1)
% improvement of the objective function is not satisfactory. Reduce the
s = gam1*s;
end

disp('new values');
r
s

iteration = iteration + 1; % one more iteration has been done

% display info on the iteration
disp(['iteration',b,b,b,b,'time elapsed so far',b,b,b,b,'f(x)',b,b,b,b,b,b,b,b,'norm of the gradient']);
disp([b,b,b,b,num2str(iteration),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(cputime-temps),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(fopt),b,b,b,b,b,b,b,b,b,b,b,num2str(norm(g)),b,b,b,b]);

end

% return the approximation to the optimal solutions
f_min = fopt;
x_min = x;

temps = cputime - temps; % temps = time to solve the problem

```