cpu_comparisons/PCA_Demo_Full.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimensionality Reduction using GPU Accelerated PCA in RAPIDS\n",
    "#### Ivana Jovicic\n",
    "-------\n",
    "\n",
    "While the world’s data doubles each year, CPU computing has hit a brick wall with the end of Moore’s law. For the same reasons, scientific computing and deep learning has turned to NVIDIA GPU acceleration, data analytics and machine learning where GPU acceleration is ideal. \n",
    "\n",
    "NVIDIA created RAPIDS – an open-source data analytics and machine learning acceleration platform that leverages GPUs to accelerate computations. RAPIDS is based on Python, has pandas-like and Scikit-Learn-like interfaces, is built on Apache Arrow in-memory data format, and can scale from 1 to multi-GPU to multi-nodes. RAPIDS integrates easily into the world’s most popular data science Python-based workflows. RAPIDS accelerates data science end-to-end – from data prep, to machine learning, to deep learning. And through Arrow, Spark users can easily move data into the RAPIDS platform for acceleration.\n",
    "\n",
    "In this notebook, we will also show how to use PCA - a popular dimensionality reduction algorithm - and how to use the GPU accelerated implementation of this algorithm in RAPIDS.\n",
    "\n",
    "**Table of Contents**\n",
    "\n",
    "* Clustering with PCA\n",
    "* Why do we use PCA?\n",
    "* Setup\n",
    "* Basic Example using Iris Dataset\n",
    "  * Importing Libraries\n",
    "  * Import Iris Dataset\n",
    "  * Dataset Overview\n",
    "  * Basic Cleaning of the Data\n",
    "  * Applying PCA\n",
    "  * Training and Making Predictions\n",
    "* PCA on Larger Dataset and CPU vs GPU comparisson\n",
    "  * Benchmark with Random Data\n",
    "* Conclusion\n",
    "\n",
    "Before going any further, let's make sure we have access to `matplotlib`, a popular Python library for data visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "try:\n",
    "    import matplotlib\n",
    "except ModuleNotFoundError:\n",
    "    os.system('conda install -y matplotlib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dimensionality Reductino with PCA\n",
    "\n",
    "*Principal component analysis (PCA)* is a technique used for identification of a smaller number of uncorrelated variables known as principal components from a larger set of data. The technique is widely used to emphasize variation and capture strong patterns in a data set. \n",
    "\n",
    "Principal component analysis is considered a useful statistical method and used in fields such as image compression, face recognition, neuroscience and computer graphics.\n",
    "\n",
    "source : https://www.techopedia.com/definition/32509/principal-component-analysis-pca"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why do we use PCA ?\n",
    "\n",
    "Large number of features in the dataset is one of the factors that affect both the training time as well as accuracy of machine learning models. You have different **options to deal with huge number of features in a dataset**.\n",
    "\n",
    "1. Try to train the models on original number of features, which take days or weeks if the number of features is too high.\n",
    "2. Reduce the number of variables by merging correlated variables.\n",
    "3. Extract the most important features from the dataset that are responsible for maximum variance in the output. Different statistical techniques are used for this purpose e.g. linear discriminant analysis, factor analysis, and principal component analysis.\n",
    "\n",
    "Principal component analysis, or PCA, is a statistical technique to **convert high dimensional data to low dimensional data** by selecting the most important features that capture maximum information about the dataset. \n",
    "\n",
    "The features are selected on the basis of variance that they cause in the output. The feature that causes highest variance is the *first principal component*. The feature that is responsible for second highest variance is considered the *second principal component*, and so on. It is important to mention that principal components do not have any correlation with each other.\n",
    "\n",
    "source: https://stackabuse.com/implementing-pca-in-python-with-scikit-learn/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook was tested using the `nvcr.io/nvidia/rapidsai/rapidsai:0.5-cuda10.0-runtime-ubuntu18.04-gcc7-py3.7` Docker container from [NVIDIA GPU Cloud](https://ngc.nvidia.com) and run on the NVIDIA Tesla V100 GPU. Please be aware that your system may be different and you may need to modify the code or install packages to run the below examples. \n",
    "\n",
    "If you think you have found a bug or an error, please file an issue here: https://github.com/rapidsai/notebooks/issues\n",
    "\n",
    "Let's check out our hardware setup by runing the `nvidia-smi` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!nvidia-smi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's see what CUDA version we have."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!nvcc --version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Example using Iris Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing  libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first, we are going to import basic python libraries needed to run a very simple exaple to explain the above mentioned algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.decomposition import PCA as skPCA\n",
    "\n",
    "\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import Iris Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://archive.ics.uci.edu/ml/datasets/iris\n",
    "\n",
    "The dataset consists of 150 records of Iris plant with **four features**: 'sepal-length', 'sepal-width', 'petal-length', and 'petal-width'. All of the features are numeric. The records have been classified into one of the **three classes** i.e. 'Iris-setosa', 'Iris-versicolor', or 'Iris-verginica'."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataset Overview"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "iris = load_iris()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#source : https://towardsdatascience.com/pca-using-python-scikit-learn-e653f8989e60\n",
    "\n",
    "# The indices of the features that we are plotting\n",
    "x_index = 0\n",
    "y_index = 1\n",
    "\n",
    "# this formatter will label the colorbar with the correct target names\n",
    "formatter = plt.FuncFormatter(lambda i, *args: iris.target_names[int(i)])\n",
    "\n",
    "plt.figure(figsize=(5, 4))\n",
    "plt.scatter(iris.data[:, x_index], iris.data[:, y_index], c=iris.target)\n",
    "plt.colorbar(ticks=[0, 1, 2], format=formatter)\n",
    "plt.xlabel(iris.feature_names[x_index])\n",
    "plt.ylabel(iris.feature_names[y_index])\n",
    "\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(iris.data.shape) ## number of rows/columns\n",
    "print(iris.feature_names) ## feature names\n",
    "print(iris.target_names) ## categories "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic Cleaning of the Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = \"https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data\"  \n",
    "names = ['sepal-length', 'sepal-width', 'petal-length', 'petal-width', 'Class']  \n",
    "df = pd.read_csv(url, names=names)  \n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#divide the dataset into a feature set and corresponding labels.\n",
    "X = df.drop('Class', 1)  \n",
    "y = df['Class']  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Splitting the dataset into the Training set and Test set\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#PCA performs best with a normalized feature set. We will perform standard scalar normalization to normalize our feature set.\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "sc = StandardScaler()  \n",
    "X_train = sc.fit_transform(X_train)  \n",
    "X_test = sc.transform(X_test)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applying PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Performing PCA using Scikit-Learn is a two-step process:\n",
    "\n",
    "1. **Initialize the PCA class** by passing the number of components to the constructor.\n",
    "2. **Call the *fit* and then *transform* methods** by passing the feature set to these methods. The *transform* method returns the specified number of principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = skPCA()  \n",
    "X_train = pca.fit_transform(X_train)  \n",
    "X_test = pca.transform(X_test)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PCA class contains `explained_variance_ratio_` which returns the variance caused by each of the principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "explained_variance = pca.explained_variance_ratio_  \n",
    "print(explained_variance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be seen that first principal component is responsible for 72.22% variance. Similarly, the second principal component causes 23.9% variance in the dataset. Collectively we can say that (72.22 + 23.9) 96.21% percent of the classification information contained in the feature set is captured by the first two principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = skPCA(n_components=1)  \n",
    "X_train = pca.fit_transform(X_train)  \n",
    "X_test = pca.transform(X_test) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training and Making Predictions\n",
    "In this case we'll use random forest classification for making the predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "classifier = RandomForestClassifier(max_depth=2, random_state=0)  \n",
    "classifier.fit(X_train, y_train)\n",
    "\n",
    "# Predicting the Test set results\n",
    "y_pred = classifier.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Performance Evaluation\n",
    "from sklearn.metrics import confusion_matrix  \n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "cm = confusion_matrix(y_test, y_pred)  \n",
    "print(cm)  \n",
    "print( accuracy_score(y_test, y_pred))  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "it can be seen from the output that with only one feature, the random forest algorithm is able to correctly predict 28 out of 30 instances, resulting in 93.33% accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA on Larger Dataset and CPU vs GPU comparisson  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will first import some helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# measure the execution time \n",
    "from timeit import default_timer\n",
    "\n",
    "\n",
    "class Timer(object):  #Class for timing execution speed of small code snippets.\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data(nrows, ncols, ):\n",
    "    print('use random data')\n",
    "    X = np.random.rand(nrows,ncols)\n",
    "    df = pd.DataFrame({'fea%d'%i:X[:,i] for i in range(X.shape[1])})\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "\n",
    "def array_equal(a,b,threshold=2e-3,with_sign=True):\n",
    "    a = to_nparray(a)\n",
    "    b = to_nparray(b)\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": [
    "### Benchmark with Random Data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the larger dataset with more features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "## here, we will load the random data with 1M rows and 400 columns. Feel free to experiment with a different dataset sizes. \n",
    "nrows = 1000000\n",
    "ncols = 400\n",
    "\n",
    "X = load_data(nrows,ncols)\n",
    "print('data',X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look into all possible parameters that we can use when applying PCA :http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html\n",
    "\n",
    "We will start here with the following : \n",
    "\n",
    "- **n_components : int, float, None or string** <br>\n",
    "Number of components to keep. if `n_components` is not set all components are kept\n",
    "- **whiten : bool, optional (default False)** <br>\n",
    "When `True` (`False` by default) the `components_` vectors are multiplied by the square root of `n_samples` and then divided by the singular values to ensure uncorrelated outputs with unit component-wise variances.<br/>Whitening will remove some information from the transformed signal (the relative variance scales of the components) but can sometime improve the predictive accuracy of the downstream estimators by making their data respect some hard-wired assumptions\n",
    "\n",
    "- **random_state : int, RandomState instance or None, optional (default None)** <br>\n",
    "If int, `random_state` is the seed used by the random number generator\n",
    "\n",
    "- **svd_solver : string {‘auto’, ‘full’, ‘arpack’, ‘randomized’}** <br>\n",
    "If `\"full\"` :run exact full SVD calling the standard LAPACK solver via scipy.linalg.svd and select the components by postprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_components = 10\n",
    "whiten = False\n",
    "random_state = 42\n",
    "svd_solver=\"full\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the time needed to execute PCA  function using standard *sklearn* library. Note: this algorithm runs on CPU only. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiprocessing\n",
    "print(multiprocessing.cpu_count()) # Return the number of CPUs in the system. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from sklearn.decomposition import PCA as skPCA \n",
    "\n",
    "pca_sk = skPCA(n_components=n_components,svd_solver=svd_solver, \n",
    "            whiten=whiten, random_state=random_state)\n",
    "result_sk = pca_sk.fit_transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, before we execute PCA function using RAPIDS *cuml* library we will first read the data in  GPU data format using *cudf* :\n",
    "- cudf - GPU DataFrame manipulation library https://github.com/rapidsai/cudf\n",
    "- cuml - suite of libraries that implements a machine learning algorithms within the RAPIDS data science ecosystem https://github.com/rapidsai/cuml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "import cudf\n",
    "import cuml #documentaton can be found here https://cuml.readthedocs.io/en/latest/\n",
    "from cuml import PCA as cumlPCA  \n",
    "\n",
    "\n",
    "X = cudf.DataFrame.from_pandas(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will execute the PCA function using cuml and check the performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "pca_cuml = cumlPCA(n_components=n_components,svd_solver=svd_solver, \n",
    "            whiten=whiten, random_state=random_state)\n",
    "result_cuml = pca_cuml.fit_transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that for dataset of size 1000000 rows and 400 columns it takes around 8X less time to execute the PCA algorithm using RAPIDS cuml library.\n",
    "\n",
    "Block of the code below compares the attributes and results of two libraries : *sklearn* PCA and *cuml* PCA. Here, we see that results and attributes are exacty the same and that user will not see any difference in the existing workflow. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for attr in ['singular_values_','components_','explained_variance_',\n",
    "             'explained_variance_ratio_']:\n",
    "    passed = array_equal(getattr(pca_sk,attr),getattr(pca_cuml,attr))\n",
    "    message = 'compare pca: cuml vs sklearn {:>25} {}'.format(attr,'equal' if passed else 'NOT equal')\n",
    "    print(message)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "passed = array_equal(result_sk,result_cuml)\n",
    "message = 'compare pca: cuml vs sklearn transformed results %s'%('equal'if passed else 'NOT equal')\n",
    "print(message)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "To learn more about RAPIDS, be sure to check out: \n",
    "\n",
    "* [Open Source Website](http://rapids.ai)\n",
    "* [GitHub](https://github.com/rapidsai/)\n",
    "* [Press Release](https://nvidianews.nvidia.com/news/nvidia-introduces-rapids-open-source-gpu-acceleration-platform-for-large-scale-data-analytics-and-machine-learning)\n",
    "* [NVIDIA Blog](https://blogs.nvidia.com/blog/2018/10/10/rapids-data-science-open-source-community/)\n",
    "* [Developer Blog](https://devblogs.nvidia.com/gpu-accelerated-analytics-rapids/)\n",
    "* [NVIDIA Data Science Webpage](https://www.nvidia.com/en-us/deep-learning-ai/solutions/data-science/)\n"
   ]
  },
  {
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}