F0D-10-SparseMatrices.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Ms8lKNI8ZNph"
   },
   "source": [
    "# FLIP (00): Data Science \n",
    "**(Module 03: Linear Algebra)**\n",
    "\n",
    "---\n",
    "- Materials in this module include resources collected from various open-source online repositories.\n",
    "- You are free to use,but NOT allowed to change and distribute this package.\n",
    "\n",
    "Prepared by and for \n",
    "**Student Members** |\n",
    "2006-2023 [TULIP Lab](http://www.tulip.org.au), Australia\n",
    "\n",
    "---\n",
    "## Session 18 Sparse Matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "iMcQx4WAZNpq"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import sparse\n",
    "import scipy.sparse.linalg as spla\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "print('aaa')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bW1FPqdkZNps"
   },
   "outputs": [],
   "source": [
    "sns.set_context('notebook', font_scale=1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uhAOd8iIZNps"
   },
   "source": [
    "## Creating a sparse matrix\n",
    "\n",
    "There are many applications in which we deal with matrices that are mostly zeros. For example, a matrix representing social networks is very sparse - there are 7 billion people, but most people are only connected to a few hundred or thousand others directly. Storing such a social network as a sparse rather than dense matrix  will offer orders of magnitude reductions in memory requirements and corresponding speed-ups in computation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DxyrCy25ZNpt"
   },
   "source": [
    "### Coordinate format\n",
    "\n",
    "The simplest sparse  matrix format is built from the coordinates and values of the non-zero entries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6_2o08RwZNpt"
   },
   "source": [
    "#### From dense matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wa-XHPnZZNpu"
   },
   "outputs": [],
   "source": [
    "A = np.random.poisson(0.2, (5,15)) * np.random.randint(0, 10, (5, 15))\n",
    "#前者生成符合泊松分布的01矩阵，后者生成随机数，乘在前者的各个数值（0和1）上\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "F4XVIsM8ZNpu"
   },
   "outputs": [],
   "source": [
    "rows, cols = np.nonzero(A)\n",
    "vals = A[rows, cols]\n",
    "#输出非0值\n",
    "#np.nonzero(a)返回数组a中值不为零的元素的下标，(rows,cols)就是一个坐标，如A[1][5]=5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "CBlHUJIFZNpu"
   },
   "outputs": [],
   "source": [
    "vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "JiLDUGiqZNpu"
   },
   "outputs": [],
   "source": [
    "rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nNX71eT9ZNpv"
   },
   "outputs": [],
   "source": [
    "cols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ApgbFu-pZNpv"
   },
   "outputs": [],
   "source": [
    "X1 = sparse.coo_matrix(A)\n",
    "X1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "GXF8VcdWZNpv"
   },
   "outputs": [],
   "source": [
    "print(X1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VHr__1D3ZNpv"
   },
   "source": [
    "#### From coordinates\n",
    "\n",
    "Note that the (values, (rows, cols)) argument is a single tuple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "BubWtc6GZNpw"
   },
   "outputs": [],
   "source": [
    "X2 = sparse.coo_matrix((vals, (rows, cols)))\n",
    "X2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5zq-r5d6ZNpw"
   },
   "outputs": [],
   "source": [
    "print(X2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Jezwot1-ZNpw"
   },
   "source": [
    "#### Convert back to dense matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9yxx0gZBZNpw"
   },
   "outputs": [],
   "source": [
    "X2.todense()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BbyEcP1JZNpx"
   },
   "source": [
    "### Compressed Sparse Row and Column formats\n",
    "\n",
    "When we have 2 or more repeated entries in the rows or cols, we can remove the redundancy by indicating the  location of the  first occurrence of a value and its increment instead of the full coordinates. These are known as CSR or CSC formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "LagQbby_ZNpx"
   },
   "outputs": [],
   "source": [
    "np.vstack([rows, cols])\n",
    "#np.vstack沿着竖直方向将矩阵堆叠起来,具体参考https://blog.csdn.net/u012609509/article/details/70319293"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Z42syQejZNpx"
   },
   "outputs": [],
   "source": [
    "indptr = np.r_[np.searchsorted(rows, np.unique(rows)), len(rows)]\n",
    "indptr\n",
    "#np.r_是按列连接两个矩阵，就是把两矩阵上下相加，要求列数相等。\n",
    "#np.searchsorted在数组a中插入数组v（并不执行插入操作），返回一个下标列表，这个列表指明了v中对应元素应该插入在a中那个位置上，参考链接：\n",
    "#https://my.oschina.net/u/3906518/blog/1844187\n",
    "#np.unique对于一维数组或者列表，unique函数去除其中重复的元素，并按元素由大到小返回一个新的无元素重复的元组或者列表\n",
    "#indptr最后一个元素是len(rows),这是一个数字\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fyjcAfoFZNpx"
   },
   "outputs": [],
   "source": [
    "X3 = sparse.csr_matrix((vals, cols, indptr))\n",
    "X3\n",
    "#print(X3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wH5uqnVcZNpy"
   },
   "outputs": [],
   "source": [
    "X3.todense()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mVjx5xvzZNpy"
   },
   "source": [
    "#### Casting from COO format\n",
    "\n",
    "Because the coordinate format is more intuitive, it is often more convenient to first create a COO matrix then cast to CSR or CSC form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "u7PCwiVwZNpy"
   },
   "outputs": [],
   "source": [
    "X4 = X2.tocsr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "04sfLD0lZNpy"
   },
   "outputs": [],
   "source": [
    "X4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BN3sEqcFZNpy"
   },
   "source": [
    "### COO summation convention\n",
    "\n",
    "When entries are repeated in a COO matrix, they are **summed**. This provides a quick way to construct confusion matrices for evaluation of multi-class  classification algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "oMcMZC7QZNpz"
   },
   "outputs": [],
   "source": [
    "rows = np.r_[np.zeros((4),dtype=np.int),np.ones((4),dtype=np.int)]\n",
    "cols = np.repeat([0,1], 4)\n",
    "vals = np.arange(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cxADiZd4ZNpz"
   },
   "outputs": [],
   "source": [
    "rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "e8t-DsZfZNpz"
   },
   "outputs": [],
   "source": [
    "cols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "MzGCy-K4ZNpz"
   },
   "outputs": [],
   "source": [
    "vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "gNL22hVnjon7"
   },
   "outputs": [],
   "source": [
    "# b = np.zeros((4),dtype=np.int)\n",
    "# print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "q7oeNUlZZNp0"
   },
   "outputs": [],
   "source": [
    "X5 = sparse.csr_matrix((vals, (rows, cols)),shape=(5,15)).toarray()\n",
    "#原命令:X5 = sparse.csr_matrix((vals, (rows, cols))\n",
    "#错误在于没有指定shape，也没有把结果转化成array形式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wAZqdvkgZNp0"
   },
   "outputs": [],
   "source": [
    "print(X5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BASxaSsqZNp0"
   },
   "source": [
    "## Application: Confusion matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-uJcPM6kZNp0"
   },
   "source": [
    "### Creating a 2 by 2 confusion matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "M_6Cq_4WZNp0"
   },
   "outputs": [],
   "source": [
    "obs = np.random.randint(0, 2, 100)\n",
    "pred = np.random.randint(0, 2, 100)\n",
    "vals = np.ones(100).astype('int')\n",
    "#print(obs)\n",
    "#print(len(obs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "S2XTVVXlZNp1"
   },
   "outputs": [],
   "source": [
    "pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6vQ1oFDuZNp1"
   },
   "outputs": [],
   "source": [
    "vals.shape, obs.shape , pred.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "uWLeawX6ZNp1"
   },
   "outputs": [],
   "source": [
    "X6 = sparse.coo_matrix((vals, (pred, obs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "DP0YQI8FZNp1"
   },
   "outputs": [],
   "source": [
    "X6.todense()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-s_-5S0lZNp1"
   },
   "source": [
    "#### Creating an $n$ by $n$ confusion matrix\n",
    "\n",
    "For classifications with a large number of classes (e.g. image segmentation),  the savings are even more dramatic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "HI7EUI4uZNp2"
   },
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.neighbors import KNeighborsClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "zhmmKIx4ZNp2"
   },
   "outputs": [],
   "source": [
    "iris = datasets.load_iris()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "56R8TwXDZNp2"
   },
   "outputs": [],
   "source": [
    "knn = KNeighborsClassifier()\n",
    "X_train, X_test, y_train, y_test = train_test_split(iris.data, iris.target, \n",
    "                                                    test_size=0.5, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "SpayU0UyZNp2"
   },
   "outputs": [],
   "source": [
    "pred = knn.fit(X_train, y_train).predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7keS7BkZZNp2"
   },
   "outputs": [],
   "source": [
    "X7 = sparse.coo_matrix((np.ones(len(pred)).astype('int'), (pred, y_test)))\n",
    "pd.DataFrame(X7.todense(), index=iris.target_names, columns=iris.target_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CeaAWy0gZNp3"
   },
   "source": [
    "## Application: PageRank\n",
    "\n",
    "SciPy provides efficient routines for solving large sparse systems as for dense matrices. We will illustrate by calculating the page rank  for airports using data from the [Bureau of Transportation Statisitcs](http://www.transtats.bts.gov/DL_SelectFields.asp?Table_ID=236)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "b0OVKtFXZNp3"
   },
   "outputs": [],
   "source": [
    "data = pd.read_csv('/content/drive/MyDrive/Data/airports.csv', usecols=[0,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1BZ7Cs3pZNp3"
   },
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "P4b_DFr2q_5N"
   },
   "outputs": [],
   "source": [
    "from google.colab import drive\n",
    "drive.mount('/content/drive')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "BPdngfjyZNp3"
   },
   "outputs": [],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ANiYnIjOZNp3"
   },
   "outputs": [],
   "source": [
    "lookup = pd.read_csv('/content/drive/MyDrive/Data/names.csv', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bBxfhU9dZNp4"
   },
   "outputs": [],
   "source": [
    "lookup.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "q575v2pXZNp4"
   },
   "outputs": [],
   "source": [
    "lookup.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "EkAdHAULZNp4"
   },
   "outputs": [],
   "source": [
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jfFLWRpAZNp4"
   },
   "source": [
    "### Construct the sparse adjacency matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "egDweZrIZNp4"
   },
   "outputs": [],
   "source": [
    "g = nx.from_pandas_edgelist(data, source='ORIGIN_AIRPORT_ID', target='DEST_AIRPORT_ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Rg0LBKYZZNp5"
   },
   "outputs": [],
   "source": [
    "airports = np.array(g.nodes())\n",
    "adj_matrix = nx.to_scipy_sparse_matrix(g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "43Oq0Q3tZNp5"
   },
   "source": [
    "### Construct the transition matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Tipy3o3yZNp6"
   },
   "outputs": [],
   "source": [
    "out_degrees = np.ravel(adj_matrix.sum(axis=1))\n",
    "diag_matrix = sparse.diags(1 / out_degrees).tocsr()\n",
    "M = (diag_matrix @ adj_matrix).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QghncrZSZNp6"
   },
   "source": [
    "### Modify the transition matrix with a damping factor\n",
    "\n",
    "The PageRank algorithm assumes that every node can be reached from every other node. To guard against case where a node has out-degree 0, we allow every node a small random chance of transitioning to any other node using a damping factor $d$. Then we solve  the linear system  to find the pagerank score $r$.\n",
    "\n",
    "$$\n",
    "r = (I - dM)^{-1}\\frac{1-d}{N}\\mathbb{1}\n",
    "$$\n",
    "\n",
    "or equivalently in the $Ax = b$ format\n",
    "\n",
    "$$\n",
    "(I - dM)r = \\frac{1-d}{N}\\mathbb{1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "JvhUF2heZNp6"
   },
   "outputs": [],
   "source": [
    "n = len(airports) \n",
    "d = 0.85 \n",
    "I = sparse.eye(n, format='csc')\n",
    "A = I - d * M\n",
    "b = (1-d) / n * np.ones(n) # so the sum of all page ranks is 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "g0bt3AjzZNp6"
   },
   "outputs": [],
   "source": [
    "A.todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ebPbeT3HZNp7"
   },
   "outputs": [],
   "source": [
    "from scipy.sparse.linalg import spsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5kK16Vz8ZNp7"
   },
   "outputs": [],
   "source": [
    "r =  spsolve(A, b)\n",
    "r.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8vGJpSgfZNp7"
   },
   "outputs": [],
   "source": [
    "idx = np.argsort(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cI186ModZNp7"
   },
   "outputs": [],
   "source": [
    "top10 = idx[-10:][::-1]\n",
    "bot10 = idx[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "AJ-Iz24RZNp7"
   },
   "outputs": [],
   "source": [
    "df = lookup.loc[airports[top10]]\n",
    "df['degree'] = out_degrees[top10]\n",
    "df['pagerank']= r[top10]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "OvFjrUBfZNp7"
   },
   "outputs": [],
   "source": [
    "df = lookup.loc[airports[bot10]]\n",
    "df['degree'] = out_degrees[bot10]\n",
    "df['pagerank']= r[bot10]\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bd5PEK5rZNp8"
   },
   "source": [
    "### Visualize the airport connections graph and label the top and bottom 5 airports by pagerank"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "SvXtzDKMZNp8"
   },
   "outputs": [],
   "source": [
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "C0xYtEsbZNp8"
   },
   "outputs": [],
   "source": [
    "labels = {airports[i]: lookup.loc[airports[i]].str.split(':').str[0].values[0] \n",
    "          for i in np.r_[top10[:5], bot10[:5]]}\n",
    "\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter('ignore')\n",
    "    nx.draw(g, pos=nx.spring_layout(g), labels=labels, \n",
    "            node_color='blue', font_color='red', alpha=0.5,\n",
    "            node_size=np.clip(5000*r, 1, 5000*r), width=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "HnKd_6ROZNp8"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "colab": {
   "collapsed_sections": [],
   "name": "F0D-10-SparseMatrices.ipynb",
   "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.6"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}