featsel_benchmark_linearmodels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature selection with vanilla linear models\n",
    "\n",
    "The multi-step feature selection approach of AutoFeat is based on multiple calls to an L1-regularized linear model to select the features in several rounds. For this, the linear model needs to be trained quickly and select reliable features in terms of precision (only relevant features, not too much junk) and recall (all relevant features). In this notebook, several sklearn models are benchmarked to see which model should be used for AutoFeat.\n",
    "\n",
    "#### Regression\n",
    "`ElasticNet` and `Lasso` are very similar, so are `Lars` and `LassoLars`. `LassoLarsIC` sometimes does not select any features. `OrthogonalMatchingPursuit` is very fast but also often selects too few features. With a good trade-off in terms of speed and selected features, we use the `LassoLarsCV` model for AutoFeat.\n",
    "\n",
    "#### Classification\n",
    "We tested `linear_model.LogisticRegressionCV` and `svm.LinearSVC` (together with a grid search), however, Logistic Regression is much faster, therefore it does not make sense to use the SVC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import colorsys\n",
    "import warnings\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sklearn.linear_model as lm\n",
    "from time import time\n",
    "from sklearn import svm\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.datasets import make_classification\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "from autofeat import AutoFeatModel\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_colors(n=100):\n",
    "    HSV_tuples = [(x * 1.0 / (n + 1), 1.0, 0.8) for x in range(n)]\n",
    "    return [colorsys.hsv_to_rgb(*x) for x in HSV_tuples]\n",
    "\n",
    "\n",
    "def create_plots(precision, recall, noise_levels, noise_feat_frac, n_train, n_feat_true):\n",
    "    colors = get_colors(len(noise_levels))\n",
    "    plt.figure()\n",
    "    for i, noise in enumerate(noise_levels):\n",
    "        plt.plot(noise_feat_frac, precision[i, :], c=colors[i], label=\"noise: %g\" % noise)\n",
    "    plt.legend(bbox_to_anchor=(1.0, 1.0))\n",
    "    plt.xlabel(\"# noise features / # training samples ($n$)\")\n",
    "    plt.title(\"Precision ($n$: %i; $d$: %i)\" % (n_train, n_feat_true))\n",
    "    plt.figure()\n",
    "    for i, noise in enumerate(noise_levels):\n",
    "        plt.plot(noise_feat_frac, recall[i, :], c=colors[i], label=\"noise: %g\" % noise)\n",
    "    plt.legend(bbox_to_anchor=(1.0, 1.0))\n",
    "    plt.xlabel(\"# noise features / # training samples ($n$)\")\n",
    "    plt.title(\"Recall ($n$: %i; $d$: %i)\" % (n_train, n_feat_true))\n",
    "    colors = get_colors(len(noise_feat_frac))\n",
    "    plt.figure()\n",
    "    for i, nfeat in enumerate(noise_feat_frac):\n",
    "        plt.plot(noise_levels, precision[:, i], c=colors[i], label=\"noise feat frac: %g\" % nfeat)\n",
    "    plt.legend(bbox_to_anchor=(1.0, 1.0))\n",
    "    plt.xlabel(\"noise level\")\n",
    "    plt.title(\"Precision ($n$: %i; $d$: %i)\" % (n_train, n_feat_true))\n",
    "    plt.figure()\n",
    "    for i, nfeat in enumerate(noise_feat_frac):\n",
    "        plt.plot(noise_levels, recall[:, i], c=colors[i], label=\"noise feat frac: %g\" % nfeat)\n",
    "    plt.legend(bbox_to_anchor=(1.0, 1.0))\n",
    "    plt.xlabel(\"noise level\")\n",
    "    plt.title(\"Recall ($n$: %i; $d$: %i)\" % (n_train, n_feat_true))\n",
    "    plt.figure()\n",
    "    plt.imshow(precision)\n",
    "    plt.xticks(list(range(len(noise_feat_frac))), noise_feat_frac)\n",
    "    plt.xlabel(\"# noise features / # training samples ($n$)\")\n",
    "    plt.yticks(list(range(len(noise_levels))), noise_levels)\n",
    "    plt.ylabel(\"noise level\")\n",
    "    plt.clim(0, 1)\n",
    "    plt.colorbar()\n",
    "    plt.title(\"Precision ($n$: %i; $d$: %i)\" % (n_train, n_feat_true))\n",
    "    plt.figure()\n",
    "    plt.imshow(recall)\n",
    "    plt.xticks(list(range(len(noise_feat_frac))), noise_feat_frac)\n",
    "    plt.xlabel(\"# noise features / # training samples ($n$)\")\n",
    "    plt.yticks(list(range(len(noise_levels))), noise_levels)\n",
    "    plt.ylabel(\"noise level\")\n",
    "    plt.clim(0, 1)\n",
    "    plt.colorbar()\n",
    "    plt.title(\"Recall ($n$: %i; $d$: %i)\" % (n_train, n_feat_true))\n",
    "\n",
    "\n",
    "def prec_rec(coefs, true_features):\n",
    "    # sort weights by absolute values\n",
    "    w_dict = dict(zip(range(len(coefs)), np.abs(coefs)))\n",
    "    sorted_features = sorted(w_dict, key=w_dict.get, reverse=True)\n",
    "    # --> how many of the true features are amongst the len(true_features) highest weights?\n",
    "    n_selected = int(min(len(true_features), np.sum(np.abs(coefs) > 0)))\n",
    "    recall = len(set(sorted_features[:n_selected]).intersection(true_features)) / len(true_features)\n",
    "    # --> how many features do we need to take to get all true features?\n",
    "    min_thr = min(w_dict[t] for t in true_features)\n",
    "    precision = len(true_features) / np.sum(np.abs(coefs) >= min_thr)\n",
    "    return precision, recall\n",
    "\n",
    "\n",
    "def prec_rec_autofeat(selected_features, true_features):\n",
    "    if not len(selected_features):\n",
    "        return 0, 0\n",
    "    TP = len(set(selected_features).intersection(true_features))\n",
    "    recall = TP / len(true_features)\n",
    "    precision = TP / len(selected_features)\n",
    "    return precision, recall\n",
    "\n",
    "\n",
    "def get_dataset(n_train=100, n_feat_noise=100, n_feat_true=10, noise=0.01, ptype=\"regression\", random_seed=None):\n",
    "    if random_seed is not None:\n",
    "        np.random.seed(random_seed)\n",
    "    # create a simple problem with 1 target variable,\n",
    "    # generated as a linear combination of random features\n",
    "    # and include some additional noise filters\n",
    "    X = np.random.randn(n_train, n_feat_noise + n_feat_true)\n",
    "    true_features = np.random.permutation(n_feat_noise + n_feat_true)[:n_feat_true]\n",
    "    y = X[:, true_features].sum(axis=1)\n",
    "    if ptype == \"regression\":\n",
    "        # add normally distributed noise\n",
    "        y -= y.mean()\n",
    "        y /= y.std()\n",
    "        y = (1 - noise) * y + noise * np.random.randn(len(y))\n",
    "        # y = y/y.std() + noise*np.random.randn(len(y))\n",
    "    else:\n",
    "        # threshold to transform into classification problem\n",
    "        y = np.array(y > y.mean(), dtype=int)\n",
    "        # randomly flip some idx\n",
    "        flip_idx = np.random.permutation(len(y))[: int(np.ceil(len(y) * noise))]\n",
    "        y[flip_idx] -= 1\n",
    "        y = np.abs(y)\n",
    "    return X, y, true_features\n",
    "\n",
    "\n",
    "def compute_autofeat_dataset(n_train=100):\n",
    "    np.random.seed(10)\n",
    "    x1 = np.random.rand(n_train)\n",
    "    x2 = np.random.randn(n_train)\n",
    "    x3 = np.random.rand(n_train)\n",
    "    y = 2 + 15 * x1 + 3 / (x2 - 1 / x3) + 5 * (x2 + np.log(x1)) ** 3\n",
    "    X = pd.DataFrame(np.vstack([x1, x2, x3]).T, columns=[\"x1\", \"x2\", \"x3\"])\n",
    "    # generate new features with autofeat\n",
    "    afreg = AutoFeatModel(verbose=1, feateng_steps=3, featsel_runs=-1, problem_type=None)\n",
    "    X = afreg.fit_transform(X, y)\n",
    "    return X, y\n",
    "\n",
    "\n",
    "def get_autofeat_dataset(n_feat_noise=100, noise=0.01, ptype=\"regression\", random_seed=None):\n",
    "    if random_seed is not None:\n",
    "        np.random.seed(random_seed)\n",
    "    true_features = [\"x1\", \"(x2 + log(x1))**3\", \"1/(x2 - 1/x3)\"]\n",
    "    noise_features = list(np.random.permutation(list(set(df.columns).difference(true_features)))[:n_feat_noise])\n",
    "    X = StandardScaler().fit_transform(df[true_features + noise_features])\n",
    "    y = target - target.mean()\n",
    "    y /= y.std()\n",
    "    if ptype == \"regression\":\n",
    "        # add normally distributed noise\n",
    "        y = (1 - noise) * y + noise * np.random.randn(len(y))\n",
    "    else:\n",
    "        # threshold to transform into classification problem\n",
    "        y = np.array(y > y.mean(), dtype=int)\n",
    "        # randomly flip some idx\n",
    "        flip_idx = np.random.permutation(len(y))[: int(np.ceil(len(y) * noise))]\n",
    "        y[flip_idx] -= 1\n",
    "        y = np.abs(y)\n",
    "    return X, y, [0, 1, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the original autofeat dataset with all features\n",
    "df, target = compute_autofeat_dataset(n_train=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# independent features\n",
    "ptype = \"regression\"\n",
    "t0 = time()\n",
    "n_train = 100\n",
    "n_feat_true = 10\n",
    "noise_levels = [0.0, 0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1.0]\n",
    "noise_feat_frac = [0.01, 0.05, 0.1, 0.25, 0.5, 1, 2, 5, 10, 25]\n",
    "precision = np.zeros((len(noise_levels), len(noise_feat_frac)))\n",
    "recall = np.zeros((len(noise_levels), len(noise_feat_frac)))\n",
    "for i, noise in enumerate(noise_levels):\n",
    "    print(\"noise level: %g\" % noise)\n",
    "    for j, nfeat in enumerate(noise_feat_frac):\n",
    "        ps, rs = [], []\n",
    "        for seed in range(10):\n",
    "            X, y, true_features = get_dataset(n_train, int(nfeat * n_train), n_feat_true, noise, ptype, random_seed=seed)\n",
    "            with warnings.catch_warnings():\n",
    "                warnings.simplefilter(\"ignore\")\n",
    "                if ptype == \"regression\":\n",
    "                    # model = lm.OrthogonalMatchingPursuitCV(cv=5, max_iter=min(X.shape[1], 20)).fit(X, y)\n",
    "                    model = lm.LassoLarsCV(cv=5).fit(X, y)\n",
    "                    coefs = model.coef_\n",
    "                else:\n",
    "                    model = lm.LogisticRegressionCV(cv=5, penalty=\"l1\", solver=\"saga\", class_weight=\"balanced\").fit(X, y)\n",
    "                    coefs = np.max(np.abs(model.coef_), axis=0)\n",
    "            p, r = prec_rec(coefs, true_features)\n",
    "            ps.append(p)\n",
    "            rs.append(r)\n",
    "        precision[i, j] = np.mean(ps)\n",
    "        recall[i, j] = np.mean(rs)\n",
    "print(\"took %.1f sec\" % (time() - t0))\n",
    "\n",
    "create_plots(precision, recall, noise_levels, noise_feat_frac, n_train, n_feat_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# autofeat features\n",
    "ptype = \"regression\"\n",
    "t0 = time()\n",
    "n_train = 100\n",
    "n_feat_true = 3\n",
    "noise_levels = [0.0, 0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1.0]\n",
    "noise_feat_frac = [0.01, 0.05, 0.1, 0.25, 0.5, 1, 2, 5, 10, 25]\n",
    "precision = np.zeros((len(noise_levels), len(noise_feat_frac)))\n",
    "recall = np.zeros((len(noise_levels), len(noise_feat_frac)))\n",
    "for i, noise in enumerate(noise_levels):\n",
    "    print(\"noise level: %g\" % noise)\n",
    "    for j, nfeat in enumerate(noise_feat_frac):\n",
    "        ps, rs = [], []\n",
    "        for seed in range(10):\n",
    "            X, y, true_features = get_autofeat_dataset(int(nfeat * n_train), noise, ptype, random_seed=seed)\n",
    "            with warnings.catch_warnings():\n",
    "                warnings.simplefilter(\"ignore\")\n",
    "                if ptype == \"regression\":\n",
    "                    # model = lm.OrthogonalMatchingPursuitCV(cv=5, max_iter=min(X.shape[1], 20)).fit(X, y)\n",
    "                    model = lm.LassoLarsCV(cv=5, eps=1e-5).fit(X, y)\n",
    "                    coefs = model.coef_\n",
    "                else:\n",
    "                    model = lm.LogisticRegressionCV(cv=5, penalty=\"l1\", solver=\"saga\", class_weight=\"balanced\").fit(X, y)\n",
    "                    coefs = np.max(np.abs(model.coef_), axis=0)\n",
    "                selected_features = np.where(coefs > 1e-8)[0]\n",
    "            p, r = prec_rec_autofeat(selected_features, true_features)\n",
    "            ps.append(p)\n",
    "            rs.append(r)\n",
    "        precision[i, j] = np.mean(ps)\n",
    "        recall[i, j] = np.mean(rs)\n",
    "print(\"took %.1f sec\" % (time() - t0))\n",
    "\n",
    "create_plots(precision, recall, noise_levels, noise_feat_frac, n_train, n_feat_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# multi-class classification\n",
    "X, y = make_classification(2000, 5000, n_informative=20, n_redundant=100, n_repeated=0, n_classes=5, random_state=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter(\"ignore\")\n",
    "    model = lm.LogisticRegressionCV(cv=5, penalty=\"l1\", solver=\"saga\", class_weight=\"balanced\").fit(X, y)\n",
    "    print(model.score(X, y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter(\"ignore\")\n",
    "    param_grid = {\"C\": np.logspace(-4, 4, 10)}\n",
    "    clf = svm.LinearSVC(penalty=\"l1\", class_weight=\"balanced\", dual=False)\n",
    "    model = GridSearchCV(clf, param_grid, cv=5)\n",
    "    model.fit(X, y)\n",
    "    print(model.score(X, y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}