intermediate_notebooks/examples/ridge_regression_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ridge Regression\n",
    "\n",
    "This notebook compares a CPU implementation and a GPU implementation of Ridge Regression.  It includes code example for doing Ridge egression using RAPIDS cuDF and cuML."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook Credits\n",
    "### Authorship\n",
    "Original Author: Unknown  \n",
    "Last Edit: Taurean Dyer, 2/20/2019  \n",
    "\n",
    "### Test System Specs\n",
    "Test System Hardware: DGX-2  \n",
    "Test System Software: Ubuntu 16.04  \n",
    "RAPIDS Version: 0.5.1 - Docker Install  \n",
    "Driver: 410.79  \n",
    "CUDA: 10.0  \n",
    "\n",
    "### Known Working Systems\n",
    "RAPIDS Versions: 0.4, 0.5, 0.5.1  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's Begin: Linear Regression\n",
    "### Imports\n",
    "Let's start with our Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn import linear_model as sklGLM\n",
    "from cuml import LinearRegression as cumlOLS\n",
    "from cuml import Ridge as cumlRidge\n",
    "import cudf\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Helper Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from timeit import default_timer\n",
    "\n",
    "class Timer(object):\n",
    "    def __init__(self):\n",
    "        self._timer = default_timer\n",
    "    \n",
    "    def __enter__(self):\n",
    "        self.start()\n",
    "        return self\n",
    "\n",
    "    def __exit__(self, *args):\n",
    "        self.stop()\n",
    "\n",
    "    def start(self):\n",
    "        \"\"\"Start the timer.\"\"\"\n",
    "        self.start = self._timer()\n",
    "\n",
    "    def stop(self):\n",
    "        \"\"\"Stop the timer. Calculate the interval in seconds.\"\"\"\n",
    "        self.end = self._timer()\n",
    "        self.interval = self.end - self.start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gzip\n",
    "def load_data(nrows, ncols, cached = '../../data/mortgage/mortgage.npy.gz'):\n",
    "    if os.path.exists(cached):\n",
    "        print('use mortgage data')\n",
    "        with gzip.open(cached) as f:\n",
    "            X = np.load(f)\n",
    "        # the 4th column is 'adj_remaining_months_to_maturity'\n",
    "        # used as the label\n",
    "        X = X[:,[i for i in range(X.shape[1]) if i!=4]]\n",
    "        y = X[:,4:5]\n",
    "        rindices = np.random.randint(0,X.shape[0]-1,nrows)\n",
    "        X = X[rindices,:ncols]\n",
    "        y = y[rindices]\n",
    "    else:\n",
    "        print('use random data')\n",
    "        X = np.random.rand(nrows,ncols)\n",
    "        \n",
    "    df_X = pd.DataFrame({'fea%d'%i:X[:,i] for i in range(X.shape[1])})\n",
    "    df_y = pd.DataFrame({'fea%d'%i:y[:,i] for i in range(y.shape[1])})\n",
    "    \n",
    "    return df_X, df_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "def array_equal(a,b,threshold=2e-3,with_sign=True):\n",
    "    a = to_nparray(a).ravel()\n",
    "    b = to_nparray(b).ravel()\n",
    "    if with_sign == False:\n",
    "        a,b = np.abs(a),np.abs(b)\n",
    "    error = mean_squared_error(a,b)\n",
    "    res = error<threshold\n",
    "    return res\n",
    "\n",
    "def to_nparray(x):\n",
    "    if isinstance(x,np.ndarray) or isinstance(x,pd.DataFrame):\n",
    "        return np.array(x)\n",
    "    elif isinstance(x,np.float64):\n",
    "        return np.array([x])\n",
    "    elif isinstance(x,cudf.DataFrame) or isinstance(x,cudf.Series):\n",
    "        return x.to_pandas().values\n",
    "    return x    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our Helper functions, lets start to compare the speed and results for SciKit Learn's CPU impletmenation versus RAPIDS cuML GPU impletementation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "use mortgage data\n",
      "data (1258291, 399)\n",
      "label (1258291, 1)\n",
      "CPU times: user 14.3 s, sys: 1.08 s, total: 15.4 s\n",
      "Wall time: 13.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "nrows = int((2**20) * 1.2)\n",
    "ncols = 399\n",
    "\n",
    "X, y = load_data(nrows,ncols)\n",
    "print('data',X.shape)\n",
    "print('label',y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though the ridge regression interface of cuML is very similar to Scikit-Learn's implemetation, cuML doesn't use some of the parameters such as \"copy\". Also, cuML includes two different implementation of ridge regression using SVD and Eigen decomposition. Eigen decomposition based implementation is very fast but causes very small errors in the coefficients which is negligible for most of the applications. SVD is stable but slower than eigen decomposition based implementation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_intercept = True\n",
    "normalize = False\n",
    "alpha = np.array([1.0])\n",
    "# eig: eigen decomposition based method, \n",
    "# svd: singular value decomposition based method,\n",
    "# cd: coordinate descend.\n",
    "solver = \"eig\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get MSE for SciKit Learn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 15.1 s, sys: 968 ms, total: 16 s\n",
      "Wall time: 4.43 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "reg_sk = sklGLM.Ridge(alpha=alpha, fit_intercept=fit_intercept, normalize=normalize, solver='auto')\n",
    "result_sk = reg_sk.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.48 s, sys: 340 ms, total: 1.82 s\n",
      "Wall time: 991 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "y_sk = reg_sk.predict(X)\n",
    "error_sk = mean_squared_error(y,y_sk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get MSE for cuML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.3 s, sys: 816 ms, total: 3.12 s\n",
      "Wall time: 2.13 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "X_cudf = cudf.DataFrame.from_pandas(X)\n",
    "y_cudf = y.values\n",
    "y_cudf = y_cudf[:,0]\n",
    "y_cudf = cudf.Series(y_cudf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 788 ms, sys: 212 ms, total: 1e+03 ms\n",
      "Wall time: 995 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "reg_cuml = cumlRidge(alpha=alpha, fit_intercept=fit_intercept, normalize=normalize, solver=solver)\n",
    "result_cuml = reg_cuml.fit(X_cudf, y_cudf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 276 ms, sys: 20 ms, total: 296 ms\n",
      "Wall time: 285 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "y_cuml = reg_cuml.predict(X_cudf)\n",
    "y_cuml = to_nparray(y_cuml).ravel()\n",
    "error_cuml = mean_squared_error(y,y_cuml)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Comparison Between SKL and cuML\n",
    "Your final output should have both MSE results close to 0 (about 1.0e-7 to 1.0e-14).  However, despite having similar answers, you should see a **massive reduction to the sys time** when using **RAPIDS cuML** versus **SciKit Learn**.  Go RAPIDS!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SKL MSE(y):\n",
      "7.84409492757235e-12\n",
      "CUML MSE(y):\n",
      "1.3055482e-10\n"
     ]
    }
   ],
   "source": [
    "print(\"SKL MSE(y):\")\n",
    "print(error_sk)\n",
    "print(\"CUML MSE(y):\")\n",
    "print(error_cuml)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}