Last active
January 17, 2023 04:28
-
-
Save youheiakimoto/ab51e88c73baf68effd95b750100aad0 to your computer and use it in GitHub Desktop.
Adversarial-CMA-ES: derivative-free min-max optimization solver using (1+1)-CMA-ES.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| from scipy import stats | |
| import matplotlib.pyplot as plt | |
| """Adversarial-CMA-ES: derivative-free min-max optimization solver using (1+1)-CMA-ES [1] | |
| It is an extention of Adversarial Evolution Strategy [2] | |
| (https://gist.github.com/youheiakimoto/15212dbf46dc546af20af38b0b48ff17) | |
| Reference | |
| --------- | |
| [1] Youhei Akimoto, Yoshiki Miyauchi, and Atsuo Maki. 2022. | |
| Saddle Point Optimization with Approximate Minimization Oracle and Its Application to Robust Berthing Control. | |
| ACM Trans. Evol. Learn. Optim. 2, 1, Article 2 (March 2022), 32 pages. | |
| DOI:https://doi.org/10.1145/3510425 | |
| [2] Youhei Akimoto. 2021. | |
| Saddle point optimization with approximate minimization oracle. | |
| In Proceedings of the Genetic and Evolutionary Computation Conference (GECCO '21). | |
| Association for Computing Machinery, New York, NY, USA, 493–501. | |
| DOI:https://doi.org/10.1145/3449639.3459266 | |
| """ | |
| class Cmaes: | |
| """(1+1)-Active-CMA-ES [Arnold 2010] with Simplified 1/5th Success Rule""" | |
| def __init__(self, dim, tau, tau_offset, tolsigma=0.0, tolh=-np.inf, verbose=False): | |
| self.verbose = verbose | |
| self.tolsigma = tolsigma | |
| self.tolh = tolh | |
| self.dim = dim | |
| self.tau = tau | |
| self.tau_offset = tau_offset | |
| self.rescale_itr = 0 | |
| self.aup = np.exp(2.0 / (2.0 + self.dim)) | |
| self.adown = self.aup ** (-0.25) | |
| self.pthresh = 0.44 | |
| self.cp = 1.0 / 12.0 | |
| self.cc = 2.0 / (self.dim + 2) | |
| self.ccovp = 2.0 / (self.dim**2 + 6) | |
| self.ccovm = 0.4 / (self.dim**1.6 + 1) | |
| self.x = np.zeros(self.dim) | |
| self.s = 1.0 | |
| self.psucc = 0.5 | |
| self.fhist = [] | |
| self.p = np.zeros(self.dim) | |
| self.chol = np.eye(self.dim) | |
| self.cholinv = np.eye(self.dim) | |
| def minimize(self, h, hx, x, sigma, chol, maxeval): | |
| """Locate approximate solutions to min_{x} h(x) | |
| Parameters | |
| ---------- | |
| h : callable | |
| objective | |
| hx : float | |
| h(x) value | |
| x : ndarray (1d) | |
| initial solution | |
| sigma : float | |
| step-size | |
| chol : ndarray (2d) | |
| cholesky factor | |
| maxeval : int | |
| maximal number of f-calls | |
| Returns | |
| ------- | |
| hz (float) : h(zt) | |
| zt (ndarray, 1d) : approximate solution | |
| theta (tuple) : parameters for z | |
| neval (int) : number of h-calls | |
| """ | |
| nsucc = 0 | |
| neval = 0 | |
| self.x = np.copy(x) | |
| self.s = sigma | |
| self.psucc = 0.5 | |
| self.fhist = [hx] | |
| self.p = np.zeros(self.dim) | |
| self.chol = np.copy(chol) | |
| self.cholinv = np.linalg.inv(self.chol) | |
| self.rescale() | |
| self.rescale_itr = 0 | |
| while nsucc < self.dim * self.tau + self.tau_offset: | |
| z_, y_, x_ = self.generate() | |
| hx = h(x_) | |
| self.update(z_, y_, x_, hx) | |
| succ = hx <= self.fhist[0] | |
| nsucc += int(succ) | |
| neval += 1 | |
| if self.rescale_itr > self.dim: | |
| self.rescale() | |
| self.rescale_itr = 0 | |
| else: | |
| self.rescale_itr += 1 | |
| if self.verbose: | |
| print(neval, hx, self.s) | |
| if self.s <= self.tolsigma or hx <= self.tolh or neval >= maxeval: | |
| break | |
| return self.fhist[0], self.x, max(self.s, self.tolsigma), self.chol, neval | |
| def generate(self): | |
| z = np.random.randn(self.dim) | |
| y = np.dot(self.chol, z) | |
| x = self.x + self.s * y | |
| return z, y, x | |
| def rescale(self): | |
| frob = np.sum(self.chol ** 2) | |
| scale = np.sqrt(frob / self.dim) | |
| self.s *= scale | |
| self.chol /= scale | |
| self.cholinv *= scale | |
| self.p /= scale | |
| def update(self, z, y, x, fx): | |
| if fx <= self.fhist[0]: | |
| self.fhist = [fx] + self.fhist[:min(4, len(self.fhist))] | |
| self.x = x | |
| self.s *= self.aup | |
| self.psucc = (1.0 - self.cp) * self.psucc + self.cp | |
| if self.psucc > self.pthresh: | |
| self.p = (1.0 - self.cc) * self.p | |
| d = self.ccovp * (1.0 - self.cc * (2.0 - self.cc)) | |
| else: | |
| self.p = (1.0 - self.cc) * self.p + np.sqrt(self.cc * (2.0 - self.cc)) * y | |
| d = self.ccovp | |
| w = np.dot(self.cholinv, self.p) | |
| wnorm2 = np.dot(w, w) | |
| if wnorm2 > 1e-6: | |
| a = np.sqrt(1.0 - d) | |
| b = np.sqrt(1.0 - d) * (np.sqrt(1.0 + self.ccovp * wnorm2 / (1.0 - d)) - 1.0) / wnorm2 | |
| self.chol = a * self.chol + b * np.outer(self.p, w) | |
| self.cholinv = self.cholinv / a - (b / (a**2 + a * b * wnorm2)) * np.outer(w, np.dot(w, self.cholinv)) | |
| else: | |
| self.s *= self.adown | |
| self.psucc = (1.0 - self.cp) * self.psucc | |
| if len(self.fhist) > 4 and fx > self.fhist[4] and self.psucc <= self.pthresh: | |
| znorm2 = np.dot(z, z) | |
| ccov = self.ccovm if 1.0 >= self.ccovm * (2.0 * znorm2 - 1.0) else 1.0 / (2.0 * znorm2 - 1.0) | |
| a = np.sqrt(1.0 + ccov) | |
| b = np.sqrt(1.0 + ccov) * (np.sqrt(1.0 - ccov * znorm2 / (1.0 + ccov)) - 1) / znorm2 | |
| self.chol = a * self.chol + b * np.outer(y, z) | |
| self.cholinv = self.cholinv / a - (b / (a**2 + a * b * znorm2)) * np.outer(z, np.dot(z, self.cholinv)) | |
| class AdversarialCmaes: | |
| def __init__(self, f, x_dim, y_dim, x_tau=5, y_tau=5, x_tau_offset=5, y_tau_offset=5, tolgap=0, aeta=1, beta=5, ceta=1.1, xsampler=None, ysampler=None, xsigma_min=0, ysigma_min=0, xlbound=None, xubound=None, ylbound=None, yubound=None, logpath=None, logchol=True): | |
| """Adversarial-CMA-ES | |
| Parameters | |
| ---------- | |
| f : callable | |
| objective function | |
| x_dim, y_dim : int | |
| dimension of x and y | |
| x_tau, y_tau, x_tau_offset, y_tau_offset : int, optional | |
| number of iterations of the internal (1+1)-CMA-ES is | |
| x_tau * x_dim + x_tau_offset (for x) and | |
| y_tau * y_dim + y_tau_offset (for y). | |
| tolgap : float >= 0, optional | |
| termination condition. | |
| stop if the estimated suboptimality error is smaller than `tolgap` | |
| aeta, beta, ceta : float, optional | |
| hyper-parameters of the learning rate adaptation | |
| aeta, beta : the number of oracle calls to estimate the convergence rate is aeta / eta + beta | |
| ceta : learning rate is increased or decreased by the factor ceta | |
| xsampler, ysampler : callable, optional | |
| x and y samplers | |
| tolxsigma, tolysigma : float >= 0, optional | |
| minimal step-size | |
| xlbound, xubound : numpy.ndarray (1d) | |
| lower and upper bound of the search interval for x | |
| ylbound, yubound : numpy.ndarray (1d) | |
| lower and upper bound of the search interval for y | |
| logpath : str, optional | |
| log file name | |
| logchol : bool, optional | |
| write out the cholesky factor of the CMA-ES into the log file if this is true | |
| """ | |
| self._f = f | |
| self.xsampler = xsampler | |
| self.ysampler = ysampler | |
| self.ysigma_min = ysigma_min | |
| self.x_cma = Cmaes(x_dim, x_tau, x_tau_offset, tolsigma=xsigma_min) | |
| self.y_cma = Cmaes(y_dim, y_tau, y_tau_offset, tolsigma=ysigma_min) | |
| self.tolgap = tolgap | |
| self.aeta = aeta | |
| self.beta = beta | |
| self.ceta = ceta | |
| self.xlbound = xlbound | |
| self.xubound = xubound | |
| self.ylbound = ylbound | |
| self.yubound = yubound | |
| self.logpath = logpath | |
| self.logchol = logchol | |
| self.neval = 0 | |
| self.yhist = [] | |
| def f(self, x, y): | |
| fxy = self._f(x, y) | |
| self.neval += 1 | |
| return fxy | |
| def fyhist(self, x, y, twooutput=False): | |
| fxy = self.f(x, y) | |
| if len(self.yhist) > 0: | |
| fxyhist = np.max([self.f(x, yy) for yy in self.yhist]) | |
| if twooutput: | |
| return max(fxyhist, fxy), fxy | |
| else: | |
| return max(fxyhist, fxy) | |
| else: | |
| if twooutput: | |
| return fxy, fxy | |
| else: | |
| return fxy | |
| def optimize_with(self, eta, x0, y0, xsigma0, ysigma0, xchol0, ychol0, nstep, maxeval, minitr): | |
| """Adversarial-CMA-ES with Fixed Learning Rate | |
| Parameters | |
| ---------- | |
| eta : float | |
| learning rate | |
| x0, y0 : ndarray (1d) | |
| initial x and x | |
| xsigma0, ysigma0 : float | |
| initial step-size | |
| xchol0, ychol0 : ndarray (2d) | |
| initial cholesky factor | |
| nstep : int | |
| maximum iterations | |
| maxeval : int | |
| maximum number of f-calls | |
| minitr : int | |
| minimum iterations | |
| Returns | |
| ------- | |
| x (ndarray, 1D) : solution | |
| y (ndarray, 1D) : solution | |
| gap (ndarray, 1D) : history of gap: f(x', yt) - f(xt, y') | |
| """ | |
| x = np.copy(x0) | |
| y = np.copy(y0) | |
| xp = np.copy(x) | |
| yp = np.copy(y) | |
| x_sigma = xsigma0 | |
| y_sigma = ysigma0 | |
| x_chol = np.copy(xchol0) | |
| y_chol = np.copy(ychol0) | |
| dim_theta_x = 1 + len(x) ** 2 if self.logchol else 1 | |
| dim_theta_y = 1 + len(y) ** 2 if self.logchol else 1 | |
| res = np.zeros(1+1+1+1+len(x)+len(y)+dim_theta_x+dim_theta_y) | |
| gap_array = np.zeros(nstep) | |
| for i in range(nstep): | |
| # Determine whether or not to start from x or xp (y or yp) | |
| fxyhist, fxy = self.fyhist(x, y, twooutput=True) | |
| fxpy = self.fyhist(xp, y) | |
| fxyp = self.f(x, yp) | |
| if fxpy > fxyhist: | |
| xp[:] = x[:] | |
| fxpy = fxyhist | |
| x_sigma = xsigma0 | |
| x_chol = np.copy(xchol0) | |
| if fxyp < fxy: | |
| yp[:] = y[:] | |
| fxyp = fxy | |
| y_sigma = ysigma0 | |
| y_chol = np.copy(ychol0) | |
| # Adversarial steps | |
| fx, xp, x_sigma, x_chol, nx = self.x_cma.minimize(lambda z: self.fyhist(z, y), fxpy, xp, x_sigma, x_chol, maxeval - self.neval) | |
| fy, yp, y_sigma, y_chol, ny = self.y_cma.minimize(lambda z: -self.f(x, z), -fxyp, yp, y_sigma, y_chol, maxeval - self.neval) | |
| fy = -fy | |
| # Mirror solutions if the domain is bounded | |
| if self.xlbound is not None: | |
| xp = mirror(xp, self.xlbound, self.xubound) | |
| if self.ylbound is not None: | |
| yp = mirror(yp, self.ylbound, self.yubound) | |
| # Try random samples | |
| if self.xsampler is not None: | |
| xtrial = self.xsampler() | |
| fxytrial = self.fyhist(xtrial, y) | |
| nx += 1 | |
| if fxytrial < fx: | |
| xp = xtrial | |
| fx = fxytrial | |
| if self.ysampler is not None: | |
| ytrial = self.ysampler() | |
| fxytrial = self.f(x, ytrial) | |
| ny += 1 | |
| if fxytrial > fy: | |
| if fy >= fxyhist: | |
| self.register(yp) | |
| yp = ytrial | |
| fy = fxytrial | |
| # Estimated Suboptimality Error | |
| gap = max(fxyhist, fy) - fx | |
| # Update | |
| x += eta * (xp - x) | |
| y += eta * (yp - y) | |
| # Output | |
| idx = 0 | |
| res[idx] = self.neval; idx += 1 | |
| res[idx] = fxyhist; idx += 1 | |
| res[idx] = gap; idx += 1 | |
| res[idx] = eta; idx += 1 | |
| res[idx:idx+len(x)] = x; idx += len(x) | |
| res[idx:idx+len(y)] = y; idx += len(y) | |
| res[idx] = x_sigma; idx += 1 | |
| if self.logchol: | |
| res[idx:idx+dim_theta_x-1] = x_chol.flatten(); idx += dim_theta_x - 1 | |
| res[idx] = y_sigma; idx += 1 | |
| if self.logchol: | |
| res[idx:idx+dim_theta_x-1] = y_chol.flatten(); idx += dim_theta_y - 1 | |
| with open(self.logpath, 'ba') as flog: | |
| np.savetxt(flog, res, newline=' ') | |
| flog.write(b"\n") | |
| # Termination condition | |
| gap_array[i] = gap | |
| if self.neval >= maxeval: | |
| break | |
| if gap_array[i] < self.tolgap: | |
| break | |
| if i >= minitr - 1: | |
| gap = gap_array[i+1-minitr:i+1] | |
| if minitr > 1 and np.all(gap[1:minitr] - gap[0:minitr-1] > 0): | |
| break | |
| return x, y, gap_array[:i+1] | |
| def optimize(self, etamin, x0, y0, xsigma0, ysigma0, xchol0, ychol0, maxeval): | |
| """Adversarial Evolution Strategy with Learning Rate Adaptation | |
| Parameters | |
| ---------- | |
| etamin : float | |
| minimal learning rate | |
| x0, y0 : ndarray (1d) | |
| initial x and x | |
| xsigma0, ysigma0 : float | |
| initial step-size | |
| xchol0, ychol0 : ndarray (2d) | |
| initial cholesky factor | |
| maxeval : int | |
| maximum number of f-calls | |
| Returns | |
| ------- | |
| x (ndarray, 1D) : solution | |
| y (ndarray, 1D) : solution | |
| """ | |
| x = np.copy(x0) | |
| y = np.copy(y0) | |
| x_sigma = xsigma0 | |
| y_sigma = ysigma0 | |
| x_chol = np.copy(xchol0) | |
| y_chol = np.copy(ychol0) | |
| eta = 1.0 | |
| slp = 0.0 | |
| while self.neval < maxeval: | |
| u = np.random.rand() | |
| if u < 1/3: | |
| eta1 = eta * self.ceta | |
| elif 1/3 <= u < 2/3: | |
| eta1 = eta | |
| else: | |
| eta1 = eta / self.ceta | |
| eta1 = np.clip(eta1, etamin, 1) | |
| nstep = self.beta + int(self.aeta/eta1) | |
| x1, y1, gap = self.optimize_with(eta1, x, y, x_sigma, y_sigma, x_chol, y_chol, nstep, maxeval, self.beta) | |
| if gap[-1] < self.tolgap or self.neval >= maxeval: | |
| break | |
| xarr = np.arange(len(gap)) | |
| slp1, _, _, _, std1 = stats.linregress(xarr, np.log(np.fmax(gap, self.tolgap))) | |
| if slp1 - 2*std1 < 0: | |
| x = x1 | |
| y = y1 | |
| x_sigma = self.x_cma.s | |
| y_sigma = self.y_cma.s | |
| x_chol = self.x_cma.chol | |
| y_chol = self.y_cma.chol | |
| if slp >= 0 and slp1 >= 0: | |
| eta = np.clip(eta / self.ceta**3, etamin, 1) | |
| slp = 0.0 | |
| elif slp1 <= slp: | |
| eta = eta1 | |
| slp = slp1 | |
| elif eta == eta1: | |
| slp = slp1 | |
| return x1, y1 | |
| def register(self, yp): | |
| if len(self.yhist) > 0: | |
| dist = np.asarray([np.linalg.norm(yh - yp) for yh in self.yhist]) | |
| if np.min(dist) > self.ysigma_min * np.sqrt(len(yp)): | |
| self.yhist += [np.copy(yp)] | |
| else: | |
| self.yhist += [np.copy(yp)] | |
| def mirror(z, lbound, ubound): | |
| width = ubound - lbound | |
| return ubound - np.abs(np.mod(z - lbound, 2 * width) - width) | |
| def generate_cmaes_parameter(lbound, ubound): | |
| dim = len(lbound) | |
| x = lbound + (ubound - lbound) * np.random.rand(dim) | |
| sigma = 1.0 | |
| chol = np.diag( (ubound - lbound) / 4. ) | |
| return x, sigma, chol | |
| def minmax(f, xlbound, xubound, ylbound, yubound, maxeval, bounded, etamin=1e-4, tolgap=1e-6, tolxsigma=1e-8, tolysigma=1e-8, logpath="advcma.csv", logchol=True, initx=None, initxsigma=None, inity=None, initysigma=None): | |
| """Optimize min-max problem f with Adversarial-CMA-ES with restart | |
| Parameters | |
| ---------- | |
| f : callable | |
| min-max objective function f(x, y) | |
| xlbound, xubound : numpy.ndarray (1d) | |
| lower and upper bound of the initialization interval for x | |
| ylbound, yubound : numpy.ndarray (1d) | |
| lower and upper bound of the initialization interval for y | |
| maxeval : int | |
| maximum number of f-calls | |
| bounded : bool | |
| the domain is bounded (== initialization interval) if it is True | |
| etamin : float >= 0, optional | |
| minimal accepted learning rate | |
| tolgap : float >= 0, optional | |
| restart is performed if the suboptimality error reach this tolerance value | |
| tolxsigma, tolysigma : float >= 0, optional | |
| minimal step-size | |
| logpath : str, optional | |
| log file name | |
| logchol : bool, optional | |
| write out the cholesky factor of the CMA-ES into the log file if this is true | |
| initx, inity : numpy.ndarray (1d), optional | |
| initial guess | |
| initxsigma, initysigma : float > 0, optional | |
| initial step size, valid only when initx and inity are given | |
| Returns | |
| ------- | |
| float : f(x, y) value | |
| numpy.ndarray (1d) : x best | |
| numpy.ndarray (1d) : y worst | |
| numpy.ndarray (2d) : f(x_i, y_j) | |
| list of numpy.ndarray (1d) : list of candidate x | |
| list of numpy.ndarray (1d) : list of candidate y | |
| """ | |
| def xsampler(): | |
| return xlbound + (xubound - xlbound) * np.random.rand(len(xlbound)) | |
| def ysampler(): | |
| return ylbound + (yubound - ylbound) * np.random.rand(len(ylbound)) | |
| xhist = [] | |
| if bounded: | |
| xl, xu, yl, yu = xlbound, xubound, ylbound, yubound | |
| def ff(x, y): | |
| xmirror = mirror(x, xlbound, xubound) | |
| ymirror = mirror(y, ylbound, yubound) | |
| return f(xmirror, ymirror) | |
| else: | |
| xl, xu, yl, yu = None, None, None, None | |
| ff = f | |
| advcma = AdversarialCmaes(ff, len(xlbound), len(ylbound), tolgap=tolgap, xsampler=xsampler, ysampler=ysampler, xlbound=xl, xubound=xu, ylbound=yl, yubound=yu, xsigma_min=tolxsigma, ysigma_min=tolysigma, logpath=logpath, logchol=logchol) | |
| while advcma.neval < maxeval: | |
| x0, xsigma0, xchol0 = generate_cmaes_parameter(xlbound, xubound) | |
| if len(xhist) == 0: | |
| if initx is not None: | |
| x0 = np.array(initx, copy=True) | |
| if initxsigma is not None: | |
| xsigma0 = initxsigma | |
| y0, ysigma0, ychol0 = generate_cmaes_parameter(ylbound, yubound) | |
| if len(xhist) == 0: | |
| if inity is not None: | |
| y0 = np.array(inity, copy=True) | |
| if initysigma is not None: | |
| ysigma0 = initysigma | |
| x, y = advcma.optimize(etamin, x0, y0, xsigma0, ysigma0, xchol0, ychol0, maxeval) | |
| xhist += [x] | |
| advcma.register(y) | |
| yhist = advcma.yhist | |
| f_arr = np.zeros((len(xhist), len(yhist))) | |
| for i, x in enumerate(xhist): | |
| for j, y in enumerate(yhist): | |
| f_arr[i, j] = f(x, y) | |
| xidx = np.argmin(np.max(f_arr, axis=1)) | |
| yidx = np.argmax(f_arr[xidx]) | |
| return f_arr[xidx, yidx], xhist[xidx], yhist[yidx], f_arr, xhist, yhist | |
| def minimize(f, lbound, ubound, bounded, maxeval, totalmaxeval, tolsigma=1e-8, tolf=-np.inf): | |
| """Locate approximate solutions to min_{x} f(x) | |
| Parameters | |
| ---------- | |
| f : callable | |
| objective | |
| lbound, ubound : numpy.ndarray (1d) | |
| lower and upper bound of the initialization interval for x | |
| bounded : bool | |
| the domain is bounded (== initialization interval) if it is True | |
| maxeval : int | |
| maximal number of f-calls for each (1+1)-CMA-ES run | |
| totalmaxeval : int | |
| maximal number of f-calls for restart | |
| tolxigma : float >= 0, optional | |
| minimal step-size | |
| tolf : float, optional | |
| minimal objective value | |
| Returns | |
| ------- | |
| float : best f(x) value | |
| numpy.ndarray (1d) : best x | |
| numpy.ndarray (1d) : list of f(x) | |
| numpy.ndarray (2d) : list of x | |
| """ | |
| if bounded: | |
| def h(z): | |
| zmirror = mirror(z, lbound, ubound) | |
| return f(zmirror) | |
| else: | |
| def h(z): | |
| return f(z) | |
| f_arr = [] | |
| z_arr = [] | |
| nfev = 0 | |
| while nfev < totalmaxeval: | |
| z, sigma, chol = generate_cmaes_parameter(lbound, ubound) | |
| cma = Cmaes(len(z), tau=0, tau_offset=maxeval, tolsigma=tolsigma, tolh=tolf) | |
| hz, z, _, _, neval = cma.minimize(h, h(z), z, sigma, chol, min(maxeval, totalmaxeval - nfev - 1)) | |
| nfev += neval + 1 | |
| f_arr += [hz] | |
| z_arr += [mirror(z, lbound, ubound) if bounded else z] | |
| zidx = np.argmin(f_arr) | |
| return f_arr[zidx], z_arr[zidx], f_arr, z_arr | |
| def worstsearch(f, x, ylbound, yubound, bounded, maxeval, tolsigma=1e-8, n_restart=100, inity=None, initysigma=None): | |
| """Worst case search by the (1+1)-CMA-ES | |
| Parameters | |
| ---------- | |
| f : callable | |
| min-max objective function | |
| x : numpy.ndarray (1d) | |
| x to be evaluated | |
| ylbound, yubound : numpy.ndarray (1d) | |
| lower and upper bound of the initialization interval for y | |
| maxeval : int | |
| maximum number of f-calls for each restart | |
| tolsigma : float >= 0, optional | |
| minimal step-size | |
| n_restart : int >= 1, optional | |
| number of restart | |
| inity : numpy.ndarray (1d), optional | |
| initial solution | |
| initysigma : float > 0, optional | |
| initial step-size, valid only when inity is given | |
| Returns | |
| ------- | |
| float : worst f(x, y) value | |
| numpy.ndarray (1d) : worst y | |
| numpy.ndarray (1d) : f(x, y_j) | |
| numpy.ndarray (2d) : list of candidate y | |
| """ | |
| if bounded: | |
| def h(y): | |
| ymirror = mirror(y, ylbound, yubound) | |
| return -f(x, ymirror) | |
| else: | |
| def h(y): | |
| return -f(x, y) | |
| f_arr = np.zeros(n_restart) | |
| y_arr = np.zeros((n_restart, len(ylbound))) | |
| for i in range(n_restart): | |
| y, sigma, chol = generate_cmaes_parameter(ylbound, yubound) | |
| if i == 0: | |
| if inity is not None: | |
| y = np.array(inity) | |
| if initysigma is not None: | |
| sigma = initysigma | |
| cma = Cmaes(len(y), tau=0, tau_offset=maxeval, tolsigma=tolsigma) | |
| hy, y, _, _, _ = cma.minimize(h, h(y), y, sigma, chol, maxeval) | |
| f_arr[i] = -hy | |
| y_arr[i] = mirror(y, ylbound, yubound) if bounded else y | |
| yidx = np.argmax(f_arr) | |
| return f_arr[yidx], y_arr[yidx], f_arr, y_arr | |
| def ex1(m=10, a=1, b=4, c=1): | |
| """experiments on f1""" | |
| n = m | |
| def f1(x, y): | |
| fx = (a/2) * np.dot(x, x) | |
| fxy = b * np.dot(x, y) | |
| fy = (c/2) * np.dot(y, y) | |
| return fx + fxy - fy | |
| def f1gap(x, y): | |
| return (1/2 + b**2 / (2 * a * c)) * (np.dot(x, x) + np.dot(y, y)) | |
| f = f1 | |
| fgap = f1gap | |
| # Initialization Interval | |
| xlbound = -np.ones(m) | |
| xubound = 5 *np.ones(m) | |
| ylbound = -np.ones(n) | |
| yubound = 5 * np.ones(n) | |
| bounded = False | |
| # Optimization | |
| maxeval = 1e6 | |
| logpre = 'advcma_f1' | |
| logpath = logpre + '.csv' | |
| logpathx = logpre + '_x.csv' | |
| logpathy = logpre + '_y.csv' | |
| logpathworst = logpre + '_worst_y.csv' | |
| print("Optimization has started.") | |
| fxy, xbest, yworst, f_arr, x_arr, y_arr = minmax(f, xlbound, xubound, ylbound, yubound, maxeval, bounded, etamin=1e-4, tolgap=1e-6, tolxsigma=1e-8, tolysigma=1e-8, logpath=logpath, logchol=False) | |
| print("Optimization has finished.") | |
| print("best x:", xbest) | |
| print("worst y:", yworst) | |
| print("f(x,y):", fxy) | |
| print(f_arr) | |
| np.savetxt(logpathx, x_arr) | |
| np.savetxt(logpathy, y_arr) | |
| # Evaluation | |
| print("Gap evaluation has started.") | |
| dat = np.loadtxt(logpath) | |
| gap_array = np.zeros(len(dat)) | |
| for i in range(len(gap_array)): | |
| x = dat[i, 4:4+m] | |
| y = dat[i, 4+m:4+m+n] | |
| gap_array[i] = fgap(x, y) | |
| print("Gap evaluation has finished.") | |
| # Plot | |
| fig = plt.figure() | |
| plt.semilogy(dat[:, 0], gap_array, label="max_y f(x, y) - min_x f(x, y)") | |
| plt.semilogy(dat[:, 0], dat[:, 2], label="max_y f(x, y) - min_x f(x, y) (estimated)") | |
| plt.semilogy(dat[:, 0], dat[:, 3], label="eta") | |
| plt.semilogy(dat[:, 0], dat[:, -2], label="sigma_x") | |
| plt.semilogy(dat[:, 0], dat[:, -1], label="sigma_y") | |
| plt.grid() | |
| plt.xlabel("#f-calls") | |
| plt.legend() | |
| plt.tight_layout() | |
| plt.savefig(logpre + ".pdf") | |
| def ex2(funcid=2, m=50, n=20): | |
| """experiments on f2--f7""" | |
| # Initialization Interval | |
| xlbound = -np.ones(m) | |
| xubound = 5 *np.ones(m) | |
| ylbound = -np.ones(n) | |
| yubound = 5 * np.ones(n) | |
| bounded = True | |
| def f2(x, y): | |
| xx = np.dot(x, x) / 2.0 | |
| xy = np.sum(x) * np.sum(y) / m | |
| yy = np.dot(y, y) / 2.0 | |
| return xx + xy - yy | |
| def f2max(x): | |
| y = np.ones(n) * (np.sum(x) / m) | |
| return f2(x, y) | |
| def f3(x, y): | |
| xx = min(np.dot(x, x), np.dot(x-4, x-4)) / 2.0 | |
| xy = np.sum(x) * np.sum(y) / m | |
| yy = np.dot(y, y) / 2 | |
| return xx + xy - yy | |
| def f3max(x): | |
| y = np.ones(n) * (np.sum(x) / m) | |
| return f3(x, y) | |
| def f4(x, y): | |
| xx = np.dot(x, x) / 2.0 | |
| xy = np.sum(x) * np.sum(y) / m | |
| yy = (np.dot(y, y) / 2.0)**2 | |
| return xx + xy - yy | |
| def f4max(x): | |
| y = np.ones(n) * (np.abs(np.sum(x)) / (m * n)) ** (1.0/3.0) * np.sign(np.sum(x)) | |
| return f4(x, y) | |
| def f5(x, y): | |
| xy = np.linalg.norm(x) * np.sum(y) / np.sqrt(m) | |
| return xy | |
| def f5max(x): | |
| y = yubound | |
| return f5(x, y) | |
| def f6(x, y): | |
| xy = np.sum(x) * np.sum(y) / m | |
| yy = np.dot(y, y) / 2.0 | |
| return xy - yy | |
| def f6max(x): | |
| y = np.ones(n) * (np.sum(x) / m) | |
| return f6(x, y) | |
| def f7(x, y): | |
| xx = np.dot(x, x) / 2.0 | |
| xy = np.sum(x) * np.sum(y) / m | |
| return xx + xy | |
| def f7max(x): | |
| if np.sum(x) > 0: | |
| y = yubound | |
| else: | |
| y = ylbound | |
| return f7(x, y) | |
| assert funcid >= 2 and funcid <= 7 | |
| f = [None, None, f2, f3, f4, f5, f6, f7][funcid] | |
| fmax = [None, None, f2max, f3max, f4max, f5max, f6max, f7max][funcid] | |
| # Optimization | |
| maxeval = 1e6 | |
| logpre = 'advcma_f{}'.format(str(funcid)) | |
| logpath = logpre + '.csv' | |
| logpathx = logpre + '_x.csv' | |
| logpathy = logpre + '_y.csv' | |
| logpathworst = logpre + '_worst_y.csv' | |
| print("Optimization has started.") | |
| fxy, xbest, yworst, f_arr, x_arr, y_arr = minmax(f, xlbound, xubound, ylbound, yubound, maxeval, bounded, etamin=1e-4, tolgap=1e-6, tolxsigma=1e-8, tolysigma=1e-8, logpath=logpath, logchol=False) | |
| print("Optimization has finished.") | |
| print("best x:", xbest) | |
| print("worst y:", yworst) | |
| print("f(x,y):", fxy) | |
| print(f_arr) | |
| np.savetxt(logpathx, x_arr) | |
| np.savetxt(logpathy, y_arr) | |
| # Evaluation | |
| print("Worst objective evaluation has started.") | |
| dat = np.loadtxt(logpath) | |
| fmax_array = np.zeros(len(dat)) | |
| for i in range(len(fmax_array)): | |
| x = dat[i, 4:4+m] | |
| fmax_array[i] = fmax(x) | |
| print("Worst objective evaluation has finished.") | |
| # Plot | |
| fig = plt.figure() | |
| plt.semilogy(dat[:, 0], fmax_array, label="max_y f(x, y)") | |
| plt.semilogy(dat[:, 0], dat[:, 2], label="max_y f(x, y) - min_x f(x, y) (estimated)") | |
| plt.semilogy(dat[:, 0], dat[:, 3], label="eta") | |
| plt.semilogy(dat[:, 0], dat[:, -2], label="sigma_x") | |
| plt.semilogy(dat[:, 0], dat[:, -1], label="sigma_y") | |
| plt.grid() | |
| plt.xlabel("#f-calls") | |
| plt.legend() | |
| plt.tight_layout() | |
| plt.savefig(logpre + ".pdf") | |
| def my_ex(): | |
| """experiment prototype""" | |
| def ftest(x, y, a=1.0, b=10.0, c=1.0): | |
| """define your objective function""" | |
| fx = (a/2) * np.dot(x, x) | |
| fxy = b * np.dot(x, y) | |
| fy = (c/2) * np.dot(y, y) | |
| return fx + fxy - fy | |
| np.random.seed(100) | |
| m = n = 10 | |
| a = 1 | |
| b = 4 | |
| c = 1 | |
| bounded = False | |
| # Initialization Interval | |
| xlbound = -np.ones(m) | |
| xubound = np.ones(m) | |
| ylbound = -np.ones(n) | |
| yubound = np.ones(n) | |
| def f(x, y): | |
| return ftest(x, y, a, b, c) | |
| # Optimization | |
| maxeval = 1e6 | |
| logpre = 'advcma_ftest' | |
| logpath = logpre + '.csv' | |
| logpathx = logpre + '_x.csv' | |
| logpathy = logpre + '_y.csv' | |
| logpathworst = logpre + '_worst_y.csv' | |
| print("Optimization has started.") | |
| fxy, xbest, yworst, f_arr, x_arr, y_arr = minmax(f, xlbound, xubound, ylbound, yubound, maxeval, bounded, etamin=1e-4, tolgap=1e-6, tolxsigma=1e-8, tolysigma=1e-8, logpath=logpath, logchol=False) | |
| print("Optimization has finished.") | |
| print("best x:", xbest) | |
| print("worst y:", yworst) | |
| print("f(x,y):", fxy) | |
| print(f_arr) | |
| np.savetxt(logpathx, x_arr) | |
| np.savetxt(logpathy, y_arr) | |
| # Evaluation | |
| print("Worst-case search has started.") | |
| fy, y, f_arr, y_arr = worstsearch(f, xbest, ylbound, yubound, bounded, maxeval=100*len(ylbound), tolsigma=1e-8, n_restart=100) | |
| print("Worst-case search has finished.") | |
| print("worst y:", y) | |
| print("f(x,y):", fy) | |
| print(f_arr) | |
| np.savetxt(logpathworst, y_arr) | |
| # Plot | |
| dat = np.loadtxt(logpath) | |
| fig = plt.figure() | |
| plt.semilogy(dat[:, 0], dat[:, 2], label="max_y f(x, y) - min_x f(x, y) (estimated)") | |
| plt.semilogy(dat[:, 0], dat[:, 3], label="eta") | |
| plt.semilogy(dat[:, 0], dat[:, -2], label="sigma_x") | |
| plt.semilogy(dat[:, 0], dat[:, -1], label="sigma_y") | |
| plt.grid() | |
| plt.xlabel("#f-calls") | |
| plt.legend() | |
| plt.tight_layout() | |
| plt.savefig("advcma_ftest.pdf") | |
| if __name__ == "__main__": | |
| import sys | |
| exid = int(sys.argv[1]) | |
| if exid == 0: | |
| my_ex() | |
| elif exid == 1: | |
| ex1() | |
| elif exid in [2, 3, 4, 5, 6, 7]: | |
| ex2(funcid=exid) | |
| else: | |
| raise ValueError("`exid` must be in [[0, 7]].") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
The
worstsearchfunction now accepts an initial search point and initial step-size to perform a local search for the first run.