notebooks/linear_regression_non_linear_link.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear regression for a non-linear features-target relationship\n",
    "\n",
    "In the previous exercise, you were asked to train a linear regression model\n",
    "on a dataset where the matrix `data` and the vector `target` do not have a\n",
    "linear link.\n",
    "\n",
    "In this notebook, we show that even if the parametrization of linear models\n",
    "is not natively adapted to the problem at hand, it is still possible to make\n",
    "linear models more expressive by engineering additional features.\n",
    "\n",
    "A machine learning pipeline that combines a non-linear feature engineering\n",
    "step followed by a linear regression step can therefore be considered\n",
    "non-linear regression model as a whole.\n",
    "\n",
    "To illustrate these concepts, we will reuse the same dataset generated in the\n",
    "previous exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "rng = np.random.RandomState(0)\n",
    "\n",
    "n_sample = 100\n",
    "data_max, data_min = 1.4, -1.4\n",
    "len_data = (data_max - data_min)\n",
    "# sort the data to make plotting easier later\n",
    "data = np.sort(rng.rand(n_sample) * len_data - len_data / 2)\n",
    "noise = rng.randn(n_sample) * .3\n",
    "target = data ** 3 - 0.5 * data ** 2 + noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition note alert alert-info\">\n",
    "<p class=\"first admonition-title\" style=\"font-weight: bold;\">Note</p>\n",
    "<p class=\"last\">To ease the plotting, we will create a pandas dataframe containing the data\n",
    "and target:</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "full_data = pd.DataFrame({\"input_feature\": data, \"target\": target})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "_ = sns.scatterplot(data=full_data, x=\"input_feature\", y=\"target\",\n",
    "                    color=\"black\", alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will highlight the limitations of fitting a linear regression model as\n",
    "done in the previous exercise.\n",
    "\n",
    "<div class=\"admonition warning alert alert-danger\">\n",
    "<p class=\"first admonition-title\" style=\"font-weight: bold;\">Warning</p>\n",
    "<p class=\"last\">In scikit-learn, by convention <tt class=\"docutils literal\">data</tt> (also called <tt class=\"docutils literal\">X</tt> in the scikit-learn\n",
    "documentation) should be a 2D matrix of shape <tt class=\"docutils literal\">(n_samples, n_features)</tt>.\n",
    "If <tt class=\"docutils literal\">data</tt> is a 1D vector, you need to reshape it into a matrix with a\n",
    "single column if the vector represents a feature or a single row if the\n",
    "vector represents a sample.</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# X should be 2D for sklearn: (n_samples, n_features)\n",
    "data = data.reshape((-1, 1))\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "linear_regression = LinearRegression()\n",
    "linear_regression.fit(data, target)\n",
    "target_predicted = linear_regression.predict(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mse = mean_squared_error(target, target_predicted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.scatterplot(data=full_data, x=\"input_feature\", y=\"target\",\n",
    "                     color=\"black\", alpha=0.5)\n",
    "ax.plot(data, target_predicted)\n",
    "_ = ax.set_title(f\"Mean squared error = {mse:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Here the coefficient and intercept learnt by `LinearRegression` define the\n",
    "best \"straight line\" that fits the data. We can inspect the coefficients\n",
    "using the attributes of the model learnt as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"weight: {linear_regression.coef_[0]:.2f}, \"\n",
    "      f\"intercept: {linear_regression.intercept_:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "It is important to note that the learnt model will not be able to handle the\n",
    "non-linear relationship between `data` and `target` since linear models\n",
    "assume the relationship between `data` and `target` to be linear.\n",
    "\n",
    "Indeed, there are 3 possibilities to solve this issue:\n",
    "\n",
    "1. choose a model that can natively deal with non-linearity,\n",
    "2. engineer a richer set of features by including expert knowledge which can\n",
    "   be directly used by a simple linear model, or\n",
    "3. use a \"kernel\" to have a locally-based decision function instead of a\n",
    "   global linear decision function.\n",
    "\n",
    "Let's illustrate quickly the first point by using a decision tree regressor\n",
    "which can natively handle non-linearity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.tree import DecisionTreeRegressor\n",
    "\n",
    "tree = DecisionTreeRegressor(max_depth=3).fit(data, target)\n",
    "target_predicted = tree.predict(data)\n",
    "mse = mean_squared_error(target, target_predicted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.scatterplot(data=full_data, x=\"input_feature\", y=\"target\",\n",
    "                     color=\"black\", alpha=0.5)\n",
    "ax.plot(data, target_predicted)\n",
    "_ = ax.set_title(f\"Mean squared error = {mse:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Instead of having a model which can natively deal with non-linearity, we\n",
    "could also modify our data: we could create new features, derived from the\n",
    "original features, using some expert knowledge. In this example, we know that\n",
    "we have a cubic and squared relationship between `data` and `target` (because\n",
    "we generated the data).\n",
    "\n",
    "Indeed, we could create two new features (`data ** 2` and `data ** 3`) using\n",
    "this information as follows. This kind of transformation is called a\n",
    "polynomial feature expansion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "data_expanded = np.concatenate([data, data ** 2, data ** 3], axis=1)\n",
    "data_expanded.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linear_regression.fit(data_expanded, target)\n",
    "target_predicted = linear_regression.predict(data_expanded)\n",
    "mse = mean_squared_error(target, target_predicted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.scatterplot(data=full_data, x=\"input_feature\", y=\"target\",\n",
    "                     color=\"black\", alpha=0.5)\n",
    "ax.plot(data, target_predicted)\n",
    "_ = ax.set_title(f\"Mean squared error = {mse:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We can see that even with a linear model, we can overcome the linearity\n",
    "limitation of the model by adding the non-linear components in the design of\n",
    "additional features. Here, we created new features by knowing the way the\n",
    "target was generated.\n",
    "\n",
    "Instead of manually creating such polynomial features one could directly use\n",
    "[sklearn.preprocessing.PolynomialFeatures](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html).\n",
    "\n",
    "To demonstrate the use of the `PolynomialFeatures` class, we use a\n",
    "scikit-learn pipeline which first transforms the features and then fit the\n",
    "regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "polynomial_regression = make_pipeline(\n",
    "    PolynomialFeatures(degree=3, include_bias=False),\n",
    "    LinearRegression(),\n",
    ")\n",
    "polynomial_regression.fit(data, target)\n",
    "target_predicted = polynomial_regression.predict(data)\n",
    "mse = mean_squared_error(target, target_predicted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "In the previous cell we had to set `include_bias=False` as otherwise we would\n",
    "create a column perfectly correlated to the `intercept_` introduced by the\n",
    "`LinearRegression`. We can verify that this procedure is equivalent to creating\n",
    "the features by hand up to numerical error by computing the maximum\n",
    "of the absolute values of the differences between the features generated by\n",
    "both methods and checking that it is close to zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.abs(polynomial_regression[0].fit_transform(data) - data_expanded).max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Then it should not be surprising that the predictions of the\n",
    "`PolynomialFeatures` pipeline match the predictions of the linear model fit on\n",
    "manually engineered features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.scatterplot(data=full_data, x=\"input_feature\", y=\"target\",\n",
    "                     color=\"black\", alpha=0.5)\n",
    "ax.plot(data, target_predicted)\n",
    "_ = ax.set_title(f\"Mean squared error = {mse:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The last possibility is to make a linear model more expressive is to use a\n",
    "\"kernel\". Instead of learning a weight per feature as we previously\n",
    "emphasized, a weight will be assigned to each sample. However, not all\n",
    "samples will be used. This is the base of the support vector machine\n",
    "algorithm.\n",
    "\n",
    "The mathematical definition of \"kernels\" and \"support vector machines\" is\n",
    "beyond the scope of this course. We encourage interested readers with a\n",
    "mathematical training to have a look at the scikit-learn [documentation on\n",
    "SVMs](https://scikit-learn.org/stable/modules/svm.html) for more details.\n",
    "\n",
    "For the rest of us, let us just develop some intuitions on the relative\n",
    "expressive power of support vector machines with linear and non-linear\n",
    "kernels by fitting them on the same dataset.\n",
    "\n",
    "First, consider a support vector machine with a linear kernel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVR\n",
    "\n",
    "svr = SVR(kernel=\"linear\")\n",
    "svr.fit(data, target)\n",
    "target_predicted = svr.predict(data)\n",
    "mse = mean_squared_error(target, target_predicted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.scatterplot(data=full_data, x=\"input_feature\", y=\"target\",\n",
    "                     color=\"black\", alpha=0.5)\n",
    "ax.plot(data, target_predicted)\n",
    "_ = ax.set_title(f\"Mean squared error = {mse:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The predictions of our SVR with a linear kernel are all aligned on a straight\n",
    "line. `SVR(kernel=\"linear\")` is indeed yet another example of a linear model.\n",
    "\n",
    "The estimator can also be configured to use a non-linear kernel. Then, it can\n",
    "learn a prediction function that computes non-linear interaction between\n",
    "samples for which we want to make a prediction and selected samples from the\n",
    "training set.\n",
    "\n",
    "The result is another kind of non-linear regression model with a similar\n",
    "expressivity as our previous polynomial regression pipeline:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svr = SVR(kernel=\"poly\", degree=3)\n",
    "svr.fit(data, target)\n",
    "target_predicted = svr.predict(data)\n",
    "mse = mean_squared_error(target, target_predicted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.scatterplot(data=full_data, x=\"input_feature\", y=\"target\",\n",
    "                     color=\"black\", alpha=0.5)\n",
    "ax.plot(data, target_predicted)\n",
    "_ = ax.set_title(f\"Mean squared error = {mse:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Kernel methods such as SVR are very efficient for small to medium datasets.\n",
    "\n",
    "For larger datasets with `n_samples >> 10_000`, it is often computationally\n",
    "more efficient to perform explicit feature expansion using\n",
    "`PolynomialFeatures` or other non-linear transformers from scikit-learn such\n",
    "as\n",
    "[KBinsDiscretizer](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.KBinsDiscretizer.html)\n",
    "or\n",
    "[Nystroem](https://scikit-learn.org/stable/modules/generated/sklearn.kernel_approximation.Nystroem.html).\n",
    "\n",
    "Here again we refer the interested reader to the documentation to get a\n",
    "proper definition of those methods. The following just gives an intuitive\n",
    "overview of the predictions we would get using those on our toy dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import KBinsDiscretizer\n",
    "\n",
    "binned_regression = make_pipeline(\n",
    "    KBinsDiscretizer(n_bins=8), LinearRegression(),\n",
    ")\n",
    "binned_regression.fit(data, target)\n",
    "target_predicted = binned_regression.predict(data)\n",
    "mse = mean_squared_error(target, target_predicted)\n",
    "\n",
    "ax = sns.scatterplot(data=full_data, x=\"input_feature\", y=\"target\",\n",
    "                     color=\"black\", alpha=0.5)\n",
    "ax.plot(data, target_predicted)\n",
    "_ = ax.set_title(f\"Mean squared error = {mse:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.kernel_approximation import Nystroem\n",
    "\n",
    "nystroem_regression = make_pipeline(\n",
    "    Nystroem(n_components=5), LinearRegression(),\n",
    ")\n",
    "nystroem_regression.fit(data, target)\n",
    "target_predicted = nystroem_regression.predict(data)\n",
    "mse = mean_squared_error(target, target_predicted)\n",
    "\n",
    "ax = sns.scatterplot(data=full_data, x=\"input_feature\", y=\"target\",\n",
    "                     color=\"black\", alpha=0.5)\n",
    "ax.plot(data, target_predicted)\n",
    "_ = ax.set_title(f\"Mean squared error = {mse:.2f}\")"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "main_language": "python"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}