7.3-Bayesian-Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "!pip install -r https://raw.githubusercontent.com/datamllab/automl-in-action-notebooks/master/requirements.txt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Load dataset\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_california_housing\n",
    "\n",
    "house_dataset = fetch_california_housing()\n",
    "\n",
    "# Import pandas package to format the data\n",
    "import pandas as pd\n",
    "\n",
    "# Extract features with their names into the a dataframe format\n",
    "data = pd.DataFrame(house_dataset.data, columns=house_dataset.feature_names)\n",
    "\n",
    "# Extract target with their names into a pd.Series object with name MEDV\n",
    "target = pd.Series(house_dataset.target, name=\"MEDV\")\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    data, target, test_size=0.2, random_state=42\n",
    ")\n",
    "X_train, X_val, y_train, y_val = train_test_split(\n",
    "    X_train, y_train, test_size=0.2, shuffle=False\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "X_train.shape, X_test.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Use LightGBM GBDT model to do regression without tuning\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "\n",
    "def build_model(hp):\n",
    "    model = lgb.LGBMRegressor(\n",
    "        boosting_type=\"gbdt\",\n",
    "        # you can also search model type such as:\n",
    "        # boosting_type=hp.Choice(\"model_type\", ['gbdt', 'goss'], default='gbdt'),\n",
    "        num_leaves=hp.Int(\"num_leaves\", 5, 50, step=1),\n",
    "        learning_rate=hp.Float(\"learning_rate\", 1e-3, 1, sampling=\"log\", default=0.01),\n",
    "        n_estimators=hp.Int(\"n_estimators\", 5, 50, step=1),\n",
    "    )\n",
    "\n",
    "    return model\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Customize tuner\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "\n",
    "import os\n",
    "import pickle\n",
    "import tensorflow as tf\n",
    "import keras_tuner as kt\n",
    "import lightgbm as lgb\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "\n",
    "class LightGBMTuner(kt.engine.base_tuner.BaseTuner):\n",
    "    def run_trial(self, trial, X, y, validation_data):\n",
    "        model = self.hypermodel.build(trial.hyperparameters)  # build the model\n",
    "        model.fit(\n",
    "            X_train,\n",
    "            y_train,\n",
    "            eval_set=[validation_data],\n",
    "            eval_metric=\"mse\",\n",
    "            early_stopping_rounds=5,\n",
    "        )  # fit the model\n",
    "        X_val, y_val = validation_data\n",
    "        y_pred = model.predict(\n",
    "            X_val, num_iteration=model.best_iteration_\n",
    "        )  # evaluate the model\n",
    "        eval_mse = mean_squared_error(y_val, y_pred)\n",
    "        self.save_model(trial.trial_id, model)  # save the model to disk\n",
    "        # inform the oracle of the eval result, the result is a dictionary with the metric names as the keys.\n",
    "        return {\"mse\": eval_mse}\n",
    "\n",
    "    def save_model(self, trial_id, model, step=0):\n",
    "        fname = os.path.join(self.get_trial_dir(trial_id), \"model.txt\")\n",
    "        model.booster_.save_model(fname, num_iteration=model.best_iteration_)\n",
    "\n",
    "    def load_model(self, trial):\n",
    "        fname = os.path.join(self.get_trial_dir(trial.trial_id), \"model.txt\")\n",
    "        model = lgb.Booster(model_file=fname)\n",
    "        return model\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Customize Bayesian Optimization search algorithm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "import numpy as np\n",
    "from scipy import optimize as scipy_optimize\n",
    "from scipy.stats import norm\n",
    "from sklearn import exceptions\n",
    "from sklearn import gaussian_process\n",
    "\n",
    "from keras_tuner.engine import hyperparameters as hp_module\n",
    "from keras_tuner.engine import multi_execution_tuner\n",
    "from keras_tuner.engine import oracle as oracle_module\n",
    "from keras_tuner.engine import trial as trial_lib\n",
    "\n",
    "\n",
    "class BayesianOptimizationOracle(oracle_module.Oracle):\n",
    "    \"\"\"Bayesian optimization oracle.\n",
    "\n",
    "    It uses Bayesian optimization with a underlying Gaussian process model.\n",
    "    The acquisition function used is upper confidence bound (UCB), which can\n",
    "    be found in the following link:\n",
    "    https://www.cse.wustl.edu/~garnett/cse515t/spring_2015/files/lecture_notes/12.pdf\n",
    "\n",
    "    # Arguments\n",
    "        objective: String or `kerastuner.Objective`. If a string,\n",
    "          the direction of the optimization (min or max) will be\n",
    "          inferred.\n",
    "        max_trials: Int. Total number of trials\n",
    "            (model configurations) to test at most.\n",
    "            Note that the oracle may interrupt the search\n",
    "            before `max_trial` models have been tested if the search space has been\n",
    "            exhausted.\n",
    "        num_initial_points: (Optional) Int. The number of randomly generated samples\n",
    "            as initial training data for Bayesian optimization. (If not specified,\n",
    "            a trick is to use the square root of the dimensionality of the\n",
    "            hyperparameter space.)\n",
    "        beta: Float. The balancing factor of exploration and exploitation.\n",
    "            The larger it is, the more explorative it is.\n",
    "        seed: Int. Random seed.\n",
    "        hyperparameters: HyperParameters class instance.\n",
    "            Can be used to override (or register in advance)\n",
    "            hyperparamters in the search space.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(\n",
    "        self,\n",
    "        objective,\n",
    "        max_trials,\n",
    "        beta=2.6,\n",
    "        acq_type=\"ucb\",\n",
    "        num_initial_points=None,\n",
    "        seed=None,\n",
    "        hyperparameters=None,\n",
    "        *args,\n",
    "        **kwargs\n",
    "    ):\n",
    "        super(BayesianOptimizationOracle, self).__init__(\n",
    "            objective=objective,\n",
    "            max_trials=max_trials,\n",
    "            hyperparameters=hyperparameters,\n",
    "            seed=seed,\n",
    "            *args,\n",
    "            **kwargs\n",
    "        )\n",
    "        # Use 2 as the initial number of random points if not presented.\n",
    "        self.num_initial_points = num_initial_points or 2\n",
    "        self.beta = beta\n",
    "        self.seed = seed or random.randint(1, 1e4)\n",
    "        self._random_state = np.random.RandomState(self.seed)\n",
    "        self.gpr = self._make_gpr()\n",
    "        self.acq_type = acq_type\n",
    "\n",
    "    def _make_gpr(self):\n",
    "        return gaussian_process.GaussianProcessRegressor(\n",
    "            kernel=gaussian_process.kernels.Matern(nu=2.5),\n",
    "            alpha=1e-4,\n",
    "            normalize_y=True,\n",
    "            random_state=self.seed,\n",
    "        )\n",
    "\n",
    "    def _vectorize_trials(self):\n",
    "        x, y = [], []\n",
    "        for trial in self.trials.values():\n",
    "            # Create a vector representation of each Trial's hyperparameters.\n",
    "            trial_hps = trial.hyperparameters\n",
    "            vector = []\n",
    "            nonfixed_hp_space = [\n",
    "                hp\n",
    "                for hp in self.hyperparameters.space\n",
    "                if not isinstance(hp, hp_module.Fixed)\n",
    "            ]\n",
    "            for hp in nonfixed_hp_space:\n",
    "                # For hyperparameters not present in the trial (either added after\n",
    "                # the trial or inactive in the trial), set to default value.\n",
    "                if trial_hps.is_active(hp):\n",
    "                    trial_value = trial_hps.values[hp.name]\n",
    "                else:\n",
    "                    trial_value = hp.default\n",
    "\n",
    "                # Embed an HP value into the continuous space [0, 1].\n",
    "                prob = hp.value_to_prob(trial_value)\n",
    "                vector.append(prob)\n",
    "\n",
    "            if trial.status == \"COMPLETED\":\n",
    "                score = trial.score\n",
    "                if self.objective.direction == \"min\":\n",
    "                    score = -1 * score\n",
    "            else:\n",
    "                continue\n",
    "\n",
    "            x.append(vector)\n",
    "            y.append(score)\n",
    "\n",
    "        x = np.array(x)\n",
    "        y = np.array(y)\n",
    "        return x, y\n",
    "\n",
    "    def _vector_to_values(self, vector):\n",
    "        hps = hp_module.HyperParameters()\n",
    "        vector_index = 0\n",
    "        for hp in self.hyperparameters.space:\n",
    "            hps.merge([hp])\n",
    "            if isinstance(hp, hp_module.Fixed):\n",
    "                value = hp.value\n",
    "            else:\n",
    "                prob = vector[vector_index]\n",
    "                vector_index += 1\n",
    "                value = hp.prob_to_value(prob)\n",
    "\n",
    "            if hps.is_active(hp):\n",
    "                hps.values[hp.name] = value\n",
    "        return hps.values\n",
    "\n",
    "    def _random_populate_space(self):\n",
    "        values = self._random_values()\n",
    "        if values is None:\n",
    "            return {\"status\": trial_lib.TrialStatus.STOPPED, \"values\": None}\n",
    "        return {\"status\": trial_lib.TrialStatus.RUNNING, \"values\": values}\n",
    "\n",
    "    def _num_completed_trials(self):\n",
    "        return len([t for t in self.trials.values() if t.status == \"COMPLETED\"])\n",
    "\n",
    "    def populate_space(self, trial_id):\n",
    "        if self._num_completed_trials() < self.num_initial_points:\n",
    "            return self._random_populate_space()\n",
    "\n",
    "        # Update Gaussian process regressor\n",
    "        x, y = self._vectorize_trials()\n",
    "        try:\n",
    "            self.gpr.fit(x, y)\n",
    "        except exceptions.ConvergenceWarning as e:\n",
    "            raise e\n",
    "\n",
    "        # Three acquisition functions\n",
    "        def _upper_confidence_bound(x):\n",
    "            x = x.reshape(1, -1)\n",
    "            mu, sigma = self.gpr.predict(x, return_std=True)\n",
    "            return -1 * (mu + self.beta * sigma)\n",
    "\n",
    "        def _probability_of_improvement(x):\n",
    "            # calculate the best surrogate score found so far\n",
    "            x_history, _ = self._vectorize_trials()\n",
    "            y_pred = self.gpr.predict(x_history, return_std=False)\n",
    "            y_best = max(y_pred)\n",
    "            # calculate mean and stdev via surrogate function\n",
    "            x = x.reshape(1, -1)\n",
    "            mu, sigma = self.gpr.predict(x, return_std=True)\n",
    "            # calculate the probability of improvement\n",
    "            z = (mu - y_best) / (sigma + 1e-9)\n",
    "            prob = norm.cdf(z)\n",
    "            return -1 * prob\n",
    "\n",
    "        def _expected_improvement(x):\n",
    "            # calculate the best surrogate score found so far\n",
    "            x_history, _ = self._vectorize_trials()\n",
    "            y_pred = self.gpr.predict(x_history, return_std=False)\n",
    "            y_best = max(y_pred)\n",
    "            # calculate mean and stdev via surrogate function\n",
    "            x = x.reshape(1, -1)\n",
    "            mu, sigma = self.gpr.predict(x, return_std=True)\n",
    "            # calculate the probability of improvement\n",
    "            z = (mu - y_best) / (sigma + 1e-9)\n",
    "            ei = (mu - y_best) * norm.cdf(z) + sigma * norm.pdf(z)\n",
    "            return -1 * ei\n",
    "\n",
    "        acq_funcs = {\n",
    "            \"ucb\": _upper_confidence_bound,\n",
    "            \"pi\": _probability_of_improvement,\n",
    "            \"ei\": _expected_improvement,\n",
    "        }\n",
    "\n",
    "        # Sampling based on acquisition functions\n",
    "        optimal_val = float(\"inf\")\n",
    "        optimal_x = None\n",
    "        num_restarts = 50\n",
    "        bounds = self._get_hp_bounds()\n",
    "        x_seeds = self._random_state.uniform(\n",
    "            bounds[:, 0], bounds[:, 1], size=(num_restarts, bounds.shape[0])\n",
    "        )\n",
    "        for x_try in x_seeds:\n",
    "            # Sign of score is flipped when maximizing.\n",
    "            result = scipy_optimize.minimize(\n",
    "                acq_funcs[self.acq_type], x0=x_try, bounds=bounds, method=\"L-BFGS-B\"\n",
    "            )\n",
    "            result_fun = result.fun if np.isscalar(result.fun) else result.fun[0]\n",
    "            if result_fun < optimal_val:\n",
    "                optimal_val = result_fun\n",
    "                optimal_x = result.x\n",
    "\n",
    "        values = self._vector_to_values(optimal_x)\n",
    "        return {\"status\": trial_lib.TrialStatus.RUNNING, \"values\": values}\n",
    "\n",
    "    def _get_hp_bounds(self):\n",
    "        nonfixed_hp_space = [\n",
    "            hp\n",
    "            for hp in self.hyperparameters.space\n",
    "            if not isinstance(hp, hp_module.Fixed)\n",
    "        ]\n",
    "        bounds = []\n",
    "        for hp in nonfixed_hp_space:\n",
    "            bounds.append([0, 1])\n",
    "        return np.array(bounds)\n",
    "\n",
    "    def get_state(self):\n",
    "        state = super(BayesianOptimizationOracle, self).get_state()\n",
    "        state.update(\n",
    "            {\n",
    "                \"num_initial_points\": self.num_initial_points,\n",
    "                \"acq_type\": self.acq_type,\n",
    "                \"beta\": self.beta,\n",
    "                \"seed\": self.seed,\n",
    "            }\n",
    "        )\n",
    "        return state\n",
    "\n",
    "    def set_state(self, state):\n",
    "        super(BayesianOptimizationOracle, self).set_state(state)\n",
    "        self.num_initial_points = state[\"num_initial_points\"]\n",
    "        self.acq_type = state[\"acq_type\"]\n",
    "        self.beta = state[\"beta\"]\n",
    "        self.seed = state[\"seed\"]\n",
    "        self._random_state = np.random.RandomState(self.seed)\n",
    "        self.gpr = self._make_gpr()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Use customized Bayesian Optimization search algorithm to tune models\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "bo_tuner = LightGBMTuner(\n",
    "    oracle=BayesianOptimizationOracle(\n",
    "        objective=kt.Objective(\"mse\", \"min\"),\n",
    "        max_trials=100,\n",
    "        acq_type=\"ucb\",  # you can switch between different acquisition functions\n",
    "        seed=42,\n",
    "    ),\n",
    "    hypermodel=build_model,\n",
    "    overwrite=True,\n",
    "    project_name=\"bo_tuner\",\n",
    ")\n",
    "\n",
    "bo_tuner.search(X_train, y_train, validation_data=(X_val, y_val))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "best_model = bo_tuner.get_best_models(1)[0]\n",
    "y_pred_test = best_model.predict(X_test)\n",
    "test_mse = mean_squared_error(y_test, y_pred_test)\n",
    "print(\"The prediction MSE on test set: {}\".format(test_mse))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "bo_tuner.results_summary(1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Plot search curves\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def plot_curve(x, y, xlabel, ylabel, title):\n",
    "    plt.plot(x, y)\n",
    "    plt.xlabel(xlabel)\n",
    "    plt.ylabel(ylabel)\n",
    "    plt.title(title)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plot_curves(\n",
    "    x, ys, xlabel, ylabel, title, ymin, ymax, legend, markers, linestyles, markevery=1\n",
    "):\n",
    "    for i, y in enumerate(ys):\n",
    "        plt.plot(x, y, marker=markers[i], linestyle=linestyles[i], markevery=markevery)\n",
    "    plt.xlabel(xlabel)\n",
    "    plt.ylabel(ylabel)\n",
    "    plt.title(title)\n",
    "    plt.ylim(ymin, ymax)\n",
    "    plt.legend(legend)\n",
    "    plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "mse_bo = [\n",
    "    bo_tuner.oracle.get_trial(trial_id).score for trial_id in bo_tuner.oracle.end_order\n",
    "]\n",
    "ids = list(range(len(mse_bo)))\n",
    "plot_curve(\n",
    "    ids, mse_bo, \"Trials in finishing order\", \"Validation MSE\", \"Searched results\"\n",
    ")\n",
    "\n",
    "high_value = float(\"inf\")\n",
    "high_mse_bo = []\n",
    "for value in mse_bo:\n",
    "    high_value = min(high_value, value)\n",
    "    high_mse_bo.append(high_value)\n",
    "plot_curve(\n",
    "    ids,\n",
    "    high_mse_bo,\n",
    "    \"Trials in finishing order\",\n",
    "    \"Highest validation MSE so far\",\n",
    "    \"Searched results\",\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "random_tuner = LightGBMTuner(\n",
    "    oracle=kt.oracles.RandomSearch(\n",
    "        objective=kt.Objective(\"mse\", \"min\"), max_trials=100, seed=42\n",
    "    ),\n",
    "    hypermodel=build_model,\n",
    "    overwrite=True,\n",
    "    project_name=\"random_tuner\",\n",
    ")\n",
    "\n",
    "random_tuner.search(X_train, y_train, validation_data=(X_val, y_val))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "best_model = random_tuner.get_best_models(1)[0]\n",
    "y_pred_test = best_model.predict(X_test)\n",
    "test_mse = mean_squared_error(y_test, y_pred_test)\n",
    "print(\"The prediction MSE on test set: {}\".format(test_mse))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "random_tuner.results_summary(1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "mse_random = [\n",
    "    random_tuner.oracle.get_trial(trial_id).score\n",
    "    for trial_id in random_tuner.oracle.end_order\n",
    "]\n",
    "mse_bo = [\n",
    "    bo_tuner.oracle.get_trial(trial_id).score for trial_id in bo_tuner.oracle.end_order\n",
    "]\n",
    "print(len(mse_random))\n",
    "print(len(mse_bo))\n",
    "\n",
    "high_value = float(\"inf\")\n",
    "high_mse_random = []\n",
    "for value in mse_random:\n",
    "    high_value = min(high_value, value)\n",
    "    high_mse_random.append(high_value)\n",
    "\n",
    "high_value = float(\"inf\")\n",
    "high_mse_bo = []\n",
    "for value in mse_bo:\n",
    "    high_value = min(high_value, value)\n",
    "    high_mse_bo.append(high_value)\n",
    "\n",
    "plot_curves(\n",
    "    ids,\n",
    "    [mse_random, mse_bo],\n",
    "    \"Trials in finishing order\",\n",
    "    \"Validation MSE\",\n",
    "    \"Searched results\",\n",
    "    0,\n",
    "    1.5,\n",
    "    markers=[\"o\", \"+\"],\n",
    "    linestyles=[\"-\", \"-.\"],\n",
    "    legend=[\"Random search\", \"Bayesian optimization\"],\n",
    ")\n",
    "plot_curves(\n",
    "    ids,\n",
    "    [high_mse_random, high_mse_bo],\n",
    "    \"Trials in finishing order\",\n",
    "    \"Highest validation MSE so far\",\n",
    "    \"Searched results\",\n",
    "    0.2,\n",
    "    0.4,\n",
    "    markers=[\"o\", \"+\"],\n",
    "    linestyles=[\"-\", \"-.\"],\n",
    "    legend=[\"Random search\", \"Bayesian optimization\"],\n",
    "    markevery=5,\n",
    ")\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "7.3-Bayesian-Optimization",
   "private_outputs": false,
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}