diff --git a/brainiak/hyperparamopt/__init__.py b/brainiak/hyperparamopt/__init__.py new file mode 100644 index 000000000..7867be74a --- /dev/null +++ b/brainiak/hyperparamopt/__init__.py @@ -0,0 +1,4 @@ +""" Hyper parameter optimization package """ + +import pyximport +pyximport.install() diff --git a/brainiak/hyperparamopt/hpo.py b/brainiak/hyperparamopt/hpo.py new file mode 100644 index 000000000..19d13130b --- /dev/null +++ b/brainiak/hyperparamopt/hpo.py @@ -0,0 +1,369 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Hyper Parameter Optimization (HPO) + +This implementation is based on the work: + +.. [Bergstra2011] "Algorithms for Hyper-Parameter Optimization", + James S. Bergstra and Bardenet, R\'{e}mi and Bengio, Yoshua + and Bal\'{a}zs K\'{e}gl. NIPS 2011 + +.. [Bergstra2013] "Making a Science of Model Search: + Hyperparameter Optimization in Hundreds of Dimensions for + Vision Architectures", James Bergstra, Daniel Yamins, David Cox. + JMLR W&CP 28 (1) : 115–123, 2013 + +""" + +# Authors: Narayanan Sundaram (Intel Labs) + +import logging +import math +import numpy as np +from scipy.special import erf +import scipy.stats as st + + +logger = logging.getLogger(__name__) + + +def get_sigma(x, min_limit=-np.inf, max_limit=np.inf): + """Compute the standard deviations around the points for a 1D GMM. + + We take the distance from the nearest left and right neighbors + for each point, then use the max as the estimate of standard + deviation for the gaussian mixture around that point. + + Arguments + --------- + x : 1D array + Set of points to create the GMM + + min_limit : Optional[float], default : -inf + Minimum limit for the distribution + + max_limit : Optional[float], default : inf + maximum limit for the distribution + + Returns + ------- + 1D array + Array of standard deviations + """ + + z = np.append(x, [min_limit, max_limit]) + sigma = np.ones(x.shape) + for i in range(x.size): + # Calculate the nearest left neighbor of x[i] + # Find the minimum of (x[i] - k) for k < x[i] + xleft = z[np.argmin([(x[i] - k) if k < x[i] else np.inf for k in z])] + + # Calculate the nearest right neighbor of x[i] + # Find the minimum of (k - x[i]) for k > x[i] + xright = z[np.argmin([(k - x[i]) if k > x[i] else np.inf for k in z])] + + sigma[i] = max(x[i] - xleft, xright - x[i]) + if sigma[i] == np.inf: + sigma[i] = min(x[i] - xleft, xright - x[i]) + if (sigma[i] == -np.inf): # should never happen + sigma[i] = 1.0 + return sigma + + +class gmm_1d_distribution: + """GMM 1D distribution. + + Given a set of points, we create this object so that we + can calculate likelihoods and generate samples from this + 1D Gaussian mixture model. + + Attributes + ---------- + points : 1D array + Set of points to create the GMM + + N : int + Number of points to create the GMM + + min_limit : Optional[float], default : -inf + Minimum limit for the distribution + + max_limit : Optional[float], default : inf + Maximum limit for the distribution + + weights : Optional[1D array], default : array of ones + Used to weight the points non-uniformly if required + """ + + def __init__(self, x, min_limit=-np.inf, max_limit=np.inf, weights=1.0): + self.points = x + self.N = x.size + self.min_limit = min_limit + self.max_limit = max_limit + self.sigma = get_sigma(x, min_limit=min_limit, max_limit=max_limit) + self.weights = (2 + / (erf((max_limit - x) / (np.sqrt(2.) * self.sigma)) + - erf((min_limit - x) / (np.sqrt(2.) * self.sigma))) + * weights) + self.W_sum = np.sum(self.weights) + + def get_gmm_pdf(self, x): + """Calculate the GMM likelihood for a single point. + + .. math:: + y = \sum_{i=1}^{N} w_i*normpdf(x, x_i, \sigma_i)/\sum_{i=1}^{N} w_i + + Arguments + --------- + x : float + Point at which likelihood needs to be computed + + Returns + ------- + float + Likelihood value at x + """ + + def my_norm_pdf(xt, mu, sigma): + z = (xt - mu) / sigma + return (math.exp(-0.5 * z * z) + / (math.sqrt(2. * np.pi) * sigma)) + + y = 0 + if (x < self.min_limit): + return 0 + if (x > self.max_limit): + return 0 + for _x in range(self.points.size): + y += (my_norm_pdf(x, self.points[_x], self.sigma[_x]) + * self.weights[_x]) / self.W_sum + return y + + def __call__(self, x): + """Return the GMM likelihood for given point(s). + + .. math:: + y = \sum_{i=1}^{N} w_i*normpdf(x, x_i, \sigma_i)/\sum_{i=1}^{N} w_i + + Arguments + --------- + x : scalar (or) 1D array of reals + Point(s) at which likelihood needs to be computed + + Returns + ------- + scalar (or) 1D array + Likelihood values at the given point(s) + """ + + if np.isscalar(x): + return self.get_gmm_pdf(x) + else: + return np.array([self.get_gmm_pdf(t) for t in x]) + + def get_samples(self, n): + """Sample the GMM distribution. + + Arguments + --------- + n : int + Number of samples needed + + Returns + ------- + 1D array + Samples from the distribution + """ + + normalized_w = self.weights / np.sum(self.weights) + get_rand_index = st.rv_discrete(values=(range(self.N), + normalized_w)).rvs(size=n) + samples = np.zeros(n) + k = 0 + j = 0 + while (k < n): + i = get_rand_index[j] + j = j + 1 + if (j == n): + get_rand_index = st.rv_discrete(values=(range(self.N), + normalized_w)).rvs(size=n) + j = 0 + v = np.random.normal(loc=self.points[i], scale=self.sigma[i]) + if (v > self.max_limit or v < self.min_limit): + continue + else: + samples[k] = v + k = k + 1 + if (k == n): + break + return samples + + +def get_next_sample(x, y, min_limit=-np.inf, max_limit=np.inf): + """Get the next point to try, given the previous samples. + + We use [Bergstra2013]_ to compute the point that gives the largest + Expected improvement (EI) in the optimization function. This model fits 2 + different GMMs - one for points that have loss values in the bottom 15% + and another for the rest. Then we sample from the former distribution + and estimate EI as the ratio of the likelihoods of the 2 distributions. + We pick the point with the best EI among the samples that is also not + very close to a point we have sampled earlier. + + Arguments + --------- + x : 1D array + Samples generated from the distribution so far + + y : 1D array + Loss values at the corresponding samples + + min_limit : float, default : -inf + Minimum limit for the distribution + + max_limit : float, default : +inf + Maximum limit for the distribution + + Returns + ------- + float + Next value to use for HPO + """ + + z = np.array(list(zip(x, y)), dtype=np.dtype([('x', float), ('y', float)])) + z = np.sort(z, order='y') + n = y.shape[0] + g = int(np.round(np.ceil(0.15 * n))) + ldata = z[0:g] + gdata = z[g:n] + lymin = ldata['y'].min() + lymax = ldata['y'].max() + weights = (lymax - ldata['y']) / (lymax - lymin) + lx = gmm_1d_distribution(ldata['x'], min_limit=min_limit, + max_limit=max_limit, weights=weights) + gx = gmm_1d_distribution(gdata['x'], min_limit=min_limit, + max_limit=max_limit) + + samples = lx.get_samples(n=1000) + ei = lx(samples) / gx(samples) + + h = (x.max() - x.min()) / (10 * x.size) + # TODO + # assumes prior of x is uniform; should ideally change for other priors + # d = np.abs(x - samples[ei.argmax()]).min() + # CDF(x+d/2) - CDF(x-d/2) < 1/(10*x.size) then reject else accept + s = 0 + while (np.abs(x - samples[ei.argmax()]).min() < h): + ei[ei.argmax()] = 0 + s = s + 1 + if (s == samples.size): + break + xnext = samples[ei.argmax()] + + return xnext + + +def fmin(loss_fn, + space, + max_evals, + trials, + init_random_evals=30, + explore_prob=0.2): + """Find the minimum of function through hyper parameter optimization. + + Arguments + --------- + loss_fn : ``function(*args) -> float`` + Function that takes in a dictionary and returns a real value. + This is the function to be minimized. + + space : dictionary + Custom dictionary specifying the range and distribution of + the hyperparamters. + E.g. ``space = {'x': {'dist':scipy.stats.uniform(0,1), + 'lo':0, 'hi':1}}`` + for a 1-dimensional space with variable x in range [0,1] + + max_evals : int + Maximum number of evaluations of loss_fn allowed + + trials : list + Holds the output of the optimization trials. + Need not be empty to begin with, new trials are appended + at the end. + + init_random_evals : Optional[int], default 30 + Number of random trials to initialize the + optimization. + + explore_prob : Optional[float], default 0.2 + Controls the exploration-vs-exploitation ratio. Value should + be in [0,1]. By default, 20% of trails are random samples. + + Returns + ------- + trial entry (dictionary of hyperparameters) + Best hyperparameter setting found. + E.g. {'x': 5.6, 'loss' : 0.5} where x is the best hyparameter + value found and loss is the value of the function for the + best hyperparameter value(s). + + Raises + ------ + ValueError + If the distribution specified in space does not support a ``rvs()`` + method to generate random numbers, a ValueError is raised. + """ + + for s in space: + if not hasattr(space[s]['dist'], 'rvs'): + raise ValueError('Unknown distribution type for variable') + if 'lo' not in space[s]: + space[s]['lo'] = -np.inf + if 'hi' not in space[s]: + space[s]['hi'] = np.inf + + if len(trials) > init_random_evals: + init_random_evals = 0 + + for t in range(max_evals): + sdict = {} + + if t >= init_random_evals and np.random.random() > explore_prob: + use_random_sampling = False + else: + use_random_sampling = True + + yarray = np.array([tr['loss'] for tr in trials]) + for s in space: + sarray = np.array([tr[s] for tr in trials]) + if use_random_sampling: + sdict[s] = space[s]['dist'].rvs() + else: + sdict[s] = get_next_sample(sarray, yarray, + min_limit=space[s]['lo'], + max_limit=space[s]['hi']) + + logger.debug('Explore' if use_random_sampling else 'Exploit') + logger.info('Next point ', t, ' = ', sdict) + + y = loss_fn(sdict) + sdict['loss'] = y + trials.append(sdict) + + yarray = np.array([tr['loss'] for tr in trials]) + yargmin = yarray.argmin() + + logger.info('Best point so far = ', trials[yargmin]) + return trials[yargmin] diff --git a/examples/hyperparamopt/hpo_example.py b/examples/hyperparamopt/hpo_example.py new file mode 100644 index 000000000..e5fc030f8 --- /dev/null +++ b/examples/hyperparamopt/hpo_example.py @@ -0,0 +1,137 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Example for using hyperparameter optimization (hpo) package. + +In this example, we will try to optimize a function of +2 variables (branin) using both hpo and grid search. + +""" + +import brainiak.hyperparamopt.hpo as hpo +import numpy as np +import scipy.stats as st +import matplotlib.pyplot as plt + +# Branin is the function we want to minimize. +# It is a function of 2 variables. +# In the range x1 in [-5, 10] and x2 in [0, 15], +# this function has 2 local minima and 1 global minima. +# Global minima of -16.6 at (-3.7, 13.7). +# This is the modified version (Branin-Hoo) of the standard branin function. +# If you want the standard version (which has 3 global minima), +# you can omit the "+5*x1" term at the end +# For more details, see http://www.sfu.ca/~ssurjano/branin.html +def branin(x1, x2): + a = 1.0 + b = 5.1/(4*np.pi*np.pi) + c = 5.0/np.pi + r = 6.0 + s = 10.0 + t = 1.0/(8*np.pi) + return a*((x2 - b*x1*x1 + c*x1 - r)**2) + s*(1-t)*np.cos(x1) + s + 5*x1 + +# This is a wrapper around branin that takes in a dictionary +def branin_wrapper(args): + x1 = args['x1'] + x2 = args['x2'] + return branin(x1,x2) + +# Define ranges for the two variables +x1lo = -5 +x1hi = 10 +x2lo = 0 +x2hi = 15 + +############################## +# Optimization through hpo +############################## + +# Define a space for hpo to use +# The space needs to define +# 1. Name of the variables +# 2. Default samplers for the variables (use scipy.stats objects) +# 3. lo and hi ranges for the variables (will use -inf, inf if not specified) +space = {'x1':{'dist': st.uniform(x1lo, x1hi-x1lo), 'lo':x1lo, 'hi':x1hi}, + 'x2':{'dist': st.uniform(x2lo, x2hi-x2lo), 'lo':x2lo, 'hi':x2hi}} + +# The trials object is just a list that stores the samples generated and the +# corresponding function values at those sample points. +trials = [] + +# Maximum number of samples that will be generated. +# This is the maximum number of function evaluations that will be performed. +n_hpo_samples = 100 + +# Call the fmin function that does the optimization. +# The function to be optimized should take in a dictionary. You will probably +# need to wrap your function to do this (see branin() and branin_wrapper()). +# You can pass in a non-empty trials object as well e.g. from a previous +# fmin run. We just append to the trials object and will use existing data +# in our optimization. +print("Starting optimization through hpo") +best = hpo.fmin(loss_fn=branin_wrapper, space=space, + max_evals=n_hpo_samples, trials=trials) + +# Print out the best value obtained through HPO +print("Best obtained through HPO (", n_hpo_samples, " samples) = ", + best['x1'], best['x2'], "; min value = ", best['loss']) + +##################################### +# Optimization through grid search +##################################### + +# Divide the space into a uniform grid (meshgrid) +n = 200 +x1 = np.linspace(x1lo, x1hi, n) +x2 = np.linspace(x2lo, x2hi, n) +x1_grid, x2_grid = np.meshgrid(x1, x2) + +# Calculate the function values along the grid +print("Starting optimization through grid search") +z = branin(x1_grid, x2_grid) + +# Print out the best value obtained through grid search +print("Best obtained through grid search (", n*n, " samples) = ", + x1_grid.flatten()[z.argmin()], x2_grid.flatten()[z.argmin()], + "; min value = ", z.min()) + +######## +# Plots +######## + +# Convert trials object data into numpy arrays +x1 = np.array([tr['x1'] for tr in trials]) +x2 = np.array([tr['x2'] for tr in trials]) +y = np.array([tr['loss'] for tr in trials]) + +# Plot the function contour using the grid search data +h = (z.max()-z.min())/25 +plt.contour(x1_grid, x2_grid, z, levels=np.linspace(z.min()-h, z.max(), 26)) + +# Mark the points that were sampled through HPO +plt.scatter(x1, x2, s=10, color='r', label='HPO Samples') + +# Mark the best points obtained through both methods +plt.scatter(best['x1'], best['x2'], s=30, color='b', label='Best HPO') +plt.scatter(x1_grid.flatten()[z.argmin()], x2_grid.flatten()[z.argmin()], + s=30, color='g', label='Best grid search') + +# Labels +plt.xlabel('x1') +plt.ylabel('x2') +plt.title('Hyperparameter optimization using HPO (Branin function)') +plt.legend() +plt.show() + diff --git a/examples/hyperparamopt/requirements.txt b/examples/hyperparamopt/requirements.txt new file mode 100644 index 000000000..6ccafc3f9 --- /dev/null +++ b/examples/hyperparamopt/requirements.txt @@ -0,0 +1 @@ +matplotlib diff --git a/requirements-dev.txt b/requirements-dev.txt index 3503c58b9..bd99b5fb3 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,7 +1,9 @@ +cython flake8 flake8-print pytest pytest-cov +pytest-cython restructuredtext-lint sphinx sphinx_rtd_theme diff --git a/setup.cfg b/setup.cfg index f024f90d2..1069aeeb5 100644 --- a/setup.cfg +++ b/setup.cfg @@ -4,6 +4,7 @@ addopts = --cov-report=xml --cov-report=html --cov-report=term + --doctest-cython [coverage:run] branch = True diff --git a/tests/hyperparamopt/test_hpo.py b/tests/hyperparamopt/test_hpo.py new file mode 100644 index 000000000..64b9b7eab --- /dev/null +++ b/tests/hyperparamopt/test_hpo.py @@ -0,0 +1,87 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import pytest +import numpy as np +import scipy.stats as st +from brainiak.hyperparamopt.hpo import gmm_1d_distribution, fmin + + +def test_simple_gmm(): + x = np.array([1., 1., 2., 3., 1.]) + d = gmm_1d_distribution(x, min_limit=0., max_limit=4.) + assert d(1.1) > d(3.5), "GMM distribution not behaving correctly" + assert d(2.0) > d(3.0), "GMM distribution not behaving correctly" + assert d(-1.0) == 0, "GMM distribution out of bounds error" + assert d(9.0) == 0, "GMM distribution out of bounds error" + + samples = d.get_samples(n=25) + np.testing.assert_array_less(samples, 4.) + np.testing.assert_array_less(0., samples) + + +def test_simple_gmm_weights(): + x = np.array([1., 1., 2., 3., 1., 3.]) + d = gmm_1d_distribution(x) + + x2 = np.array([1., 2., 3.]) + w = np.array([3., 1., 2.]) + d2 = gmm_1d_distribution(x2, weights=w) + y2 = d2(np.array([1.1, 2.0])) + + assert d2(1.1) == y2[0],\ + "GMM distribution array & scalar results don't match" + assert np.abs(d(1.1) - d2(1.1)) < 1e-5,\ + "GMM distribution weights not handled correctly" + assert np.abs(d(2.0) - d2(2.0)) < 1e-5,\ + "GMM distribution weights not handled correctly" + + +def test_simple_hpo(): + + def f(args): + x = args['x'] + return x*x + + s = {'x': {'dist': st.uniform(loc=-10., scale=20), 'lo': -10., 'hi': 10.}} + trials = [] + + #Test fmin and ability to continue adding to trials + best = fmin(loss_fn=f, space=s, max_evals=40, trials=trials) + best = fmin(loss_fn=f, space=s, max_evals=10, trials=trials) + + assert len(trials) == 50, "HPO continuation trials not working" + + # Test verbose flag + best = fmin(loss_fn=f, space=s, max_evals=10, trials=trials) + + yarray = np.array([tr['loss'] for tr in trials]) + np.testing.assert_array_less(yarray, 100.) + + xarray = np.array([tr['x'] for tr in trials]) + np.testing.assert_array_less(np.abs(xarray), 10.) + + assert best['loss'] < 100., "HPO out of range" + assert np.abs(best['x']) < 10., "HPO out of range" + + #Test unknown distributions + s2 = {'x': {'dist': 'normal', 'mu': 0., 'sigma': 1.}} + trials2 = [] + with pytest.raises(ValueError) as excinfo: + best2 = fmin(loss_fn=f, space=s2, max_evals=40, trials=trials2) + assert "Unknown distribution type for variable" in str(excinfo.value) + + s3 = {'x': {'dist': st.norm(loc=0., scale=1.)}} + trials3 = [] + best3 = fmin(loss_fn=f, space=s3, max_evals=40, trials=trials3) +