function [xopt,fopt] = r1nes(f,d,x,timeout) % Written by Sun Yi (yi@idsia.ch). % parameters etax = 1; % learning rate for the mean etaa = 0.1; % learning rate for the scale factor etac = 0.1; % primary representation is L = floor(8*log(d)); % population size shape = max(0.0, log(L/2+1.0)-log(1:L)); % fitness shaping shape = shape / sum(shape); % initialize normalize = @(v) v/sqrt(v'*v); xopt = x; fopt = f(x); a = 0; % fixed diagonal strength c = 0; v = normalize(randn(d,1)); r = exp(c); u = r*v; weights = zeros(1,L); fit = zeros(1,L); tm = cputime; while cputime - tm < timeout % step 1: sampling W = randn(d,L) + u*randn(1,L); X = repmat(x,1,L) + exp(a)*W; for i = 1 : L, fit(i) = f(X(:,i)); end % step 2: fitness reshaping [~, idx] = sort(fit); weights(idx) = shape; if fit(idx(1)) < fopt xopt = X(:,idx(1)); fopt = fit(idx(1)); end % step 3: compute the gradient dx = exp(a)*sum(repmat(weights,d,1).*W,2); v = u/r; W = W(:,weights ~= 0); weights = weights(weights~=0); wws = sum(W.^2,1); wvs = v'*W; kp1 = ((r^2-d+2)*(wvs.^2) - (r^2+1)*wws)/(2*r*(d-1)); kp2 = wvs/r; da = 1/(2*(d-1))*((wws-d) - (wvs.^2-1))*weights'; du = (kp1*weights')*v + W*(kp2.*weights)'; dc = (du'*v)/r; dv = du/r - dc*v; % step 4: compute the update x = x + etax * dx; a = a + etaa * da; r = sqrt(u'*u); eps = min(etac,2*sqrt(r^2/(du'*du))); % prevent the update from getting crazy if dc > 0 % additive u = u + eps*du; else % multiplicative r = r * exp(eps*dc); v = normalize(u/r + eps*dv); u = r*v; end if a < -10*log(10), break; end end