%Some aspects of code tailored to Example 14.9 %on CRN vs. Non-CRN. This code implements CRN version %of SPSA. global p z lambda; %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) format long; p=10; N=20000; %total no. of loss measurements (2 x no. of iterations) cases=50; alpha =1; gamma =.49; a=.7; c=.5; %chosen by standard guidelines gavg=1; %no. of grad. estimates averaged in update A=0; lambda=[1.1025,1.6945,1.4789,1.9262,0.7505,1.3267,0.8428,0.7247,0.7693,1.3986]'; theta_star=[0.286,0.229,0.248,0.211,0.325,0.263,0.315,0.327,0.323,0.256]'; lossfinalsq=0; %variable for cum.(over 'cases')squared loss values lossfinal=0; %variable for cum. loss values theta_lo=zeros(10,1); %lower bounds on theta theta_hi=Inf; %upper bounds on theta tolerancetheta=Inf; %max. allowable change in any element of theta toleranceloss=0; %tolerance in loss-based blocking step theta_0=ones(p,1); loss='loss_exponential_CRN'; % loss function used in algorithm operations rand('seed',35111137) losses=zeros(cases,1); thetanorm=zeros(cases,1); for i=1:cases theta=theta_0; for k=1:N/(2*gavg) ak = a/(k+A)^alpha; ck = c/k^gamma; ghat=0; for j=1:gavg delta = 2*round(rand(p,1))-1; thetaplus = theta + ck*delta; thetaminus = theta - ck*delta; % Two lines below invoke constraints thetaplus=min(thetaplus,theta_hi); thetaminus=max(thetaminus,theta_lo); yplusminus=feval(loss,thetaplus,thetaminus); ghat = yplusminus./(2*ck*delta)+ghat; end theta=theta-ak*ghat/gavg; % Two lines below invoke constraints theta=min(theta,theta_hi); theta=max(theta,theta_lo); end % statements below produce loss value for final evaluation lossvalue=0; for j=1:p lossvalue=lossvalue+lambda(j)/(lambda(j)+theta(j)); end lossvalue=lossvalue+theta'*theta; losses(i)=lossvalue; %used in Example 14.9 lossfinalsq=lossfinalsq+lossvalue^2; lossfinal=lossfinal+lossvalue; thetanorm(i)=((theta-theta_star)'*(theta-theta_star))^0.5; 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 losses %statement for use with example 14.9 thetanorm/((theta_0-theta_star)'*(theta_0-theta_star))^0.5 %for use with example 14.9