%J.C. Spall, Jan. 2000 %This code implements SPSA with constraints for theta to lie in %a specified hypercube. %Code includes more "bells and whistles" than required for basic SPSA (as %available at the SPSA web site and code 'spsa_basic_constrained'). %Code includes loss-based and theta-based %blocking of iterates (both of these blocking steps are optional). %Code also allows for averaging SP gradients at each iteration. This code includes %some evaluation material that is not essential to a basic run of SPSA (i.e., when %cases > 1, the code is performing multiple runs on the same problem). clear all global p z sigma; %declaration of random var. (z) used for normal noise %generation in loss fn. calls given seed above; %also sigma in noise (noise may be dependent on theta) p=10; n=1000; %total no. of loss measurements cases=1; alpha =.602; gamma =.101; a=.0017; c=1.9; %chosen by standard guidelines gavg=1; %no. of grad. estimates averaged in update A=50; lossfinalsq=0; %variable for cum.(over 'cases')squared loss values lossfinal=0; %variable for cum. loss values theta_lo=-1000*ones(p,1); %lower bounds on theta theta_hi=1000*ones(p,1); %upper bounds on theta tolerancetheta=Inf; %max. allowable change in any element of theta toleranceloss=0; %tolerance in loss-based blocking step avg=0; %no. of loss evaluations used in loss-based %blocking (0 implies no blocking) theta_0=6*ones(p,1); lossfinaleval='loss_test'; % choice of loss function for final perf. % evaluation loss='loss_test'; % loss function used in algorithm operations rand('seed',31415927) randn('seed',3111113) sigma=.1; %multiplier in loss noise (lower bd. to standard dev.) feval(lossfinaleval,theta_0) for i=1:cases theta=theta_0; lossold=0; %lossold initialization is for use in loss-based blocking step below for j=1:avg lossold=lossold+feval(loss,theta); end for k=0:(n-avg)/(2*gavg+avg)-1 ak = a/(k+1+A)^alpha; ck = c/(k+1)^gamma; ghat=0; for j=1:gavg delta = 2*round(rand(p,1))-1; thetaplus = theta + ck*delta; thetaminus = theta - ck*delta; yplus=feval(loss,thetaplus); yminus=feval(loss,thetaminus); ghat = (yplus - yminus)./(2*ck*delta)+ghat; end thetalag=theta; theta=theta-ak*ghat/gavg; % Two lines below invoke constraints theta=min(theta,theta_hi); theta=max(theta,theta_lo); % Steps below perform "blocking" step with "avg" no. of loss evaluations % This blocking is based on extra loss evaluation(s) lossnew=0; for j=1:avg lossnew=lossnew+feval(loss,theta); end if lossnew > lossold-avg*toleranceloss; %if avg=0, this statement is always false theta=thetalag; else %statements to follow are harmless when avg=0 lossold1=lossold; lossold=lossnew; end % Blocking below is based on magnitude of theta change (no loss evals.) if max(abs(thetalag-theta)) > tolerancetheta; theta=thetalag; lossold=lossold1; %only relevant if also using loss-based blocking end end lossvalue=feval(lossfinaleval,theta); lossfinalsq=lossfinalsq+lossvalue^2; lossfinal=lossfinal+lossvalue; end % Display results: Mean loss value and standard deviation % disp('mean loss value over "cases" runs') lossfinal/cases % if cases > 1 disp('sample standard deviation of mean loss value') sd=((cases/(cases-1))^.5)*(lossfinalsq/cases-(lossfinal/cases)^2)^.5; sd=sd/(cases^.5) else end