Advanced Numerical Linear Algebra/week_3_topic_modeling_with_NMF_and_SVD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Topic Modeling with NMF and SVD\n",
    "        Topic Modeling is a great way to get started with matrix factorizations. We'll start with term-document matrix:\n",
    "\n",
    "        We can decompose this into one tall thin matrix times one wide short matrix.\n",
    "\n",
    "        This representation does not take into account word order or sentence structure. It's an example of a bag of words approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Motiation \n",
    "\n",
    "        Consider the most extreme case - reconstructing the matrix using an outer product of two vectors. Clearly, in most cases we won't be able to reconstruct the matrix exactly. But if we had one vector with the relative frequency of each vocabulary word out of the total word count, and one with the average number of words per document, then that outer product would be as close as we can get.\n",
    "\n",
    "        Now consider increasing that matrices to two columns and two rows. The optimal decomposition would now be to cluster the documents into two groups, each of which has as different a distribution of words as possible to each other, but as similar as possible amongst the documents in the cluster. We will call those two groups \"topics\". And we would cluster the words into two groups, based on those which most frequently appear in each of the topics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In this week's lecture,\n",
    "        We'll take a dataset of documents in several different categories, and find topics for them. Knowing the actual categories help us evaluate if the topics we find make sense.\n",
    "\n",
    "\n",
    "        We'll try this with : Singular Value Decomposition and Non-Negative Matrix Factorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.datasets import fetch_20newsgroups\n",
    "from sklearn import decomposition\n",
    "from scipy import linalg\n",
    "import matplotlib.pyplot as plt  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Up Data\n",
    "        Scikit Learn comes with a number of built in datasets, as well as loading utilities to load several standard external datasets.  \n",
    "\n",
    "        Newsgroups are discussion groups on USenet, which were popular in the 80s and 90s before the web really took off. This dataset includes 18,000 newsgroups posts with 20 topics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2034,), (2034,))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "categories = ['alt.atheism', 'talk.religion.misc', 'comp.graphics', 'sci.space']\n",
    "remove = ('headers', 'footers', 'quotes')\n",
    "newsgroups_train = fetch_20newsgroups(subset='train', categories=categories, remove=remove)\n",
    "newsgroups_test = fetch_20newsgroups(subset='test', categories=categories, remove=remove)\n",
    "newsgroups_train.filenames.shape, newsgroups_train.target.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        Look at some of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hi,\n",
      "\n",
      "I've noticed that if you only save a model (with all your mapping planes\n",
      "positioned carefully) to a .3DS file that when you reload it after restarting\n",
      "3DS, they are given a default position and orientation.  But if you save\n",
      "to a .PRJ file their positions/orientation are preserved.  Does anyone\n",
      "know why this information is not stored in the .3DS file?  Nothing is\n",
      "explicitly said in the manual about saving texture rules in the .PRJ file. \n",
      "I'd like to be able to read the texture rule information, does anyone have \n",
      "the format for the .PRJ file?\n",
      "\n",
      "Is the .CEL file format available from somewhere?\n",
      "\n",
      "Rych\n",
      "\n",
      "\n",
      "Seems to be, barring evidence to the contrary, that Koresh was simply\n",
      "another deranged fanatic who thought it neccessary to take a whole bunch of\n",
      "folks with him, children and all, to satisfy his delusional mania. Jim\n",
      "Jones, circa 1993.\n",
      "\n",
      "\n",
      "Nope - fruitcakes like Koresh have been demonstrating such evil corruption\n",
      "for centuries.\n",
      "\n",
      " >In article <1993Apr19.020359.26996@sq.sq.com>, msb@sq.sq.com (Mark Brader) \n",
      "\n",
      "MB>                                                             So the\n",
      "MB> 1970 figure seems unlikely to actually be anything but a perijove.\n",
      "\n",
      "JG>Sorry, _perijoves_...I'm not used to talking this language.\n",
      "\n",
      "Couldn't we just say periapsis or apoapsis?\n",
      "\n",
      " \n",
      "I have a request for those who would like to see Charley Wingate\n",
      "respond to the \"Charley Challenges\" (and judging from my e-mail, there\n",
      "appear to be quite a few of you.)  \n",
      "\n",
      "It is clear that Mr. Wingate intends to continue to post tangential or\n",
      "unrelated articles while ingoring the Challenges themselves.  Between\n",
      "the last two re-postings of the Challenges, I noted perhaps a dozen or\n",
      "more posts by Mr. Wingate, none of which answered a single Challenge.  \n",
      "\n",
      "It seems unmistakable to me that Mr. Wingate hopes that the questions\n",
      "will just go away, and he is doing his level best to change the\n",
      "subject.  Given that this seems a rather common net.theist tactic, I\n",
      "would like to suggest that we impress upon him our desire for answers,\n",
      "in the following manner:\n",
      "\n",
      "1. Ignore any future articles by Mr. Wingate that do not address the\n",
      "Challenges, until he answers them or explictly announces that he\n",
      "refuses to do so.\n",
      "\n",
      "--or--\n",
      "\n",
      "2. If you must respond to one of his articles, include within it\n",
      "something similar to the following:\n",
      "\n",
      "    \"Please answer the questions posed to you in the Charley Challenges.\"\n",
      "\n",
      "Really, I'm not looking to humiliate anyone here, I just want some\n",
      "honest answers.  You wouldn't think that honesty would be too much to\n",
      "ask from a devout Christian, would you?  \n",
      "\n",
      "Nevermind, that was a rhetorical question.\n"
     ]
    }
   ],
   "source": [
    "print(\"\\n\".join(newsgroups_train.data[:4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        In what categories the data belong to?\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['comp.graphics', 'talk.religion.misc', 'sci.space', 'alt.atheism'],\n",
       "      dtype='<U18')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(newsgroups_train.target_names)[newsgroups_train.target[:4]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        We assign target index to each of them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 3, 2, 0, 2, 0, 2, 1, 2, 1], dtype=int64)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.target[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_topics, num_top_words = 6, 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        Scikit Learn has a method that will extract all the words count for us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2034, 26576)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectorizer = CountVectorizer(stop_words='english')\n",
    "vectors = vectorizer.fit_transform(newsgroups_train.data).todense() # (documents, vocab)\n",
    "vectors.shape #, vectors.nnz / vectors.shape[0], row_means.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2034 (2034, 26576)\n"
     ]
    }
   ],
   "source": [
    "print(len(newsgroups_train.data), vectors.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "vocab = np.array(vectorizer.get_feature_names_out())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['cosmonauts', 'cosmos', 'cosponsored', 'cost', 'costa', 'costar',\n",
       "       'costing', 'costly', 'costruction', 'costs', 'cosy', 'cote',\n",
       "       'couched', 'couldn', 'council', 'councils', 'counsel',\n",
       "       'counselees', 'counselor', 'count', 'countdown', 'counted',\n",
       "       'counter', 'counter_clockwise', 'counterargument',\n",
       "       'counterclockwise', 'countered', 'counterexamples',\n",
       "       'counterfactual', 'counterpart', 'counterproductive', 'counters',\n",
       "       'counting', 'countless', 'countries', 'country', 'countryside',\n",
       "       'counts', 'county', 'coup', 'couple', 'coupled', 'couples',\n",
       "       'courage', 'courageous', 'courant', 'cournoyer', 'course',\n",
       "       'courses', 'court'], dtype=object)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab[7000:7050]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Singular Value Decomposition\n",
    "\n",
    "    We would expect that the words that appear most frequently in one topic would appear less frequently in the other. Therefore the topics should be orthogonal\n",
    "\n",
    "    The SVD algorithm factorizes a matrix into one matrix with orthogonal columns and one with orthogonal rows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        SVD is an exact decompostion, since the matrices it creates are big enough to fully cover the original matrix. SVD is extremely widely used in linear algebra, and specifically in data sciences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 1min 33s\n"
     ]
    }
   ],
   "source": [
    "%time U, s, Vh = linalg.svd(vectors, full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2034, 2034) (2034,) (2034, 26576)\n"
     ]
    }
   ],
   "source": [
    "print(U.shape, s.shape, Vh.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZH0lEQVR4nO3de5Ac51nv8e/TPXuTdnVfy4qkRHaiXFwUsR0lCEhSEB+M40BkLkmZwzkWOT5xUSSQVKCCIXU4UEUVMZyTHFxQBoNdkVOGJCSkrEoZEuM4mFMVX1a+23KstWNFkmVpraulvc7Mwx/9zmzPbq/2orn16Pepmpqet3tmHvWsfv3O2z3d5u6IiEhniVpdgIiI1J/CXUSkAyncRUQ6kMJdRKQDKdxFRDpQodUFAKxbt863bNnS6jJERHJlz549r7n7YNa8tgj3LVu2MDQ01OoyRERyxcz2zzVPwzIiIh1I4S4i0oEU7iIiHUjhLiLSgRTuIiIdSOEuItKBFO4iIh0o1+H+6MvH+cJ3fsBksdzqUkRE2kquw33P/hPc+t1himWFu4hIWq7D3cK9rjciIlIr3+Ee0l3ZLiJSK9/hXu27i4hIWq7DvULXgRURqZXrcNewjIhItlyHe4U67iIitXId7qauu4hIpnyHe7h3pbuISI18h3ul465sFxGpke9wD/fKdhGRWvkOd9Nx7iIiWXId7hU6zl1EpNaCw93MYjN73My+FR5fYmYPm9mwmX3VzLpDe094PBzmb2lQ7TrOXURkDovpuX8K2Jt6fAvwRXd/C3ACuDG03wicCO1fDMs1hE4cJiKSbUHhbmabgA8Bfx8eG/AB4OthkV3AdWF6R3hMmH+VNWpwPLysDoUUEam10J77/wM+C1ROnL4WOOnuxfD4ILAxTG8EDgCE+afC8jXM7CYzGzKzoZGRkSUVX91iKNtFRGrMG+5m9gvAUXffU883dvfb3X2bu28bHBxc0mtozF1EJFthAcv8NPBhM7sW6AVWAH8JrDKzQuidbwIOheUPAZuBg2ZWAFYCx+peOdOn/NWYu4hIrXl77u7+B+6+yd23ANcD33X3XwceAH41LLYTuCdM7w6PCfO/6w06VnG65650FxFJO5/j3H8f+IyZDZOMqd8R2u8A1ob2zwA3n1+Jc9NPmEREsi1kWKbK3b8HfC9MvwS8J2OZceAjdahtEXU1891ERNpfrn+hqh2qIiLZ8h3u1R2qincRkbRchzs65a+ISKZch7t2qIqIZMt3uJuOcxcRyZLvcA/3Os5dRKRWvsNd4zIiIplyHe4VGpYREamV63DXce4iItnyHe46zl1EJFO+w109dxGRTLkO9wp13EVEauU63Kev3qd0FxFJy3e4h3v13EVEauU73HWcu4hIplyHe4U67iIitXId7rqGqohItnyHu66hKiKSKd/hHu7VcxcRqZXvcNfFOkREMuU63Ct9dw3LiIjUynW4q+cuIpIt3+He6gJERNpUrsNdRESy5TrcdQ1VEZFs+Q73cK8dqiIitfId7tqhKiKSqTPCvbVliIi0nXyHuy6zJyKSKdfhjnruIiKZch3uOreMiEi2fIe7rtYhIpIp1+E+TV13EZG0XIe7hmVERLLlO9y1Q1VEJFO+w12X2RMRyTRvuJtZr5k9YmZPmtmzZvYnof0SM3vYzIbN7Ktm1h3ae8Lj4TB/S6OKn/6FqtJdRCRtIT33CeAD7v5O4HLgGjPbDtwCfNHd3wKcAG4My98InAjtXwzLNcT0uWVERCRt3nD3xJnwsCvcHPgA8PXQvgu4LkzvCI8J86+yRh2zqHPLiIhkWtCYu5nFZvYEcBS4D3gROOnuxbDIQWBjmN4IHAAI808BazNe8yYzGzKzoZGRkSUVb7pch4hIpgWFu7uX3P1yYBPwHuDt5/vG7n67u29z922Dg4Pn91oamBERqbGoo2Xc/STwAPCTwCozK4RZm4BDYfoQsBkgzF8JHKtHsTOZBt1FRDIt5GiZQTNbFab7gJ8D9pKE/K+GxXYC94Tp3eExYf53vUGHsyjbRUSyFeZfhA3ALjOLSTYGX3P3b5nZc8BXzOxPgceBO8LydwBfNrNh4DhwfQPqBnSZPRGRucwb7u7+FHBFRvtLJOPvM9vHgY/Upbp5TP9CVekuIpKW81+oJtRzFxGple9w17llREQy5Trc0XHuIiKZch7uCZ1bRkSkVq7DXcMyIiLZ8h3ulQmlu4hIjXyHe+U4d6W7iEiNfId7uNeQu4hIrXyHu075KyKSKd/hXrnMXovrEBFpN/kOdx3mLiKSKdfhXqHj3EVEanVGuLe6ABGRNpPrcNcOVRGRbPkOd12uQ0QkU77DXT13EZFMnRHurS1DRKTt5Dvc0WX2RESy5DvcdZk9EZFM+Q73VhcgItKmch3uFRqWERGpletwr5zyt6x0FxGpketwL0RJuJfKCncRkbRch3scwr2ocBcRqZHrcC/E6rmLiGTJd7hHSfnquYuI1Mp5uIeee6nc4kpERNpLrsM9jjXmLiKSJdfhrqNlRESy5TrcdbSMiEi2fId7+BFTsaRwFxFJy3e4V4Zl9AtVEZEauQ53MyMyKGtYRkSkRq7DHZLeu3ruIiK1ch/ukZl67iIiM+Q+3AuR6WgZEZEZch/uUWQ6zl1EZIZ5w93MNpvZA2b2nJk9a2afCu1rzOw+M9sX7leHdjOzW81s2MyeMrMrG/kPiCPT+dxFRGZYSM+9CPyuu18GbAc+YWaXATcD97v7VuD+8Bjgg8DWcLsJuK3uVacU1HMXEZll3nB398Pu/liYfh3YC2wEdgC7wmK7gOvC9A7gLk88BKwysw31LrwiMoW7iMhMixpzN7MtwBXAw8B6dz8cZr0KrA/TG4EDqacdDG0zX+smMxsys6GRkZHF1l0Vq+cuIjLLgsPdzPqBbwCfdvfT6Xnu7sCiEtbdb3f3be6+bXBwcDFPrRGZjnMXEZlpQeFuZl0kwX63u/9zaD5SGW4J90dD+yFgc+rpm0JbQxRiHecuIjLTQo6WMeAOYK+7fyE1azewM0zvBO5Jtd8QjprZDpxKDd/UXWw6zl1EZKbCApb5aeC/A0+b2ROh7Q+BzwNfM7Mbgf3AR8O8e4FrgWFgFPhYPQueKdKhkCIis8wb7u7+/wGbY/ZVGcs78InzrGvBYh0tIyIyS+5/oZocLdPqKkRE2kuHhLvSXUQkLffhHkWGLsQkIlIr9+Ee62IdIiKz5D7cC1GkHaoiIjPkPtyjCIW7iMgMuQ93XWZPRGS23Ie7zgopIjJb7sNd53MXEZkt9+EeR5HOLSMiMkPuw70r1o+YRERmyn24x5HOCikiMlPuw11j7iIis+U+3OMooqjzD4iI1Mh9uHfFRlFj7iIiNXIf7rpAtojIbLkP94J2qIqIzJL7cNeYu4jIbLkP94LG3EVEZsl/uGvMXURklo4Id425i4jUyn24x1GEu87pLiKSlvtwL8QGoHF3EZGU/Id7lIS7eu4iItNyH+5xVOm5K9xFRCpyH+6VnruOdRcRmZb7cF/WUwDg9NhUiysREWkfuQ/3jav6ADh8arzFlYiItI/ch/tAb9JzPztRbHElIiLtI/fhvjwMy5ydVLiLiFTkPtz7Q7ifUc9dRKQq9+Fe7bkr3EVEqnIf7su6YgDOTJRaXImISPvIfbhHkbF+RQ8/fO1sq0sREWkbuQ93gE2rl3Hi7GSryxARaRsdEe7LumNGdbSMiEjVvOFuZnea2VEzeybVtsbM7jOzfeF+dWg3M7vVzIbN7Ckzu7KRxVck4a4xdxGRioX03L8EXDOj7WbgfnffCtwfHgN8ENgabjcBt9WnzHNb3l3Qce4iIinzhru7Pwgcn9G8A9gVpncB16Xa7/LEQ8AqM9tQp1rn1NcdM6aeu4hI1VLH3Ne7++Ew/SqwPkxvBA6kljsY2mYxs5vMbMjMhkZGRpZYRmJ5T4GzOhRSRKTqvHeoursDiz7frrvf7u7b3H3b4ODgedWwrDtmbKqkC3aIiARLDfcjleGWcH80tB8CNqeW2xTaGmpZd/JDprEp9d5FRGDp4b4b2BmmdwL3pNpvCEfNbAdOpYZvGmZZd3IKAh0OKSKSKMy3gJn9I/AzwDozOwj8b+DzwNfM7EZgP/DRsPi9wLXAMDAKfKwBNc+yvCfpuY9OlGCgGe8oItLe5g13d/+1OWZdlbGsA58436IWq69Lp/0VEUnriF+oVnruOhxSRCTREeFeGXN/Xaf9FREBOiTc37CqF4BXTo61uBIRkfbQEeG+fqCX7kLEj46NtroUEZG20BHhHkXG5tV97Fe4i4gAHRLuABcN9HLs7ESryxARaQsdE+467a+IyLTOCfeegsJdRCTomHBf19/N4VNjjOv8MiIinRPuP/u2ixifKvP9F4+1uhQRkZbrmHD/8U0rAXhx5EyLKxERab2OCfeVfV0M9BT40XEdDiki0jHhbmZsXrNM4S4iQgeFO8A7NqzgkR8ep6wrMonIBa6jwn37pWsYnSzx6Mszr+ctInJh6ahw/8DbLwLgiQMnW1uIiEiLdVS4r+3v4bINK/j6noMk1w0REbkwdVS4A3z8/Zew7+gZdj/5SqtLERFpmY4L9w+/cyPv2LCCW/7leV2ZSUQuWB0X7nFk/K8PvYNXTo3zR/c80+pyRERaouPCHeCn3rKOj7/vEv5pz0EeeP5oq8sREWm6jgx3gN+9+m1svaif37r7Mb797KutLkdEpKk6Ntx7u2Lu/p8/wVsvHuC3/+Fx7nvuSKtLEhFpmo4Nd4CLVvTypd94N5cOLufjdw1x2/debHVJIiJN0dHhDrB6eTe7P/lerr5sPbf86/P82b17GZ0strosEZGG6vhwB+guRNz2397F9e/ezN8++BI/+3++x9cePUBJ56ARkQ51QYQ7JIdIfv5Xfpyv/+ZPsmFlH5/9xlN86Nb/4MEXRlpdmohI3V0w4V6xbcsavvlbP8Vf/dcrODtZ5IY7H+GGOx/h+y8eo1gqt7o8EZG6sHY4B8u2bdt8aGio6e87USzx5e/v59b793F6vMja5d3suHwjH9m2ibdfPICZNb0mEZGFMrM97r4tc96FHO4VZyaKPPjCCN966hXue+4IUyVncKCHqy9bz7vetJrtl65lw8pehb2ItBWF+yIcOzPBd547wn/sG+GB50cYm0rOT7Ouv5t3bFgRbgO8/eIVvHmwn+7CBTeyJSJtQuG+RKWys/fwafbsP8HTh07x/KuneeHIGSaLydh8V2y8ebCfresH2LS6L9yWsWl1HxtX9dHbFbf4XyAinexc4V5odjF5EkfGj21cyY9tXFltK5bKvPTaWfYePs3ew6+z9/Bpnjxwkn995jBTpdoN5drl3Vy8speLV/SyfmUv6wd6uWhFD+tX9HDRQC8XDfSwenk3XbF6/yJSXwr3RSrEEW9dP8Bb1w+w4/Lp9lLZOXJ6nEMnxzh4YpQDx8c4fGqcV0+N8cqpcR4/cJLjZyczX3NFb4G1/T2sXd7NmuXdrO1P7lcv62ZlXxerwn3ltqKvQF9XrH0AIjInhXudxJHxhlV9vGFVH+/esiZzmclimZEzExw9Pc6R0xOMvD7OsbOTHD87mdyfmWT/sVEe+9EJToxOnfNHVl2x0d9ToL+3QH9PFwNhenlPgf6eAgO9yQZgWXdy6+suhPuYZV0xy7oL9HVHSXtX0t5TiLTBEOkQCvcm6i5EbFyVjMfPp1x2Xp8ocmp0ilNjye3k2CSnx4qcHk8en50o8vp4cjszMcXI6xP88LWznJkocma8WN0ZvFBxZPQWIvq6Y3q7krDviiN6ChHdYbq7ENEdR3QVInriGe1hXvq+dp6F+zjMSx73ZLx2d5zcokgbG5GlULi3qSiy6jDMUpXLzthUidHJEmOTJUanitPTkyVGJ4vV6WS5IuNTZcamSoxPlpgolpkslZksTt/OTBST6VT7VGW6VJ613+F8dcVWG/xzbGy6Z22IpjckXQWjJzy3EEcUIqMQG4XIiKOoOl2ZF0dGV5zM6wqPk2Ui4srjyIgiqy4/3R4RmxGH14xselmRZmpIuJvZNcBfAjHw9+7++Ua8j5xbFBnLe5KhmmYpl52pcjr4PQR/icmiz9ooVDYgU6kNRrU9Y96sjUp4PDpaZLLkTBZL4Tk+6zmtZEYS+rM2DBFxRM2Go7qMJRuVOLLqc9Pz0xuUmg1NxsYlDu8Z2/R9HJGaTpZNLxdHVNvi1GvXPKfaxqy29GumXyv9mtXpymvMatNGcanq/r/ezGLgr4GfAw4Cj5rZbnd/rt7vJe0nioyeKKan0F6Hgbo7xbJTKjtTpXK4Tx4Xy2WKpdnzK+2lsjNVdoqhvZx6rWS56enK43K1vUypDKVyOXnsTim8V/V1SqG95rnlzNefLJannxtqm35umXKZWc+t1FIOy+XtfHlxxgZk5sYhythwVjaO6W9P8Yy29LeyQmx0Ve7j5NtfIU6+vRXi2vmV9q44SvZr9RTo70n2ZfWHDtXqZV0t3YfViC7de4Bhd38JwMy+AuwAFO7SMmYWhni44H9/4J4EfCkV+CVPNgLT01Q3RKXqRmF6I1OuaZv7tarzZ73+zPek2lYqz5gfXiP9/tk1J69TSm3cqhu4sDEslsuMF6dft7JRL4YhxWI5uZ8qlavLL3WoccPKXvoX8K35d67ayi++8w1Leo9zaUS4bwQOpB4fBH5i5kJmdhNwE8Ab3/jGBpQhIlnMjNiSHrHMz1PfqiqhPxW+1U2VytX9V2cmSoxOFDk7WeLE2UmePHiS8gJ+JHo++9XOpWU7VN39duB2SH6h2qo6RETOxcJQTiFn3/oa8dPIQ8Dm1ONNoU1ERJqkEeH+KLDVzC4xs27gemB3A95HRETmUPdhGXcvmtkngW+THAp5p7s/W+/3ERGRuTVkzN3d7wXubcRri4jI/HQ6QhGRDqRwFxHpQAp3EZEOpHAXEelAbXGZPTMbAfYv8enrgNfqWE69tGNd7VgTtGddqmnh2rGudqwJ6l/Xm9x9MGtGW4T7+TCzobmuIdhK7VhXO9YE7VmXalq4dqyrHWuC5talYRkRkQ6kcBcR6UCdEO63t7qAObRjXe1YE7RnXapp4dqxrnasCZpYV+7H3EVEZLZO6LmLiMgMCncRkQ6U63A3s2vM7AdmNmxmNzfxfTeb2QNm9pyZPWtmnwrtf2xmh8zsiXC7NvWcPwh1/sDMfr6Btb1sZk+H9x8KbWvM7D4z2xfuV4d2M7NbQ11PmdmVDajnban18YSZnTazT7diXZnZnWZ21MyeSbUtet2Y2c6w/D4z29mAmv7CzJ4P7/tNM1sV2reY2Vhqnf1N6jnvCp/7cKh7yZdZmqOmRX9e9f7/OUddX03V9LKZPRHam7Wu5sqClv5dAcklpPJ4Izmd8IvApUA38CRwWZPeewNwZZgeAF4ALgP+GPi9jOUvC/X1AJeEuuMG1fYysG5G258DN4fpm4FbwvS1wL8ABmwHHm7CZ/Yq8KZWrCvg/cCVwDNLXTfAGuClcL86TK+uc01XA4UwfUuqpi3p5Wa8ziOhTgt1f7DONS3q82rE/8+sumbM/7/AHzV5Xc2VBS39u3L3XPfcqxfidvdJoHIh7oZz98Pu/liYfh3YS3Lt2LnsAL7i7hPu/kNgmKT+ZtkB7ArTu4DrUu13eeIhYJWZbWhgHVcBL7r7uX6N3LB15e4PAscz3m8x6+bngfvc/bi7nwDuA66pZ03u/h13L4aHD5FczWxOoa4V7v6QJ0lxV+rfUZeazmGuz6vu/z/PVVfofX8U+MdzvUYD1tVcWdDSvyvI97BM1oW4zxWwDWFmW4ArgIdD0yfD1607K1/FaG6tDnzHzPZYchFygPXufjhMvwqsb0FdkFyVK/2fr9XrCha/bppd3/8g6elVXGJmj5vZv5vZ+1K1HmxCTYv5vJq9nt4HHHH3fam2pq6rGVnQ8r+rPId7y5lZP/AN4NPufhq4DXgzcDlwmORrYrO9192vBD4IfMLM3p+eGXorTT/+1ZJLLn4Y+KfQ1A7rqkar1s1czOxzQBG4OzQdBt7o7lcAnwH+wcxWNKmctvu8Zvg1ajsOTV1XGVlQ1aq/qzyHe0svxG1mXSQf5t3u/s8A7n7E3UvuXgb+junhhKbV6u6Hwv1R4JuhhiOV4ZZwf7TZdZFsbB5z9yOhvpavq2Cx66Yp9ZnZbwC/APx6CAfC0MexML2HZEz7reH900M3da9pCZ9X0z5HMysAvwx8NVVv09ZVVhbQBn9XeQ73ll2IO4zv3QHsdfcvpNrT49W/BFT26u8GrjezHjO7BNhKslOn3nUtN7OByjTJjrlnwvtX9r7vBO5J1XVD2IO/HTiV+ipZbzU9q1avq5TFrptvA1eb2eowNHF1aKsbM7sG+CzwYXcfTbUPmlkcpi8lWTcvhbpOm9n28Ld5Q+rfUa+aFvt5NfP/538Bnnf36nBLs9bVXFlAO/xdnc/e2FbfSPY8v0CyVf5cE9/3vSRfs54Cngi3a4EvA0+H9t3AhtRzPhfq/AHnsXd+nrouJTkq4Ung2co6AdYC9wP7gH8D1oR2A/461PU0sK1BdS0HjgErU21NX1ckG5fDwBTJmOaNS1k3JOPgw+H2sQbUNEwy/lr52/qbsOyvhM/1CeAx4BdTr7ONJHBfBP6K8OvzOta06M+r3v8/s+oK7V8CfnPGss1aV3NlQUv/rtxdpx8QEelEeR6WERGROSjcRUQ6kMJdRKQDKdxFRDqQwl1EpAMp3EVEOpDCXUSkA/0nv01Kyro8oL4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot for the singular values\n",
    "plt.plot(s);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2644c7470b8>]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhR0lEQVR4nO3deXxV9Z3/8dcnudkT1iQQcoMBARFByQLihlsdtdiCYqidrmOn1pl2apffTGt/8+vMdGYe/XV+M9o63dzaaR9trWwCUqzaghUsgglhkUWJBCELkAABkkCWm+/vj3vAgEgSspzce9/PxyOPnPM959x8ch/kfQ/f8z3fY845REQkusT5XYCIiPQ9hbuISBRSuIuIRCGFu4hIFFK4i4hEoYDfBQBkZma6/Px8v8sQEYkoZWVl9c65rPNtGxThnp+fT2lpqd9liIhEFDN794O2qVtGRCQKKdxFRKKQwl1EJAop3EVEopDCXUQkCincRUSikMJdRCQKRXS4Vxw6wXee30Fre4ffpYiIDCoRHe77jjTzs9cqWb3roN+liIgMKhEd7rMnZpGdkcSi0iq/SxERGVQiOtwD8XHcXZjLK2/XcejEKb/LEREZNCI63AFKivIIdTie21TtdykiIoNGxIf7hOx0CsYOY1FZFXoerIhIWMSHO4TP3isONbJ5f4PfpYiIDApREe53XZVDckIci8p0YVVEBKIk3IckJ3Dn1Bye31LDqbaQ3+WIiPguKsIdoKQoyIlT7by4/YDfpYiI+C5qwn3W+JEEh6dozLuICFEU7nFxxvzCIK+9U091w0m/yxER8VXUhDvAvUVBnIMlurAqIjEuqsI9b0Qq14wfyeKyKjo6NOZdRGJXVIU7QElxkH1Hmtm494jfpYiI+Cbqwv3OqTmkJwV0YVVEYlrUhXtKYjx3XZnDqm21NLa0+12OiIgvoi7cIdw1c7ItxKqttX6XIiLii6gM98Kxwxmflcaisv1+lyIi4otuh7uZxZtZuZmt9NbHmdkGM6sws2fNLNFrT/LWK7zt+f1U+4Vq5d6iIG/sPUplfdNA/3gREd/15Mz9IWBnp/XvAY865yYAR4HPee2fA4567Y96+w24+YVB4gwW6+xdRGJQt8LdzILAHOApb92AW4DF3i6/AOZ5y3O9dbztt3r7D6hRQ5KZPSmLJWXVhDTmXURiTHfP3L8P/APQ4a2PBBqcc6eHo1QBud5yLrAfwNt+zNv/LGb2gJmVmllpXV3dxVXfhZKiPA4cP8W6ivp+eX0RkcGqy3A3s7uAQ865sr78wc65J5xzxc654qysrL586TM+NCWbYakJLCxV14yIxJbunLlfB3zUzPYCvyXcHfMDYJiZBbx9gsDph5hWA3kA3vahwOE+rLnbkgLxzL1qDC9vP0hDc6sfJYiI+KLLcHfOPeycCzrn8oH7gNXOuU8Aa4B7vd0+Ayz3lld463jbVzsfH25aUpxHa6iDFVtq/CpBRGTA9Wac+zeAr5lZBeE+9ae99qeBkV7714Bv9q7E3rlizBAuzxmi6QhEJKYEut7lPc65V4BXvOU9wMzz7HMKKOmD2vqEmVFSFOQ7K3ew68BxJo8e4ndJIiL9LirvUD3XvIJcEuJNZ+8iEjNiItxHpCVy6+RRLCuvpi3U0fUBIiIRLibCHcKTiR1uamX1rkN+lyIi0u9iJtxvnJRFVkaSumZEJCbETLgH4uO4pyCXNW8dou5Ei9/liIj0q5gJdwh3zYQ6HMvKq7veWUQkgsVUuE/IzmB63jAWle3Hx/uqRET6XUyFO4TP3t8+2MjWqmN+lyIi0m9iLtw/ctUYkgJxekqTiES1mAv3IckJ3DF1NCs213CqLeR3OSIi/SLmwh3C87wfP9XOSzsO+l2KiEi/iMlwv/bSkeQOS2GR5nkXkSgVk+EeF2fML8xlXUU9NQ0n/S5HRKTPxWS4A9xblIdzsHST7lgVkegTs+E+dmQqV48bweKyKo15F5GoE7PhDrCgOI+9h5t5Y+9Rv0sREelTMR3ud04bTXpSQBdWRSTqxHS4pyYGmDMth99tq6Wppd3vckRE+kxMhzuEpyNobg2xalut36WIiPSZmA/3okuGMz4zjUVlGjUjItEj5sPdzJhfFGRj5RH21jf5XY6ISJ+I+XAHmF8YJM5gsc7eRSRKKNyB0UOTuWFiFks2VRHq0Jh3EYl8XYa7mSWb2UYz22Jm283sX7z2/zGzSjPb7H1N99rNzB4zswoz22pmhf38O/SJkuIgtcdO8VpFvd+liIj0WqAb+7QAtzjnGs0sAVhnZi942/7eObf4nP3vBCZ6X1cDP/G+D2ofunwUQ1MSWFRWxexJWX6XIyLSK12eubuwRm81wfu6UN/FXOCX3nGvA8PMLKf3pfav5IR45k4fw4vbD3Csuc3vckREeqVbfe5mFm9mm4FDwMvOuQ3epn/3ul4eNbMkry0X6HzLZ5XXdu5rPmBmpWZWWldXd/G/QR8qKcqjtb2DFVtr/C5FRKRXuhXuzrmQc246EARmmtlU4GFgMjADGAF8oyc/2Dn3hHOu2DlXnJU1OLpBpuYOYfLoDBZrOgIRiXA9Gi3jnGsA1gB3OOdqva6XFuDnwExvt2ogr9NhQa9t0DMz7i0KsqXqGG8fPOF3OSIiF607o2WyzGyYt5wC3AbsOt2PbmYGzAPe9A5ZAXzaGzUzCzjmnIuYe/vvLsglEGeaTExEIlp3ztxzgDVmthV4g3Cf+0rg12a2DdgGZAL/5u2/CtgDVABPAn/b51X3o5HpSdwyOZvnyqtpC3X4XY6IyEXpciikc24rUHCe9ls+YH8HfLH3pfmnpDiPl3Yc5JW36rhtyii/yxER6THdoXoeN12WRWZ6orpmRCRiKdzPIyE+jnsKg6zedYj6xha/yxER6TGF+wcoKQrS3uFYVh4RA31ERM6icP8AE0dlcFXeMD1AW0QiksL9AkqKguw6cII3q4/7XYqISI8o3C/gI1eNISkQx6IyXVgVkciicL+AoSkJ3H7FaJZvruFUW8jvckREuk3h3oWS4iDHTrbxh50H/S5FRKTbFO5duPbSTMYMTWZRqR7BJyKRQ+Hehfi48AO01+6u48CxU36XIyLSLQr3bri3KEiHgyWbdPYuIpFB4d4Nl4xMY+a4ERrzLiIRQ+HeTSVFQSrrmyh796jfpYiIdEnh3k0fnpZDamI8CzWZmIhEAIV7N6UlBZgzLYffba2lubXd73JERC5I4d4DJcV5NLWGWLXtgN+liIhckMK9B2bkDyd/ZKrmeReRQU/h3gOnH6C9ofII+w43+12OiMgHUrj30D2FQcxgsSYTE5FBTOHeQ2OGpXDDxCyWbKqmo0Nj3kVkcFK4X4SSoiDVDSf58zuH/S5FROS8FO4X4bYpoxiSHNA87yIyaHUZ7maWbGYbzWyLmW03s3/x2seZ2QYzqzCzZ80s0WtP8tYrvO35/fw7DLjkhHjmTs/l928e4NjJNr/LERF5n+6cubcAtzjnrgKmA3eY2Szge8CjzrkJwFHgc97+nwOOeu2PevtFnZLiIC3tHazcWuN3KSIi79NluLuwRm81wftywC3AYq/9F8A8b3mut463/VYzs74qeLCYljuUy0ZlaJ53ERmUutXnbmbxZrYZOAS8DLwDNDjnTt+HXwXkesu5wH4Ab/sxYOR5XvMBMys1s9K6urpe/RJ+MDNKioNs3t9AxaETfpcjInKWboW7cy7knJsOBIGZwOTe/mDn3BPOuWLnXHFWVlZvX84X8wpyCcSZzt5FZNDp0WgZ51wDsAa4BhhmZgFvUxCo9pargTwAb/tQICrHDGamJ3Hz5GyWllfTHurwuxwRkTO6M1omy8yGecspwG3ATsIhf6+322eA5d7yCm8db/tqF8VPuCgpClJ3ooU/vR15XUsiEr26c+aeA6wxs63AG8DLzrmVwDeAr5lZBeE+9ae9/Z8GRnrtXwO+2fdlDx43T84mMz1RXTMiMqgEutrBObcVKDhP+x7C/e/ntp8CSvqkugiQEB/H/MIgT6zdwzMb9/HxmWP9LklEpOtwl6599bZJvH3wBA8v3UbjqXY+P3u83yWJSIzT9AN9IDkhnsc/VcycK3P491U7eeSlt/QgbRHxlc7c+0hiII7H7isgIynAY6srOH6qnW/fNYW4uKi7f0tEIoDCvQ/FxxnfvWca6UkBnlpXSWNLO//3nmkE4vUfJBEZWAr3PmZm/O85l5ORnMCjf3ibppZ2vn/fdJIC8X6XJiIxRKeU/cDMeOhDE/n2XVN44c0D/PUvSmlube/6QBGRPqJw70f3Xz+O/5h/Ja9V1PPppzdy/JSmBxaRgaFw72cLZuTx3x8vZEtVAx9/4nUON7b4XZKIxACF+wCYc2UOT366mHfqGlnw+Hpqj530uyQRiXIK9wFy02XZ/PL+qzl4vIV7f7KevfVNfpckIlFM4T6AZo4bwTOfn0Vzazslj6/nrQOaB15E+ofCfYBNCw5l4ReuIc5gwePr2by/we+SRCQKKdx9MHFUBosfvJahKQl84snX+fM79X6XJCJRRuHuk7wRqSx68BrGDEvhsz9/gz/sOOh3SSISRRTuPho1JJlnv3ANk0dn8OCvyli+ubrrg0REukHh7rMRaYn8+q+vpuiS4Xzl2c38ZsM+v0sSkSigcB8EMpIT+MX9M7n5smy+9dw2Hv/TO36XJCIRTuE+SCQnxPPTTxZx15U5fPeFXfzni5oTXkQunmaFHEQSA3H84L4C0pMC/HBNBY0tmhNeRC6Own2QOT0nfEZygCfXVnLiVDvfm6854UWkZxTug5CZ8a0Ph+eEf+Tlt2lsaeOxjxdoTngR6TadDg5SZsaXbw3PCf/i9oOaE15EeqTLcDezPDNbY2Y7zGy7mT3ktf+zmVWb2Wbv68OdjnnYzCrM7C0zu70/f4Fod//14/iPe8Nzwn/q6Y0cO6k54UWka905c28Hvu6cmwLMAr5oZlO8bY8656Z7X6sAvG33AVcAdwA/NjP1J/TCguI8fviXhWz15oSv15zwItKFLsPdOVfrnNvkLZ8AdgK5FzhkLvBb51yLc64SqABm9kWxsezD08Jzwu+pD88JX9OgOeFF5IP1qM/dzPKBAmCD1/QlM9tqZj8zs+FeWy6wv9NhVZznw8DMHjCzUjMrraur63nlMej0nPB1x1so+el6KjUnvIh8gG6Hu5mlA0uArzjnjgM/AS4FpgO1wH/15Ac7555wzhU754qzsrJ6cmhMmzluBM88MIuTbSFKfrqeXQeO+12SiAxC3Qp3M0sgHOy/ds4tBXDOHXTOhZxzHcCTvNf1Ug3kdTo86LVJH5maO5SFX5hFfBx87PHXKd931O+SRGSQ6c5oGQOeBnY65x7p1J7Tabe7gTe95RXAfWaWZGbjgInAxr4rWQAmZHeaE/6pDZoTXkTO0p0z9+uATwG3nDPs8T/MbJuZbQVuBr4K4JzbDiwEdgC/B77onAv1T/mx7fSc8MHhmhNeRM5mg2FyquLiYldaWup3GRHraFMrn/n5RrbXHOf/zLmcj80YS0qiRp+KRDszK3POFZ9vm+5QjQLDvTnhrx43gn9+fgcz/v0PPLx0G+X7jmpmSZEYpTP3KOKcY2PlERaWVrFqWy0n20JMzE5nQXEedxfmkpme5HeJItKHLnTmrnCPUidOtbFyay0LS/dTvq+BQJxx6+XZLCjO48ZJWZplUiQKKNxj3O6DJ1hUVsXSTVXUN7aSnZHE/KIgJUVBxmel+12eiFwkhbsA0BbqYPWuQywq3c+at+oIdThm5A+npDiPOdNySEvSDNAikUThLu9z6PgplpZXs/CN/eypbyItMZ67rhzDghlBCscOJ3x7g4gMZgp3+UDOOcrePcrC0v2s3FpLc2uI8VlpLCjO457CXLIzkv0uUUQ+gMJduqWppZ3fbatlUel+3th7lPg44+bLsllQHOTmydkk6CKsyKCicJcee6eukUWlVSzZVEXdiRYy05OYX5hLSXGQCdkZfpcnIijcpRfaQx386e06nn1jP6t3HaK9w1E4dhgLivOYc2UOGckJfpcoErMU7tIn6k60sKy8mmdL91NxqJGUhHg+PC2Hj83IY0a+LsKKDDSFu/Qp5xyb9zewsHQ/z2+ppbGlnfyRqZQU5zG/MMjooboIKzIQFO7Sb5pb23lh2wEWlu5nQ+UR4iz8xKgHb7yUmeNG+F2eSFRTuMuA2FvfxKKy/SwqraKusYW/vn4cX/+Ly0hO0AyVIv1Bs0LKgMjPTOPvb5/MK39/E5+4eixPrq3koz9cx5vVx/wuTSTmKNylz6UmBvi3edP4n7+aQUNzG3f/+DV+tKaC9lCH36WJxAyFu/Sbmy7L5sWvzOYvrhjN/3vxLRY8vp699U1+lyUSExTu0q+GpyXyo78s5Af3TafiUCN3/mAtv3r9XT1ERKSfKdxlQMydnstLX72R4vzh/OOyN/nsz9/g4PFTfpclErUU7jJgRg9N5pf3z+Q7c69gQ+Vhbv/+q6zcWuN3WSJRSeEuA8rM+PQ1+az68g1cMjKNL/2mnC8/U86x5ja/SxOJKgp38cX4rHSWPHgNX7ttEqu21XL7919l7e46v8sSiRoKd/FNID6OL986kef+9jrSkwN86umNfHv5m5xsDfldmkjE6zLczSzPzNaY2Q4z225mD3ntI8zsZTPb7X0f7rWbmT1mZhVmttXMCvv7l5DINi04lJV/dz33XzeOX65/lzmPraV831G/yxKJaN05c28Hvu6cmwLMAr5oZlOAbwJ/dM5NBP7orQPcCUz0vh4AftLnVUvUSU6I59sfmcJvPn81p9pC3PvT9Tzy0lu06cYnkYvSZbg752qdc5u85RPATiAXmAv8wtvtF8A8b3ku8EsX9jowzMxy+rpwiU7XXprJ7786m7nTx/DY6gru/vFr7D54wu+yRCJOj/rczSwfKAA2AKOcc7XepgPAKG85F9jf6bAqr+3c13rAzErNrLSuThfS5D1DkhN4ZMF0fvrJQmoaTjHnv9fx9LpKOjp045NId3U73M0sHVgCfMU5d7zzNhe+3bBHf3nOuSecc8XOueKsrKyeHCox4o6pOfz+Kzdww4RM/nXlDj7x1AaqG076XZZIROhWuJtZAuFg/7VzbqnXfPB0d4v3/ZDXXg3kdTo86LWJ9Fh2RjJPfaaY782fxtaqBu549FWWlFVp+gKRLnRntIwBTwM7nXOPdNq0AviMt/wZYHmn9k97o2ZmAcc6dd+I9JiZ8bEZY3nhodlMzsng64u28De/2sThxha/SxMZtLp8WIeZXQ+sBbYBp4cufItwv/tCYCzwLrDAOXfE+zD4IXAH0Az8lXPugk/i0MM6pLtCHY6n1u7hv156myEpCXxv/jRuvXxU1weKRCE9iUmizs7a43z12c3sOnCC+2bk8Y93TSE9KeB3WSIDSk9ikqhzec4Qln/pOv7mpktZWLqfO3/wKhsrj/hdlsigoXCXiJUUiOcbd0xm4ReuwTA+9sR6vvvCTlraNX2BiMJdIl5x/gheeOgG7psxlsf/tIe5P3yNHTXHuz5QJIop3CUqpCUF+O490/jZZ4upb2xl7o/W8eNXKnQWLzFLF1Ql6hxpauUfl21j1bYDpCTEM3PcCG6YmMkNE7OYNCqd8IAukcin0TISc5xzrN1dzx93HmTt7nr2eA/mzs5I4voJmVw/MfyVnZHsc6UiF+9C4a6xYxKVzIzZk7KYPSk8tUV1w0nW7a5j7e561rx1iKXl4ZumJ4/OOBP2V48bSUpivJ9li/QZnblLzOnocOyoPc6ru+tYt7ue0r1HaQ11kBgfR3H+cK6fmMkNE7K4YswQ4uLUhSODl7plRC7gZGuIjXuPnDmz33UgPMXw8NQErp2QyeyJmVw/MYvcYSk+VypyNnXLiFxASmI8N07K4kavC+fQiVO8VlHP2t31rNtdz++2hqdGGp+ZFj6rn5jFrPEjyEhO8LNskQvSmbvIBTjn2H2okVffrmNdRT0b9hzhZFuI+DijIG+YF/aZXBUcRiBeI4tlYKlbRqSPtLSH2PRuA+sqwv31W6uP4RxkJAW45tKR3OB14eSPTNWQS+l3CneRfnK0qZU/v3OYdRXh/vqqo+GHiQSHp4SDfkIW100YybDURJ8rlWikcBcZAM453j3czFrvwuz6dw5zoqWdOIPrJ2YxvzCX268YTXKChltK31C4i/igPdTBlqpjrN51kGXlNVQ3nCQjKcCcK3OYXxSk+JLh6rqRXlG4i/iso8PxeuVhlpRV88KbtTS3hhg7IpV7CnOZXxgkb0Sq3yVKBFK4iwwiTS3t/P7NAyzZVMX6PYdxDmaOG8G9hUHunDZaQyyl2xTuIoNUdcNJnttUxZJN1VTWN5GcEMftV4xmfmGQ6yZkEq87ZOUCFO4ig5xzjvL9DSwpq+L5LTUcP9XO6CHJzCvI5d6iXCZkZ/hdogxCCneRCHKqLcTqXYdYUlbFK2/XEepwXBUcyj2FQT561RiGp2lYpYQp3EUiVN2JFpZvrmbJpmp21h4nId64ZXI28wuD3HRZNokB3RUbyxTuIlFgR81xlm6qYtnmGuobWxiRlshHrxrD/MIgU3OHaFhlDOpVuJvZz4C7gEPOuale2z8DnwfqvN2+5Zxb5W17GPgcEAK+7Jx7sasCFe4i3dce6uDV3XUsKavm5R0HaQ11MGlUOvMLg8wryGXUED2AJFb0NtxnA43AL88J90bn3H+es+8U4BlgJjAG+AMwyTl3wQdZKtxFLs6x5jae31rD0k1VbNrXoLthY0yvpvx1zr1qZvnd/Flzgd8651qASjOrIBz067tbrIh039DUBD456xI+OesS9tQ1snRTNUs3VfHQbzefuRv2nsIgM/J1N2ys6c187l8ys08DpcDXnXNHgVzg9U77VHlt72NmDwAPAIwdO7YXZYgIwPisdP7X7Zfxtdsm8fqewyzeVMWKLTX89o39Z+6Gvbsgl0tGpvldqgyAbl1Q9c7cV3bqlhkF1AMO+Fcgxzl3v5n9EHjdOfcrb7+ngRecc4sv9PrqlhHpH+e7G7Zg7DDuLshlzrQcRqYn+V2i9EKfP4nJOXew04s/Caz0VquBvE67Br02EfFBWlKA+UVB5hcFqWk4yYotNTy3qZpvL9/Od57fwY2TsphXkMuHLh+lh4NHmYsKdzPLcc7Veqt3A296yyuA35jZI4QvqE4ENva6ShHptTHDUnjwxkt58MZL2Vl7nGXl1SzfXMMfdx0iPSnAHVNHM296LtdcOlLTHkSBLsPdzJ4BbgIyzawK+CfgJjObTrhbZi/wBQDn3HYzWwjsANqBL3Y1UkZEBt7lOUO4PGcI/3DHZDZUHmZZeTUvbDvA4rIqsjOSmDt9DPMKcpmSo/HzkUo3MYkIEJ724I87D/FceTV/evsQbSHHpFHpzCvIZe70XHKHpfhdopxDd6iKSI8cbWrld9tqWVZeTem7R4HwtMR3F+Ty4ak5DE3VtMSDgcJdRC7avsPNLN9czXObq9lT10RifBy3TM5mXkEuN0/OIimgC7F+UbiLSK8559hWfYxl5TWs2BKe32ZIcoA5V45h3vQxzMgfQZwuxA4ohbuI9Kn2UAevvRO+EPvi9gM0t4bIHZbC3OljuLsgl4mjNP/8QFC4i0i/aW5t5+UdB3muvJq1u+sJdTiuGDOEuwty+chVYzSRWT9SuIvIgKg70cLKrTUsK69mS9Ux4gyuvTSTeQW53DF1NOlJvZnxRM6lcBeRAfdOXSPLy6tZtrmGfUeaSU6I47Ypo5k9MZOM5AApiQFSE+NJSYgnNTGe1MQAKd66HkLSPQp3EfGNc45N+xpYVl7Nyq01HG1u6/KYQJyRktgp9L0PgNNtKQnxZz4czrQnvPcB8d4+4bbOxyYH4qPmwq/CXUQGhbZQB9VHT9LcGuJkWzvNraHw8unvbSFOtp7Tfk7bqbZQp+3tNLeF6GmMpSTEk5WRxLjMNMZnpTE+M43xWemMz0pj9JDkiLkrt88nDhMRuRgJ8XHkZ/btlMPOOVraO7zAbz/ngyL0Xvs5HwoHjrewp66RN/Yeobn1vVlSUhLi3xf6p9czkiPn5i2Fu4hENDMjOSGe5IR4RqQl9vh45xwHvaB/p76Jyrom9tQ3sq36GKu21dLR6X8Fp8/2L81KCwd+ZvhsP29EKgnxg+s6gcJdRGKamTF6aDKjhyZz7YTMs7a1tIfYd7iZPfVN7KlrorK+kT11Tby4/SBHmlrP7BeIM8aOSD1zhj/OC/3xWWlkpSf50s2jcBcR+QBJgXgmjso4701ZDc2t7wv9yvom1lXU09LecWa/jKQA48450z/9IZCa2H8RrHAXEbkIw1ITKRybSOHY4We1d3Q4qhtOUlnfxJ66xvD3+iZK9x5l+eaas/YdPSSZz10/js/PHt/n9SncRUT6UFyckTcilbwRqcyelHXWtlNtISrrm84E/566JrKH9M+jDhXuIiIDJDkh/syDUvrb4Lq8KyIifULhLiIShRTuIiJRSOEuIhKFFO4iIlFI4S4iEoUU7iIiUUjhLiIShQbFfO5mVge8e5GHZwL1fVhOpNP7cTa9H+/Re3G2aHg/LnHOZZ1vw6AI994ws9IPmqw+Fun9OJvej/fovThbtL8f6pYREYlCCncRkSgUDeH+hN8FDDJ6P86m9+M9ei/OFtXvR8T3uYuIyPtFw5m7iIicQ+EuIhKFIjrczewOM3vLzCrM7Jt+1+MnM8szszVmtsPMtpvZQ37X5DczizezcjNb6XctfjOzYWa22Mx2mdlOM7vG75r8YmZf9f5G3jSzZ8ws2e+a+kPEhruZxQM/Au4EpgAfN7Mp/lblq3bg6865KcAs4Isx/n4APATs9LuIQeIHwO+dc5OBq4jR98XMcoEvA8XOualAPHCfv1X1j4gNd2AmUOGc2+OcawV+C8z1uSbfOOdqnXObvOUThP94c/2tyj9mFgTmAE/5XYvfzGwoMBt4GsA51+qca/C1KH8FgBQzCwCpQE0X+0ekSA73XGB/p/UqYjjMOjOzfKAA2OBzKX76PvAPQIfPdQwG44A64OdeN9VTZpbmd1F+cM5VA/8J7ANqgWPOuZf8rap/RHK4y3mYWTqwBPiKc+643/X4wczuAg4558r8rmWQCACFwE+ccwVAExCT16jMbDjh/+GPA8YAaWb2SX+r6h+RHO7VQF6n9aDXFrPMLIFwsP/aObfU73p8dB3wUTPbS7i77hYz+5W/JfmqCqhyzp3+n9xiwmEfiz4EVDrn6pxzbcBS4Fqfa+oXkRzubwATzWycmSUSviiywueafGNmRrhPdadz7hG/6/GTc+5h51zQOZdP+N/FaudcVJ6ddYdz7gCw38wu85puBXb4WJKf9gGzzCzV+5u5lSi9uBzwu4CL5ZxrN7MvAS8SvuL9M+fcdp/L8tN1wKeAbWa22Wv7lnNulX8lySDyd8CvvROhPcBf+VyPL5xzG8xsMbCJ8AizcqJ0GgJNPyAiEoUiuVtGREQ+gMJdRCQKKdxFRKKQwl1EJAop3EVEopDCXUQkCincRUSi0P8HsP/BeOM8FwoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s[:10])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_top_words=8\n",
    "\n",
    "def show_topics(a):\n",
    "    top_words = lambda t: [vocab[i] for i in np.argsort(t)[:-num_top_words-1:-1]]\n",
    "    topic_words = ([top_words(t) for t in a])\n",
    "    return [' '.join(t) for t in topic_words]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['critus ditto propagandist surname galacticentric kindergarten surreal imaginative',\n",
       " 'jpeg gif file color quality image jfif format',\n",
       " 'graphics edu pub mail 128 3d ray ftp',\n",
       " 'jesus god matthew people atheists atheism does graphics',\n",
       " 'image data processing analysis software available tools display',\n",
       " 'god atheists atheism religious believe religion argument true',\n",
       " 'space nasa lunar mars probe moon missions probes',\n",
       " 'image probe surface lunar mars probes moon orbit',\n",
       " 'argument fallacy conclusion example true ad argumentum premises',\n",
       " 'space larson image theory universe physical nasa material']"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(Vh[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        We get topics that match the kinds of clusters we would expect! This is despite the fact that this is an unsupervised algorithm - which is to say, we never actually told the algorithm how our documents are grouped."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-negative Matrix Factorization (NMF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Idea\n",
    "        Rather than constraining our factors to be orthogonal, another idea would be to constrain them to be non-negative. NMF is a factorization of a non-negative data set V:\n",
    "            V = WH\n",
    "\n",
    "        into non-negative matrices W, H. Often positive factors will be more reasily iterpretable.\n",
    "        \n",
    "        NMF is a non-exact factorization that factors into one skinny positive matrix and one short positive matrix. NMF is NP-hard and non-unique. There are number of variations on it, created by adding different constraints.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NMF from Sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "m,n = vectors.shape\n",
    "d = 5 # number of topics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = decomposition.NMF(n_components  = d, random_state = 135)\n",
    "W1 = clf.fit_transform(vectors)\n",
    "H1 = clf.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['jpeg image gif file color images format quality',\n",
       " 'edu graphics pub mail 128 ray ftp send',\n",
       " 'space launch satellite nasa commercial satellites year market',\n",
       " 'jesus god people matthew atheists does atheism said',\n",
       " 'image data available software processing ftp edu analysis']"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(H1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TF-IDF\n",
    "\n",
    "        Topic Frequency-Inverse Document Frequency is a way to normalize term counts by taking into account how often they apear in a document, how long the document is and how common/rare they are in a document.\n",
    "\n",
    "        TF = (# occurences of terms t in document) / (# of words in documents)\n",
    "\n",
    "        IDF = log(# of documents / # documents with term t in it)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorizer_tfidf = TfidfVectorizer(stop_words='english')\n",
    "vectors_tfidf = vectorizer_tfidf.fit_transform(newsgroups_train.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\Monit Sharma\\AppData\\Local\\Programs\\Python\\Python37\\lib\\site-packages\\sklearn\\decomposition\\_nmf.py:294: FutureWarning: The 'init' value, when 'init=None' and n_components is less than n_samples and n_features, will be changed from 'nndsvd' to 'nndsvda' in 1.1 (renaming of 0.26).\n",
      "  FutureWarning,\n"
     ]
    }
   ],
   "source": [
    "W1 = clf.fit_transform(vectors_tfidf)\n",
    "H1 = clf.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['people don think just like objective say morality',\n",
       " 'graphics thanks files image file program windows know',\n",
       " 'space nasa launch shuttle orbit moon lunar earth',\n",
       " 'ico bobbe tek beauchaine bronx manhattan sank queens',\n",
       " 'god jesus bible believe christian atheism does belief']"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(H1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2644cd38d68>]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlEklEQVR4nO3de5gcZZ0v8O+PYNhHxCOY6PEBJMiJx40rgo7IWS/rDQmoiXgjeFnw6HL0mF336OoGkIgBBBFQgeASSBBBCAEkDEkkJCFALkwyk9tkMpPJTCbXyWUmM5OZZO498zt/dHVS01PV/VZ3VVd19ffzPHnSXV1T9VZX9a/eeq+iqiAiouJ3StgJICIifzCgExHFBAM6EVFMMKATEcUEAzoRUUycGtaOx40bpxMmTAhr90RERWnDhg1HVHW802ehBfQJEyagqqoqrN0TERUlEdnj9hmLXIiIYoIBnYgoJhjQiYhiggGdiCgmGNCJiGKCAZ2IKCYY0ImIYsIooIvIZBGpF5FGEZnhss43RKRWRLaJyJP+JpOI8rVpbwdqmjvDTgYFKGvHIhEZA2A2gMsA7AdQKSLlqlprW2cigBsAfExVO0TkHUElmIhyc9WDawEAu+/8QsgpoaCY5NAvAdCoqk2qOgBgPoCpaev8C4DZqtoBAKra4m8yiYgoG5OAfjaAfbb3+61ldu8F8F4RWSMiFSIy2WlDInK9iFSJSFVra2tuKSYiIkd+VYqeCmAigE8BuAbAwyLytvSVVHWOqpapatn48Y5jyxARUY5MAnozgHNt78+xltntB1CuqoOqugvADiQDPMWUquLupfXYfaQ77KQQkcUkoFcCmCgi54vIWADTAJSnrbMQydw5RGQckkUwTf4lk6Jmf0cvHljZiO/+qTLspBCRJWtAV9UEgOkAlgKoA7BAVbeJyCwRmWKtthRAm4jUAlgJ4Geq2hZUoik6BoeGw04CEVmMxkNX1SUAlqQtm2l7rQB+Yv0jIqIQsKcoEVFMMKATUckZHlZ09Q2GnQzfMaATUcn53fIduPCWl9HRPRB2UnzFgE5EJWdx9UEAQHsPAzoREUUQAzoRUUwwoBMRxQQDOhFRTDCgExHFBAM6EVFMMKATEcUEAzoRUUwwoBNRyVINOwX+YkAnotIjYScgGAzoREQxwYBO5EFNcyc0bs/pFBsM6ESGXt/Rii/evxpPVOwJOylEjhjQiQztae8BAGw/dCzklBA5Y0AnIooJBnQiophgQCciigkGdCKimGBAJyKKCQZ0IqKYMAroIjJZROpFpFFEZjh8fp2ItIrIZuvf9/1PKhERZXJqthVEZAyA2QAuA7AfQKWIlKtqbdqqT6vq9ADSSEREBkxy6JcAaFTVJlUdADAfwNRgk0VERF6ZBPSzAeyzvd9vLUv3VRGpFpFnReRcpw2JyPUiUiUiVa2trTkkl4iI3PhVKfoigAmqeiGAZQAec1pJVeeoapmqlo0fP96nXRMREWAW0JsB2HPc51jLTlDVNlXtt94+AuDD/iSPiIhMmQT0SgATReR8ERkLYBqAcvsKIvIu29spAOr8SyJRtHDwXIqqrK1cVDUhItMBLAUwBsA8Vd0mIrMAVKlqOYB/E5EpABIA2gFcF2CaiUIR00luKEayBnQAUNUlAJakLZtpe30DgBv8TRoREXnBnqJEVMLiVYDGgE5EJSeuxWcM6EREMcGATkQUEwzoREQxwYBORBQTDOhERDHBgE5EFBMM6EREMcGATkQUEwzoREQxwYBORBQTDOhEHmm8hv+gGGFAJzIkcR0AhGKDAZ2IKCYY0ImIYoIBnYgoJhjQiYhiggGdiCgmGNCJiGKCAZ2IKCYY0ImoZMWtkxgDOhGVHIlpLzGjgC4ik0WkXkQaRWRGhvW+KiIqImX+JZGIiExkDegiMgbAbABXAJgE4BoRmeSw3hkAfgxgnd+JJCKi7Exy6JcAaFTVJlUdADAfwFSH9W4F8BsAfT6mj4iIDJkE9LMB7LO9328tO0FEPgTgXFVd7GPaiCIqZjVpFBt5V4qKyCkA7gXwU4N1rxeRKhGpam1tzXfXRAUliGdFGsWHSUBvBnCu7f051rKUMwD8A4BXRWQ3gEsBlDtVjKrqHFUtU9Wy8ePH555qIiIaxSSgVwKYKCLni8hYANMAlKc+VNVOVR2nqhNUdQKACgBTVLUqkBQTEZGjrAFdVRMApgNYCqAOwAJV3SYis0RkStAJJCIiM6earKSqSwAsSVs202XdT+WfLCIi8oo9RYmIYoIBnYgoJhjQiYhiggGdiCgmGNCJiGKCAZ2IKCYY0ImIYoIBnYgoJhjQiYhiggGdiCgmGNCJqGTFbWR7BnQiKjlxHdmeAZ2IKCYY0ImIYoIBnchHQ8OKf7xjBcq3HAg7KVSCGNCJfNQ9kMCBzj7c9NetYSeFShADOhFRTDCgExHFBAM6EVFMMKATBSBuHVaoODCgE/mo0B1Wbl5YgwkzFhd4rxRVDOhERezxij1hJ4EihAGdiCgmGNCJPFIWkFNEGQV0EZksIvUi0igiMxw+/4GIbBWRzSKyWkQm+Z9UonCJhwJyZdSnEGQN6CIyBsBsAFcAmATgGoeA/aSqfkBVLwJwF4B7/U4oUTEQL1GfyGcmOfRLADSqapOqDgCYD2CqfQVV7bK9PR1stUVEVHCnGqxzNoB9tvf7AXw0fSUR+RGAnwAYC+AzThsSkesBXA8A7373u72mlYiIMvCtUlRVZ6vqBQD+E8AvXNaZo6plqlo2fvx4v3ZNREQwC+jNAM61vT/HWuZmPoAv55EmoqLHMkcKg0lArwQwUUTOF5GxAKYBKLevICITbW+/AKDBvyQSFQ9WiVKYsgZ0VU0AmA5gKYA6AAtUdZuIzBKRKdZq00Vkm4hsRrIc/dqgEkxE8bCvvQdzV+8KNQ1xa11qUikKVV0CYEnaspm21z/2OV1EFHPfnrsOe9p6cNXFZ+Os08cWdN9xbV3KnqJEAYhbzi8Ix/oSANgJy08M6EQ+imvOj4oDAzoRUUwwoBNF0P6OHtQfOhZ2MgqCBS7+MaoUJaLC+vhvVgIAdt/5hZBTEhyWTvmPOXQij0zq8JT5TgoBAzqRIXuOctPeDpTdtgydPYNp6zDfSeFhQCfKwX0rGnDk+AA27G0POylFj60W/cOATkShYBNP/zGgExHFBAM6UQBKtRihuz+B7v5E2MkoWWy2SOSjUi9GeP8vlwLw1tySLYL8wxw6EYWkxO9+AWBAp6xWNxzBpb9egd6BobCTQkQZMKBTVr9eUodDXX3Y2Xo87KQQUQYM6EQBYKmwB/yyfMOATkShKPUK5CAwoBMRxQQDOhGRR4uqD+DTd7+K4eFolRexHToVpcGhYQDAm8YUPk9i0m6apQnmwgyJubaB/9kz1egdHEJfYghvHhudMMocOhWlj9y+HBfe8nJB9+mlzDda+bZoCvOmlxoVczChSFiZgzhgQKeidLRnEL2D4beLL9Uu/nHxpQdW4/O/ez3sZPiGAZ0oB8ImGrHRdKQ77CT4hgGdqACO9yfwyKomKLP0o8TpK1nX1Ia5q3eFtn+jgC4ik0WkXkQaRWSGw+c/EZFaEakWkRUicp7/SSUqImlB6tYXa3Hb4jqsrG8JJz0Bueze1/C9P1Xm9LdxfMi5ek4Fbl1UG9r+swZ0ERkDYDaAKwBMAnCNiExKW20TgDJVvRDAswDu8juhRMWsszc5VV3/YHwq4ACgoeU4VmyP102qmJnk0C8B0KiqTao6AGA+gKn2FVR1par2WG8rAJzjbzKJ4iFGpQvGGluO45FVTWEnoySYBPSzAeyzvd9vLXPzPQB/c/pARK4XkSoRqWptbTVPJVGRi2PxgqmrHlyD2xbXuTYP5Hjo/vG1UlREvg2gDMBvnT5X1TmqWqaqZePHj/dz10SREKcKPr+kZjBKbxkk7H7lO5MuTs0AzrW9P8daNoKIfA7ATQD+SVX7/UkeUXGZX7kXADDglhst4YCfbOHDIB4kkxx6JYCJInK+iIwFMA1AuX0FEbkYwEMApqgqa0ioZC2vO+y4vJSLXNza7LOoxX9ZA7qqJgBMB7AUQB2ABaq6TURmicgUa7XfAngLgGdEZLOIlLtsjqgklXLOvNj0DQ6hLwK9kHNhNKqMqi4BsCRt2Uzb68/5nC4iipn0e1pUy9Dfd/NLOPPNb8KmmZ8POymesacokaF8ctklXeQSdgJy0NEzaLRe1J68GNApJweO9oadhEiKaq6T/BXVGzQDeg56B4awwqXyq1RcPacCAJAYilgWJUB+/IijWBF4sLMXbceDb5jmlpuNWi63mDGg5+DmF2rwvceqUHugK+ykFJX7VzTghc2jWryWhCjn3P/XHa/gw7ctD2z7bjfCqOZyi1l0ptooInvaksNtHrc6TJCZe5btAABMvShTR+PiZg9Sr+1oxd+/6wy844y/Cy9BFLp97T0YGlZMGHd64PtiDj0PHAqVMrl23np847/eGLGslC+ZKBY3BW3l9hZ84q6V+NTdrxZkfwzoOYjy4zP5p3dgCM9t2J/XjXt3W0/2lWLO6ffyT79diYOdfQDiPWBZoYfSZUAncjFrUS1++swWVDS1h52UWLDfF/eEfKMLo/y+9kAXXtxyINB9MKDnwS1nceR4P74zdx06ugcKmh7yV0tXMgfZ7UddSSk/1JXysdtced8q/OtTmwLdBwN6LrJcoHNX78KqhiN4cv3ewqSHqEBUFWt3Hin5+qOoHj4Deh6ielLJX36e5tS2VBUVTW1FFxifWr8P33x4HV6sPujbNoP6DtqO92Nfe7BFO/aimyiM/8KAnoNsT5CF/I2u3XkEP/rLxtACQym2XMhF+jWzcHMzps2pwHMbi6td/p72ZJPd5g7znsJhlbhcescKfOKulYFsu9cheN/4162B7MsLBvQ8ZAtmhah4uW5eJRZvPeg6/jblz8tpdBsqNl2qUnBvwDnIUjYYUC/mVQ3Os61t3NsRyP68YEDPAXu4kW+KrMglH3E51PpDx8JOgisG9CJRd7ALE2YsxuZ9Rx0/j8uPpVREvS+DquL/PF6F13Y450a9FLVlywDx2vUPA7qh9u4BPLKqaWRZtdtgQ4YX+8r6FkyYsRiNLcezrvvK9uREUEu3HRr5QbTjQixlCkBxOR3DCizddhjffXT9iOX2G5FpvU2cA7b92EyL24LEsVwM/XTBZqysb0XZhLOMc1fZ1ltstRTYuLcD/+MdbzHaZpx/HFHndD69no/0IFjMp7M84E4y5B1z6Ia6+pKdSxK2ykfXH6PXX6nB+q4j1nncVVSoKm7469ZIVCRl43R63L53t1M5asb7EE/cfSsafNlO6zGzIXcjkHEtGQzoOfDrAs1lM27FOcWWc+8dHMJT6/fiWw+vCzspoQrjvN1rjXqZj/nr93n+m7g0cY1C0YobBvSQrGk8gmc27AdgdqG7Fd8U4tpKDCefSvZ7aHtsKi4/8hS305G+fGtzZ9BJCUTqevPS3NKvCuCXag7haE90h9OIQphnQM9DttxVpmA7d/Uuf9MSYGDccThZaXvT8/51nIh6Kw8/OQWhZbXJGa/idkMLSktXH37wxAb88ImNYScl0lcuA3oATH6i9ovC02N32rrFHhiLragok8TQMI72jp5cOKgOLsUmn3Pdn0g+JUahI1aES1wY0IOU6bzbLwqjG0CEL6JcpI4nTqHuFwtrsMWhn0Dczp1X2Y7/E3etRFff6Buhn7r7E+gZSKD5aC86HW66vojAeTYK6CIyWUTqRaRRRGY4fP5JEdkoIgkR+Zr/yYwOhfkPVAEsqNyHj935isOnuZ19twBYTDndA0f9L4uPglzGui6m8wYEF7M27T2aeb957vj9v1yKi361DB+78xV87t7XMDSc+xcfgbjtKmtAF5ExAGYDuALAJADXiMiktNX2ArgOwJN+JzAKZr1Yiw17Mjev23agE1+6fzV6BkaOnf3z56rR7BDARuTQTZotui3P4er6zD2v4tMFmhLLyaB93JkiC2i5MCleq2nuxBMVewqSHjfzVu/CjsPBdWvPdKpT7fMfeKUBE2YszrpeLlLjHbUe68dwsd1JDZnk0C8B0KiqTao6AGA+gKn2FVR1t6pWA4jlCFHz1jhXYNortO5Ysh1bmztHBP5MsfYU24c3Pr911NyTbgaHhrG4+qBrB5VvPlyB69J696Vrau3GriPdRvuz86voQNVe5OL9h9XZE+zjuRczX6gxWs/tu0sd/RfvX41fLDTbVlBmLarFF+9bnXGdXMKgyWWT2u7dLzs3qUw1FTzQ2TdicvbE0DAWVO7DsJXjvnVRLS7/3es5pNIlXap4tb7F6EYShZy7SUA/G4C90el+a5lnInK9iFSJSFVrq/MYEcXAqSIydcdXdc9F1B3swpdnr0HPQGLUNtbvzjzNWSogPLpmN3705Ea8VHPISstIa3e24dX6YL5bPzM1uVbm7mw9jg/Oetm/hGTgFoS7+gaxwhqK4YA1L2ZcDAwNo3dgCBfcuCTspIxgPxXf+1Plidfz1uzCz5+rxtNVyRA1d/Uu1Pv4lPH8pmZc92glnrK1u2c7dIuqzlHVMlUtGz9+fCF3HYhUgEvO4tKWfG37PP283764Dpv3HcWGPR1553bbApzerndgCOt3mc2jebjLrLegXUVTGzZZPUS93iR2tY58smg51ocJMxZjvjU71LceqcCFtyz1nKZsVta3YNBqj79022FPf2sPAMXQTLEjQ1vvfC7bjLlcD1+L/Sk49Ts4GtBTW6q+p/losnXN3NW78LrrgGXhMxnLpRnAubb351jLSlZ6MM6lfsVrQM+neeL2Q10Y95bTMO4tpxmtP+Ov1Xhh8wGs+vmnce5Zb855v+7bP9meXZH8gZ5+2hi877+/1fO2dh9J/tCe27gfF737bVjT2Oa43qLqA+joGcSahiP4xHvH4VsfPc94H2t3HsGja3Zj/Blm31+66ObnnPkdmExytApFTS6drWyJTXiYEyDXXPati2o9rd+UQ9FmPkxy6JUAJorI+SIyFsA0AOXBJqs4OF34JmVtqvm3H1ckZxHvHhjKut/Jv1+Fz97zmvG2tx9MPrJ2p1XwBvGkOTSs+Oof12Ly71fl9Pf24860jelPbsLNC2vw0rZDuOl5b2XVR44nc4FuY5esamjFpb9egb5E9oDi9bz/8IkNntZ/bsP+vMfHsX+nQUy/5+aL948uv39kVRP2tfcYXXs9AU4B5/bzGvFEHtjezWUN6KqaADAdwFIAdQAWqOo2EZklIlMAQEQ+IiL7AXwdwEMisi3IREfNiFYMcD/59ovScw49bf2DR3tx5X3mQTCwtrcFlv49vG7NHhNkB6tsW759cR0OdfUZNYVLL3LJdv//W82hzCuk+ekzW/CVB9d6+pugGVWKavp7RUtXH25bXIdrH10fWge6KJeXOzEqQ1fVJar6XlW9QFVvt5bNVNVy63Wlqp6jqqer6ttV9f1BJjrq1lnlz5meAPO9UB58dafnv/nXpzZ5Wr8YWnbttqZyC/L3nu1UmZzLfJP3tT+uxQ0e5qy0twTxyu28v1rfghbDERbz3afqyaLM7v5E1o54QcfdIvgpAGBP0awyFWWkPhuxhp4ceOnIceeLX+F//DG54Ew7vkQlU7JxbwdW1reMWDYqbdaBuxWHrNzekrFdsx+yfV2Z+hyYVJI+XbkXVXs68NT6vejsHTQqK15c7f9Y5dc9WolnrQHlcuFXBmFoWPOeBs7tt5mPKOTmGdCz+H9Pbx61LNOJM23F8EbT6Mq7/sTQqI5JpewrD67Fdx+tzLhO6vt2a1e/cLNz/f3wsBp3UvESsJ3/Xjz92NNz1//53Mmc+Qd/9fKISuWUY31mgT4UBoeefibS36dv4vLfJ9uab9l/NLm+esskBdUqJmwM6Fks3Oye08kWDtxGVLx23nrHHOUVf1iFSTOdm9xla0boR+7n3mU78Pc3v5T/hgoo03Fv2NOOF1zO33tuXGJchJEtGHvNmGUrDnnglcaMny/cNPom9YFbXsZPn9niLSEuolDUZnqzrWhK/i6Czgi51osFulfvGNBzkOkkZpxvMssvv6nVvYnTy7Xe2j7n4r4VDei1tRSIwg87XXrlWKYu3Fc/VJFxW/MrT3YWOd6fwOqGIy77zOyULOfV/tSmCtz/SoP9w1FyzWm73by8CqytfKZm6FnHovY3KaYiUIriCQN6TFQadgRyk88YGWHKVJbq5Yh+8vRmfHvuupwGDjP5zaeesF7b0YrBRHF+17ky+X4ee2P3iPdev6EojM0ShdjPgO7RiOsm/GvohO//uSqvvz/UdbIL+542tzGno3DJjrTbNa3eblINLclJPHqd2jLnW4gOnBig7dX6lhFDxeZyCSWGFR/99XK8sj2/p7ZsA2EBmQ+93bC3sskxrmsamSFJP3XFPu5/oTCg5+Bw18jxO0Z0xHC5eoeH1bXLcBTY0+0Y1JAcw+SCG5dgVUPhj+Mjty/Hh29d5ulvTAKJSVDKFkyyhRr7d9vVl8irpUjK4a5+3Lqozmjdqt3taGwZ/STzuMvojqb3wUwjM+5zmIgiVZRzyGH8m/Qctr3Y53BXP9q6M7dKiUAGPRJMuv5Tmu3WY/6BztGP507LAOCV7S2Oy6PCpKxwwOoJma3SLgipSuTq/Wbdw9u7B4x+5B09Azjr9LH5JC2vctZ8irpM//ZrhiN5et+/+2d9Vqags3dwxA1v/a52fOMhs/TUHeo68XrGc5krsIOO50bz/kbgIYI59DwMDytuXlhzoukUAMx8wbmTbL9Bt/BiEXSv06sfegMvuDQ3fHWH2Y3RdFjbVIVmKjjm8pvcmuUmE1SwyXe7bk8efqRXAayoO4wP/upldPUlW6CoJkccdZLeyVYVWFx98MT7zQ4zQaWvn49vPuxcgZ7t6WzmC9sC7+fgBQN6HnoGhvB4xR5cM2dd2EnJm9OF65Yr2Z5WEen3TOzrdrXjx/M3O36WbWabFNMb6Clphz3lgTUAgBc2N2N5XfLmka3CLZHH7DcPr9qF3Wlt6E23lm8Qcx2j3Yfyi2HVEz2mc9n2ztbjo85NJgrNq2NParRUr57bmH/xmZ8Y0D1qOTa6/K8YhkR1ssjWm9Ctp6WJi2YtQ38iuIGRcmE6xdhLNYcw/cmTM8mn2ojbbyjPO7T79tO/5FihXcgJk00C8YDtJurW0cf03vebl+qzNge1e+i1JvzDL/0dNnnCjMUnhnkulp84A7pH0588OR7KiVl3IniyH1nV5LjcPtuP/ViW1XobBCrdQFqOeNPejtyGQ/WJaZ3FHX/bjkW2R/sgZMvhp1rYpBSqKNZtP7nMZgUA5bahJVSB2rTiFYWXp4/8ctzp28pVqv/H85uaPc2UNTSsrsWGQWJAz0OqmMKkDWyu1+bQsOKWcu+DV/5heYPjctfZfvL88aTnvK56cK3jcKil6LG1u8NOgie/W35yGrjUaXW6xNMXpXeIWpXWUatvcMg4uDYf7UWDTzMP3WxYn9LS5T77VMuxfvz70+6D26Uf1zcfrnAtNgwSA7oP8ihCzWrtziP4Uw4BwTRJB4724qUa5xyqp4xNhnW7+xO4zePEAIVmL3Ja5nOvXKcmfFGWa4bW/jtwKoa86Xnz0SKbWrtRlWVi9mHDH94TFXuN1vt6ltY3mWbnSk+KU/1BIbDZYh6CLjsv33IA/+ZxyFu7TXs7Tsw96uYrD67Foa4+/PizE3PeD5B8Slm4qRmNLcfxH5f/zxGf/WJhTV7l0LNeDP5mkJooBMi9TNtNPkPZ5mPDng7XViWA+1AUTnHS5Eq3P6k63RTqDh7DP14wzmBLZn6/wvkpNLl/77/NPW09OfUUBpJP0lEoei25gL5+VzvKzjsTp3ipQg/Jk+ucO36YON6fwFUGEx2keojm27LhYlunn/SAbhrMD3f14Z1v/btRy+etcR7krFikWssU2lf/mNtEF+kd5QYSwycrB10MJIaNriE/u+jflyGg5+pXL+Y2N8+waiSGHyipIpdVDa34xkNvYI5LhaFXXroju3end3bbotoTI8mZOtaXe/vw+xw6C1XtbseVf8htargULzeKj/56RV77iovmHHOJdtnG91pQtc94Pz97dguunjO6nbb93F47b/2IlkVOp/1QVx8WVO0b/UEAdhw+nn0lB5la1qRPyWinGo3GESUV0A9aXY4bW3I72em83JF/89J243XvXlqPR1yG3s3kL+vMygpN3fJi7aiWCl6df8MSn1JTOrxOO+fkxizl1T9/ttr1s/R+Bm6jOK6sPzkExBtNbVnL0AFgZ4YRRf2UGi/dq0wBPVOmbEij0Xi56AL60Z4BLM+x0ip1svx6NLrjb+ZB2osHVubWtf7OgNJTSFMeKO2WMRe7tUIqAt22uoL0Fi7FIp/WaKZ9H4JUdAH9h09sxPf/XJXTFFKnRLjdOCWZjtUSVx1FPJPOPctONnf87dL6EFOSOy+dmexU1bijl2nrnFwUXUDf05Z8ZKtp7vTcAcKeQ992oBMTZizGOoep4IDk2BHPFKi8L26qbWPbEBWToRxze15y5++5cUlgufmia+WSap1ynTXX5O47v2D0d1v2HT3Re3BYgTWNyUfC5XWH8dH3vH3U+l+evcaP5Jak1HgoRMVmcY69hk3buqcMDg1jzCljctpXJkWXQx+ToblhYmgYj1fscZzCa+rsNSe6Jg8P64kKHKe2uE5jRxMRubH3rjUR1HDaRRfQ08u4jvUNYsnWg2g51ocfPLEBNy+sOTE5s6pi15HuUY83XX2DJ8rR57zedKJLcqps63P35lZDTkRkIp/B8DIxKnIRkckA/gBgDIBHVPXOtM9PA/BnAB8G0AbgalXd7W9Sk9Iz6B+4ZXSrgM37jqJ8ywG8vO0QFlUfxMXvftuIz1c1HEHtgZPN8d5nm+n+li9N8jW9RETp8hlyORPJ1vFDRMYA2AHgMgD7AVQCuEZVa23r/F8AF6rqD0RkGoCrVPXqTNstKyvTqirvXawvu/e1UaPTEREVk7NOH4uNN1+W09+KyAZVLXP6zKTI5RIAjarapKoDAOYDmJq2zlQAj1mvnwXwWfFr7Ms0mcrQiYiKgekE216ZBPSzAdjb7+23ljmuo6oJAJ0ARjUdEZHrRaRKRKpaW3ObaPhbl56X098REUXFPV//YCDbLWizRVWdA2AOkCxyyWUb37n0PHyHQZ2IaBSTHHozgHNt78+xljmuIyKnAvhvSFaOEhFRgZgE9EoAE0XkfBEZC2AagPK0dcoBXGu9/hqAV9SPmWaJiMhY1iIXVU2IyHQAS5FstjhPVbeJyCwAVapaDmAugMdFpBFAO5JBn4iICsioDF1VlwBYkrZspu11H4Cv+5s0IiLyouh6ihIRkTMGdCKimGBAJyKKCQZ0IqKYyDqWS2A7FmkFkOu09uMAFOccV96UwnHyGOOhFI4RiMZxnqeq450+CC2g50NEqtwGp4mTUjhOHmM8lMIxAtE/Tha5EBHFBAM6EVFMFGtAnxN2AgqkFI6TxxgPpXCMQMSPsyjL0ImIaLRizaETEVEaBnQiopgouoAuIpNFpF5EGkVkRtjp8UpEdovIVhHZLCJV1rKzRGSZiDRY/59pLRcRuc861moR+ZBtO9da6zeIyLVu+ysEEZknIi0iUmNb5tsxiciHre+s0frbgs9D6HKMt4hIs3UuN4vIlbbPbrDSWy8il9uWO16/1vDU66zlT1tDVReUiJwrIitFpFZEtonIj63lcTuXbsdZ/OdTVYvmH5LD9+4E8B4AYwFsATAp7HR5PIbdAMalLbsLwAzr9QwAv7FeXwngbwAEwKUA1lnLzwLQZP1/pvX6zBCP6ZMAPgSgJohjArDeWlesv70iIsd4C4D/cFh3knVtngbgfOuaHZPp+gWwAMA06/V/AfhhCMf4LgAfsl6fgeTk8JNieC7djrPoz2ex5dBNJqwuRvZJth8D8GXb8j9rUgWAt4nIuwBcDmCZqrarageAZQAmFzjNJ6jq60iOg2/nyzFZn71VVSs0+ev4s21bBeNyjG6mApivqv2qugtAI5LXruP1a+VSP4PkBOvAyO+rYFT1oKputF4fA1CH5HzBcTuXbsfppmjOZ7EFdJMJq6NOAbwsIhtE5Hpr2TtV9aD1+hCAd1qv3Y63GL4Hv47pbOt1+vKomG4VN8xLFUXA+zG+HcBRTU6wbl8eGhGZAOBiAOsQ43OZdpxAkZ/PYgvocfBxVf0QgCsA/EhEPmn/0Mq5xKotaRyPyfJHABcAuAjAQQD3hJoan4jIWwA8B+DfVbXL/lmczqXDcRb9+Sy2gG4yYXWkqWqz9X8LgOeRfGw7bD2Owvq/xVrd7XiL4Xvw65iardfpy0OnqodVdUhVhwE8jOS5BLwfYxuSxRWnpi0vOBF5E5JB7i+q+ldrcezOpdNxxuF8FltAN5mwOrJE5HQROSP1GsDnAdRg5CTb1wJ4wXpdDuCfrdYElwLotB59lwL4vIicaT0Wft5aFiW+HJP1WZeIXGqVTf6zbVuhSgU5y1VInksgeYzTROQ0ETkfwEQkKwMdr18r17sSyQnWgZHfV8FY3+9cAHWqeq/to1idS7fjjMX5LETNq5//kKxZ34Fk7fJNYafHY9rfg2RN+BYA21LpR7LMbQWABgDLAZxlLRcAs61j3QqgzLat/41k5UwjgO+GfFxPIfmIOohkeeH3/DwmAGVI/rh2AngAVg/nCBzj49YxVCP5o3+Xbf2brPTWw9aSw+36ta6N9daxPwPgtBCO8eNIFqdUA9hs/bsyhufS7TiL/nyy6z8RUUwUW5ELERG5YEAnIooJBnQiophgQCciigkGdCKimGBAJyKKCQZ0IqKY+P8qGT8aw1GR1QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(clf.components_[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43.7129260580871"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.reconstruction_err_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NMF in summary\n",
    "\n",
    "    Benefits : fast and easy to use\n",
    "\n",
    "    Downside : took years of research and expertise to create\n",
    "\n",
    "    Notes:\n",
    "        1. For NMF, matrix needs to be at least as tall as it is wide, or we get error with fit_transform\n",
    "        2. Can use df_min in CountVectorizer to only look ar words that were in at least k of the split texts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NMF from scratch in  NumPy ,  using Stochastic Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        The key idea of standard gradient descent \n",
    "        1. Randomly choose some weights to start\n",
    "        2. Loop:\n",
    "           1. Use Weights to calculate a prediciton\n",
    "           2. Calculate the derivative of the loss\n",
    "           3. Update the weights\n",
    "        3. Repeat step 2 lots of times, eventually end up with some descent weights.\n",
    "        \n",
    "        We want to decrease our loss and the derivative tells us the direction of steepest descent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic Gradient Descent (SGD)\n",
    "\n",
    "    It is an incredibly useful optimization method ( it is also the heart of deep learning, since used for backpropoagation)\n",
    "\n",
    "    We evaluate the loss using all our data for standard gradient descent, which is slow.\n",
    "\n",
    "    In stochastic gradient descent , we evaluate our loss function on just a sample of our data. We would get differet loss values on different samples of the data. It turns out this is still effective way to optimize.\n",
    "\n",
    "\n",
    "#### Applying SGD to NMF\n",
    "        Goal: Decompose V(m * n) into \n",
    "                V = WH\n",
    "            where W and H >=0 , and we've minimized the Frobenius norm of V- WH\n",
    "\n",
    "        Approach: We will pick random positive W & H , and then use SGD to optimize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "lam = 1e3\n",
    "lr = 1e-2 \n",
    "m, n = vectors_tfidf.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\Monit Sharma\\AppData\\Local\\Programs\\Python\\Python37\\lib\\site-packages\\sklearn\\utils\\validation.py:598: FutureWarning: np.matrix usage is deprecated in 1.0 and will raise a TypeError in 1.2. Please convert to a numpy array with np.asarray. For more information see: https://numpy.org/doc/stable/reference/generated/numpy.matrix.html\n",
      "  FutureWarning,\n",
      "c:\\Users\\Monit Sharma\\AppData\\Local\\Programs\\Python\\Python37\\lib\\site-packages\\sklearn\\decomposition\\_nmf.py:294: FutureWarning: The 'init' value, when 'init=None' and n_components is less than n_samples and n_features, will be changed from 'nndsvd' to 'nndsvda' in 1.1 (renaming of 0.26).\n",
      "  FutureWarning,\n"
     ]
    }
   ],
   "source": [
    "W1= clf.fit_transform(vectors)\n",
    "H1 = clf.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['jpeg image gif file color images format quality',\n",
       " 'edu graphics pub mail 128 ray ftp send',\n",
       " 'space launch satellite nasa commercial satellites year market',\n",
       " 'jesus god people matthew atheists does atheism said',\n",
       " 'image data available software processing ftp edu analysis']"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(H1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 1e-6\n",
    "def grads(M, W, H):\n",
    "    R = W@H-M\n",
    "    return R@H.T + penalty(W, mu)*lam, W.T@R + penalty(H, mu)*lam # dW, dH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def penalty(M, mu):\n",
    "    return np.where(M>=mu,0, np.min(M - mu, 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def upd(M, W, H, lr):\n",
    "    dW,dH = grads(M,W,H)\n",
    "    W -= lr*dW; H -= lr*dH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def report(M,W,H): \n",
    "    print(np.linalg.norm(M-W@H), W.min(), H.min(), (W<0).sum(), (H<0).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = np.abs(np.random.normal(scale=0.01, size=(m,d)))\n",
    "H = np.abs(np.random.normal(scale=0.01, size=(d,n)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "44.42710988660707 1.4206941402729504e-06 6.220086649529336e-08 0 0\n"
     ]
    }
   ],
   "source": [
    "report(vectors_tfidf, W, H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "upd(vectors_tfidf,W,H,lr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "44.419184642449935 -0.0007765716529699559 -7.190500098920772e-05 134 284\n"
     ]
    }
   ],
   "source": [
    "report(vectors_tfidf, W, H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "44.41318263157269 -0.0008430536728499327 -7.167740576098002e-05 130 261\n",
      "44.375136145889066 -0.0003328731899605962 -4.81782057093232e-05 52 522\n",
      "44.34628347718999 -0.0002643278292372564 -6.728869170547078e-05 41 958\n",
      "44.31425263334867 -0.00014943471301808218 -9.658154475246471e-05 26 1518\n",
      "44.279008899552316 -0.0001003564775838509 -0.00012961678356928944 25 2244\n"
     ]
    }
   ],
   "source": [
    "for i in range(50): \n",
    "    upd(vectors_tfidf,W,H,lr)\n",
    "    if i % 10 == 0: report(vectors_tfidf,W,H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['god space like people don just think know',\n",
       " 'space people don just think know god like',\n",
       " 'space god don people think just know like',\n",
       " 'god don people space just like think know',\n",
       " 'god space don people just know like think']"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        This is painfully slow to train! Lots  o fparameter fiddling and still slow to train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyTorch\n",
    " \n",
    "        It is a Python framework for tensors and dynamic neural networks with GPU accelaration. Similar to NumPy , only with increases parallelization of using GPU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        If you are not using a GpPU , you will need to remove the .cuda() from this method below                                ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch as tc\n",
    "from torch.autograd import Variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "def V(M): return Variable(M, requires_grad=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "v=vectors_tfidf.todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_vectors = torch.Tensor(v.astype(np.float32))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 1e-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grads_t(M, W, H):\n",
    "    R = W.mm(H)-M\n",
    "    return (R.mm(H.t()) + penalty_t(W, mu)*lam, \n",
    "        W.t().mm(R) + penalty_t(H, mu)*lam) # dW, dH\n",
    "\n",
    "def penalty_t(M, mu):\n",
    "    return (M<mu).type(tc.FloatTensor)*torch.clamp(M - mu, max=0.)\n",
    "\n",
    "def upd_t(M, W, H, lr):\n",
    "    dW,dH = grads_t(M,W,H)\n",
    "    W.sub_(lr*dW); H.sub_(lr*dH)\n",
    "\n",
    "def report_t(M,W,H): \n",
    "    print((M-W.mm(H)).norm(2), W.min(), H.min(), (W<0).sum(), (H<0).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_W = tc.FloatTensor(m,d)\n",
    "t_H = tc.FloatTensor(d,n)\n",
    "t_W.normal_(std=0.01).abs_(); \n",
    "t_H.normal_(std=0.01).abs_();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "d=6; lam=100; lr=0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(44.3427) tensor(-0.0051) tensor(-0.0004) tensor(783) tensor(1398)\n",
      "tensor(43.7378) tensor(-0.0065) tensor(-0.0100) tensor(1473) tensor(18494)\n",
      "tensor(43.7144) tensor(-0.0064) tensor(-0.0042) tensor(1633) tensor(19758)\n",
      "tensor(43.7058) tensor(-0.0044) tensor(-0.0055) tensor(1688) tensor(19458)\n",
      "tensor(43.7021) tensor(-0.0055) tensor(-0.0052) tensor(1856) tensor(19207)\n",
      "tensor(43.7011) tensor(-0.0041) tensor(-0.0040) tensor(2211) tensor(20863)\n",
      "tensor(43.7006) tensor(-0.0034) tensor(-0.0062) tensor(2395) tensor(21466)\n",
      "tensor(43.7003) tensor(-0.0041) tensor(-0.0039) tensor(2256) tensor(22396)\n",
      "tensor(43.7001) tensor(-0.0043) tensor(-0.0055) tensor(2204) tensor(22454)\n",
      "tensor(43.6999) tensor(-0.0039) tensor(-0.0050) tensor(2513) tensor(28230)\n"
     ]
    }
   ],
   "source": [
    "for i in range(1000): \n",
    "    upd_t(t_vectors,t_W,t_H,lr)\n",
    "    if i % 100 == 0: \n",
    "        report_t(t_vectors,t_W,t_H)\n",
    "        lr *= 0.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['people don think just objective like say morality',\n",
       " 'god jesus bible believe atheism christian belief does',\n",
       " 'graphics comp software group 3d just aspects amiga',\n",
       " 'thanks files file image program windows format know',\n",
       " 'space nasa launch shuttle orbit moon lunar earth']"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(t_H.cpu().numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x26483662748>]"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAl+klEQVR4nO3de3gU9b0/8PcHEPsc66m20nNasQUt7fnZak9t6vH8jm09rRfUHtGjrehzvPUcObby9O6v4K0W73driygqKt4Q71GiCMpFBCQBwiVAQgiBJEASCCEJkMtmP78/djZMNjM7M7uzO7OT9+t5eNidncx+Z2fmM9/5XkVVQUREhW9I0AkgIiJ/MKATEUUEAzoRUUQwoBMRRQQDOhFRRAwL6ouPOeYYHTVqVFBfT0RUkFauXLlbVUdYfRZYQB81ahTKysqC+noiooIkItvsPmORCxFRRDCgExFFBAM6EVFEMKATEUUEAzoRUUS4CugiMlZEKkWkWkQm2azzMxHZICIVIvKSv8kkIiInjs0WRWQogKkAzgJQD6BURIpVdYNpnTEAJgP4N1XdKyJfzFWCiYjImpsc+qkAqlW1RlW7AcwCMC5lnWsBTFXVvQCgqk3+JpOIsrV6+16sb9gXdDIoh9wE9GMB1Jne1xvLzL4O4Osi8omILBeRsVYbEpEJIlImImXNzc2ZpZiIMnLRY0vxk78uCToZlEN+VYoOAzAGwBkALgPwpIgclbqSqk5X1SJVLRoxwrLnKhERZchNQG8AcJzp/UhjmVk9gGJV7VHVrQCqkAjwRESUJ24CeimAMSIyWkSGAxgPoDhlnbeQyJ1DRI5Bogimxr9kEhGRE8eArqoxABMBzAWwEcBsVa0QkSkicoGx2lwAe0RkA4AFAG5Q1T25SjQREQ3karRFVS0BUJKy7FbTawXwO+MfEREFgD1FiYgiggGdiCgiGNCJiCKCAZ2IKCIY0ImIIoIBnTLS2dOLUZPm4O3y1D5mRBQUBnTKSHN7FwDg/rmVAaeEiJIY0ImIIoIBnYgoIhjQiYgiggGdiCgiGNCJiCKCAZ2IKCIY0ImIIoIBnYgGne5YHJsb24NOhu8Y0Ilo0LnlrfU46+HFaGrvDDopvmJAJ6JBp7S2BQDQ3hkLOCX+YkAnIooIBnQioohgQCciiggGdCKiiGBAJyKKCAZ0IqKIYEAnIooIVwFdRMaKSKWIVIvIJIvPrxaRZhEpN/79j/9JJSKidIY5rSAiQwFMBXAWgHoApSJSrKobUlZ9RVUn5iCNREQ5oRp0CvzlJod+KoBqVa1R1W4AswCMy22yiIhySIJOQG64CejHAqgzva83lqW6WETWishrInKc1YZEZIKIlIlIWXNzcwbJJSIiO35Vir4DYJSqngxgHoDnrFZS1emqWqSqRSNGjPDpq4mICHAX0BsAmHPcI41lfVR1j6p2GW+fAvBdf5JHRERuuQnopQDGiMhoERkOYDyAYvMKIvIl09sLAGz0L4lEROSGYysXVY2JyEQAcwEMBTBDVStEZAqAMlUtBvArEbkAQAxAC4Crc5hmIiKy4BjQAUBVSwCUpCy71fR6MoDJ/iaNiIi8YE9RIqKIYEAnIooIBnQioohgQCciiggGdCKiiGBAJ3KpoyuGCTPL0NTWGXRSiCwxoBO59ObqBnywoRF/+XBz0EkhssSATkQUEQzoREQRwYBORBQRDOhERBHBgE5Eg1i05qBjQCeiQSeiM9AxoBMRRQUDOhFRRDCgExFFBAM6EVFEMKATEUUEAzqRR9Fq6EZRwoBO5FJUm7pRdDCgExFFBAM6EVFEMKATEUUEAzoRUUS4CugiMlZEKkWkWkQmpVnvYhFRESnyL4lEROSGY0AXkaEApgI4F8CJAC4TkRMt1jsSwK8BfOp3IomIyJmbHPqpAKpVtUZVuwHMAjDOYr3bAdwLgDPoEhEFwE1APxZAnel9vbGsj4icAuA4VZ2TbkMiMkFEykSkrLm52XNiiYjIXtaVoiIyBMBDAH7vtK6qTlfVIlUtGjFiRLZfTUREJm4CegOA40zvRxrLko4E8C0AC0WkFsBpAIpZMUpElF9uAnopgDEiMlpEhgMYD6A4+aGq7lPVY1R1lKqOArAcwAWqWpaTFBMRkSXHgK6qMQATAcwFsBHAbFWtEJEpInJBrhNIFDbK0bkiI2rHcpiblVS1BEBJyrJbbdY9I/tkEYWPcHSuyJCIHkz2FCUiiggGdCKiiGBAJyKKCAZ0IqKIYEAnIooIBnQioohgQCciiggGdCKiiGBAJyKKCAZ0IqKIYEAnIooIBnQizyI2ohNFBgM6kUuCaA7oRNHBgE5EFBEM6EREEcGATkQUEQzoREQRwYBORINW1NorMaAT0aAT1fZKDOhERBHBgE5EFBEM6EQF7OUV2zH2kcVBJ4NCYljQCSCizE1+Y13QSaAQYQ6dyGePLazGluaOoJNBg5CrgC4iY0WkUkSqRWSSxefXicg6ESkXkSUicqL/SSUKv/bOHtz3fiUufWJ50EmhQcgxoIvIUABTAZwL4EQAl1kE7JdU9SRV/WcA9wF4yO+EEhWCZLvmrp7eQNNBg5ObHPqpAKpVtUZVuwHMAjDOvIKqtpneHoHotdcnIgo9N5WixwKoM72vB/AvqSuJyPUAfgdgOIAfWW1IRCYAmAAAX/nKV7ymlYiI0vCtUlRVp6rqCQD+COBmm3Wmq2qRqhaNGDHCr68mCh0+olIQ3AT0BgDHmd6PNJbZmQXgwizSRFSwotqlnAqDm4BeCmCMiIwWkeEAxgMoNq8gImNMb88HsNm/JBIRkRuOZeiqGhORiQDmAhgKYIaqVojIFABlqloMYKKInAmgB8BeAFflMtFERDSQq56iqloCoCRl2a2m17/2OV1EBU2VpehOKne14/VV9Zh87j9BhIVVfmBPUSIfMTC5d/mTyzF9cQ327O8OOimRwYBORIFIPsPwFugfBnQioohgQCfyiMXjFFYM6EQueSkezzbmb9jRhqVbdme5lcIQ5P0xajdnBnSiDLQe6MbUBdWIx/tHBL/Kg8979GNc/uSnPm0tnIIsO49q3TUDOlEGbn5rPe6fW4kl1YMjF02FgQGdKAP7u2IAgFg8HnBKCl/Uij2CxIBOlAMMUs6iWuwRJAZ0Ih8xSFGQGNCJyDdFd8zH6fd+FHQyBi1XY7kQEbmxu6PL898oR4/3DXPoRBQQlk/5jQGdHHXFelFe1xp0MgoKc50UBAZ0cnTrWxW4cOonqGs5EHRSQk+Y66QAMaCTo3UN+wAA+w72BJwSiiQ+zPiGAZ3II7Yx9webePqPAZ3IJS/xh0GfgsCATuQj5jopSAzoRBQoPsz4hwGdKAssWslcIT/MLNm8Gz9/tnTA8MlBY09RcsRihIGcJoMO12VOfrt2ZhkO9vSiM9aLvxsenjDKHDoVpHkbGvH++l1BJ4MoVFwFdBEZKyKVIlItIpMsPv+diGwQkbUi8qGIfNX/pBIdcu3MMlz3wsqgk0E+CLLYasaSrShesyO4BPjMMaCLyFAAUwGcC+BEAJeJyIkpq60GUKSqJwN4DcB9fieUiKIlyKK8ZI/eV8rq8KuXVweXEJ+5yaGfCqBaVWtUtRvALADjzCuo6gJVTfYLXw5gpL/JJCosqbGqoyuGp5dshbIWlXLITUA/FkCd6X29sczOfwN4z+oDEZkgImUiUtbc3Ow+lUQFJjVs3/7OBtz+7gYsqGwKJD25sriqGau27w06GWTwtVJURP4LQBGA+60+V9XpqlqkqkUjRozw86spD5i5dNbY1gkAiPX2n2s0OQ5OV0+05iC9csYK/OdjS7PaBkem9I+bgN4A4DjT+5HGsn5E5EwANwG4QFW9j3JPoWVV1lnd3JH/hIREugB081vrAQAha54cqMa2TsytGNgiKYojU27dvR+LqoIrfXAT0EsBjBGR0SIyHMB4AMXmFUTkOwCeQCKYR+uZkixd80wpAKCnN1o5znTcVOLF+RgzwKVPLMP/Pr8SvRG8y6Ue7n9/YCGumrEimMTARUBX1RiAiQDmAtgIYLaqVojIFBG5wFjtfgCfBfCqiJSLSLHN5ihiopjLIn9ttxlHP6xFLfsO9qCjK5Z2nbB2tnPVxUlVSwCUpCy71fT6TJ/TRQUirBdlWPHXCr9v//kDDBsiqL7rvKCT4hl7ihL5yO6JJaw5uiCF+ekuVqDFQwzo5Bpz45lj0TrlAwM6OUqXk/ISqOJxHfQda8KbJ6UoYECnvDn+xhJcPC27NssUPYP8Hu8rBvQMdMV6UVbbEnQyQsHrtbhqe2sukhEahVhWvrCyCeV1rXn/3kL8rcKOAT0Dd7y7EZc8vgybG9uDTgoVmDBmRq9+phQXTv0k598z2Ivb8oEBPQObdrUBAPYe6Ak4JRRm3bF434w2gzk36jQZSCFzukWp6oBhIHKJAT0LzHGw/DOdr9/8Hm4yhgIge4V4Crm9RU1btAVfu+k97MtT5o8BPQNhbj9L/onHte9pDMjs5vXyiu0Z/23UDYar6LWyegBA7Z79efk+BvQcKK9rxahJc1CxY1/QSaEs/G1BNcY+8jHWN/h3HAdDEKOBxuWhjgJgQM+KXaYrObLcwsrBMOZ7dLOea4yWH7v2JYbE9aMoOLq/FoUBA3omHC5sPl5TqgjXC2ZtsNRFdfb05rwsnQE9C07nYT4u4h2tB/Hm6vrcf5GNQXItuua2RUdvXPHyiu15bQHhh3hcMX9DY0ZB2PwXD8+rwg7jyScI2V6bVnu/xWKOAPN6lz6xDN+e8kF2X+yAAT0DYcpsXTp9GX77yhp0xXqDTkrk5PJe9cLybZj8xjrMXLYth9/ivxc+3Yb/mVmGN1cPmOPGltX18pcPN/uXKBvPfrIVk99Yl5NtH+hOXG/mfbv2ubK0f7OmPvd1agzoWQjDYFWN+/I3OdRgzY3n4klr74FuAEDrwcLqy9DQehAA0NgW/knJbntnQ18rIz8tr9ljuTwMk5swoGcgiPJQVcUHFbv6OqoM/Dx3351uf4M/hcPF86kRgiAQtEL7CcytnsKWdAb0HMhFzv31VQ2Y8PxKvPBpyiN6mMp/yLUo9GU42M1iPrMw9IhlQM+GU6Wow0Wrqq7bOCdnk98ZYEWSn+JxxUWPfYJ5GxqDTkpWUnOXXm/lYcvhuaVQvJSD4gzKDgO6Swe7e7FgU2L+a8fclYurVFVx93ub8JO/LsGyLdZlctZ/1/998HmCzHTGerF6eyt+9fLqoJOSkWx/9yAzcz9+cKHjOqqKUZPm4P65m/otN5/7blu6FOpNqxAxoLt0y9vrcc2zpf27gmexvaeXbMX0xTUAgLq91pPomjkFgEIrh4yqQrjBbml27oaerKqZtnCLb99bCL9NNsKwfwzoLm3dnbgIOjpjvuSuPqnefeiNh2CcWj4fdLFdtp1CwtBSKEiFeiOu2uV96OhsdzUsM155KSvPd3oZ0HPI7XF3E9TCUInm56kZhv0JQvL6LvS9f6t8h+t1/djXjq4Yjr+xBI/5+MQQRQzoWbC7+boJfBnXiNt+Z+5yArkMPiHIcNmyy11lkma7w11oTyhBPRHu3Z9ot//Sp8FXxNr+BCG4SzOge+Tl8nN7fN0ECLsLKeicbibhqPVAd+BFRV4k0+r0Wy+sbMKiKusB2QZUZhfI/ufjdvPx5t3OK4VImI+dq4AuImNFpFJEqkVkksXnPxCRVSISE5FL/E9m8DI5hopE2bufuQrzBfb8sloc7Em0BQ5zTjfVPlPvyAJKtqPrX1zluE4yGDzwQRWA8B63fMasG99ch7X1rXn5rhteXeNr79EwlOmbOQZ0ERkKYCqAcwGcCOAyETkxZbXtAK4G8JLfCQyDp5dsRaVp/lAvd+gLp36CG99MP56EqyIai2W3vF0xYFljWyea2tK3Vf/h/Qtw+r0fufhWZ5mczyG7BrLy21fKg05CXpnPQz870iSnc+zpjfcVr/T7XuOrsg2gr66sx+Q31iFm0+PaDbu9DkPG3U0O/VQA1apao6rdAGYBGGdeQVVrVXUtgMIaOs6l29/dgPbO2IDlTuWfgv65UT/09Mbx3rqdA07s5Lt/uetDnHrXh2m3sW3PAdTvPej6O5Pb9ntfoqC9a+B5kYnJb6zDqElzfNlWpmqaO3Cg25/9ydSvZ63Gd26fN2B58uaResWpKkprW/quh5Xb9uItFwOHhS1n7Rc3Af1YAHWm9/XGMs9EZIKIlIlIWXNz4U7+YFWWunr7XvzbPR+hvbPH88ny4cYmvLs2fauBZA7lmU9q8YsXV/VNopEPa41R4m70aeQ6tX3jrCvWG7pccV3LAdz93kZXu2LXszQXg0h59aMHF+Hnz5b2W5ZtT1ivStZZn9d2ud/ZZXX46ePLMGfdTgDAxdOW4jcW58emXW2obvLe1HJ3RxcunrbU8akXCEfX/2H5/DJVnQ5gOgAUFRUV5C3SLlZf+fQKtHfFUG7McgMMLJqJ9caxcWc7Thr5uX4n6PyNjZi/sRE/OfnLrtOxx+KxNNfaOgfm0DPN6WR67n9ctdvT0K258OrK/uPP//LFVVjnMIRDCK51V5bXtAAA3s9BhiEXmeIao39IXUv6J86xj3yc0fZnrdiOldv24rlltbjhnH/C1c+sQMWONst1w5Drd5NDbwBwnOn9SGPZoJc8fvG49j16pzumD3xQhf/42xJs3Gl9QqST+lTg5dzZ0txhWS5pp62zB+8ZOZ7+X+r+O9O57vmV+ONra41N5uciuGz6cvzr3emLouyoAq+trLcdc76nwCapcNLc3oVf2lTwZnJfcnsz293hPCTvzn2deLvcFH7yHEMXVjajud390MG1e5x7gfvJTUAvBTBGREaLyHAA4wEU5zZZ4ZZ6gqrN637rmAbicnPiOvEy082PH1yEsx9Z7Hr9384qxy9eXIXa3bmZqbyysb2vY0pPr2KbhxnRM83pLqvZg537OqHqvrdhcq1FVc34w6trcN/7lTZpSp+oRDGcl9QGqzugG9QZ9y/s976nN46zHlqEjzY19jvuv5+9Jr8Jg9umxcE/hjkGdFWNAZgIYC6AjQBmq2qFiEwRkQsAQES+JyL1AH4K4AkRGdj8IoLsg7fD3/nQMeW2dzakbDP9Rpvbu1w/GSQrTJNNItPxI079MOVCztSNbzpXLI6eXIL/fX6lp+0mK8QzrQDt7DkUIAdkBhx+wPve35R+hRQXT1uKR+ZXefqbsOhI+X33dHRjc1MHJr+xzlV/i14PLVfcBt8wBGkvXLVDV9USVf26qp6gqncay25V1WLjdamqjlTVI1T1C6r6zVwmOmzMh9wcWFNPQvO5EcR5cu5f3JUjhvkcTk1baW1L32u37f0/8HnIXqefy5xmrzfzxxZuQcm6nVhQ2eRqqOWV2/bikfnZTe9mlznY3dGVt+KldBkUq09EgP0Bt9AJg7xWikaN1UmnONRiwa7MNRdP3znZpqtmGzn4YsOB7hhiccXff+Yw23U27WzD90Z93vbznt44PtyY2zHX3dwAs7lJmsuza+8539XftHf24Mg0v1s6dse96I75GW3PD0FnMmyfxk2vw5APYtd/B0u3DOyWnPoYltoMb78xk0uys4RbbZ09nipc/NbZ09uvfD/ocUZOv3cBTr6t/yzpAyqHAezp6EKlxeh/Lfu7ccnjy3DdC849ONPJeuzzNJ95/Y1nLqvFfhdFPyVWldoFxHxTUe3/G/bGFTOWbM1q+9VNHVn9fVgxoDu4/MlPfd+m3ePk9+9dgO/dmVkuyI9Kt6ufWYGiO+Yf6sThYpt+dayx0uKyZc6ZDy3CORaVvv/11KdYY2pGavbUxzWo2OHPLOxec4+rt+9N+3m6suBb367wXAdQqNL9rlPeTdQhtXrMNGWqUCq1GdCzYHWMzTmu/mXrh3L2Vz9Tivkbmwb8bbqemGW16YOAH5JtkO2uo1Cc0xaViqlPQuc8nAjum9N0JLljzkac/+gS03bUU6WaF6lPdO+sMeWeLb7y2aW1abdnN+u8X3IVvLw8jbhd85WyRJ/HXM1v6nSz7j8UQk6S4AkDukeqmv4R2nQmOl2YXjh29AhFtM1cu0WnJSupv73V005y3B0vgemJxTU44cYS65uq40WdfgWv1/nO1vSdZBSJno/mOppRk+b4OruQn9y0UFlU2b/nuPm4amIjacULJQudYwzo2XA4h1IHAMrlDfx3s8t929ahgZB826Sjk1LKyt1Kbb6ZqdmliZxevzoEn/bfnHO7q2Sjp/Wt9MYVYx/5GDe9ub7f8ns9NnHc3NiOO+cM/P2CqDt5dWXdgGX9cr8urp6gM8hBD2UNsJVLRga2JQ4+d/DhpoFFOF40msaqsOvaDCTKf7/55c9l9V2ZWN+wD585bKinv3FzVDp7evGZw4amXdePC3VzY6ISrqH1IDbsPFR2n82Zs3JbdsVwVzy9ArtcjFFiJ92v0tMbx2FDE/lF8w2iqrEdZz88sL4jnpL5URz6bZrbuzinrkvMoWcg3cnT2BZcK5VsuHlk7eiK4aLHluLu95xzmdk4/9GPMbusf47tJ39dgjMfWuR73nF7S/+u2bno2i4QbDa1qkjWVdhvz10q3GYkfvb4Mtz69voBy3t9nJEplV0v42Vb3JX/qwIHTOXizzkUX2b7VOGl57XZ9pYDmLM2PC2KGNA9UqBvVprUnm0AHMc+DyurXKjdRZLJWDReVOxow/8zxnpJ9ZFFZbKV2aV1eXtyyib/no80rqhtwcxl2yy+O3ffmW4IDCv7LSo1py+u6Xv914+q/UiWrYunLU37ud21cP6jS3D9S4lmsawULXD7Dvagoytmm9MxU4TjgNvxkrYgywoPuGzNMM9lZ6LUQ5d88u/oiqFubyL3Hounz72tqU/f/NFr7tH11IWetup+C37E+Wzm202s537MHT/YHcMwlIt7wYCehVhvHN/601zcVuw8dE1PbxwLK8M1BnxN86FiAKvWHXbX07KUZnN1LfkbUa5lv7siLVV3wWNtfSueW1rbFzzOfGgRAOCSaUtR05woNnjbwwz3VtINxpZNzNqW5Uh+dq00LXtAe0yo3U3MbcvQzY0dnsZRmbpgS8YV63bOfnhRzgaoyxUGdI/Ms6EkW7G8ljI+tpV5Po8f4sRNT7gfPbio7/W7FuWAbi/h79+3AJ0pA3kt2NSEpdX+T/67wPVNUV0FyxteW4s/WdyQN1n0PM3UtIU1tp89tWSrqxl28sntU1Aq8/FWTcx5+7rp2lB1f2N44INKDA04OlU1dvS1c393zU5PM3apKqoa/TuH3GJA92hW6aHKumQOwk2uI9OMWMv+7oymJrv0iWWWy6fYNPMbYpEZ8pIrSx206ZpnS3H5U/73snXLquNWUJx+xwEz7OTpKd8uXVbNH92cCskiquT6t7xdgd+/uiajYgtVYEieyyhTMyVmDa0H8XsPTYPfLt9h2Zon1xjQQ87NCHtWumNxtOzvHpBLnvGJ9RgY2ZYVOt3Uwj4JhDn5fpfdet1avspt7dLldsiFVOZDbC5ySdZBVDd1eCpi8iug73M5PMDEl9KP+dPQ6r6Jp9U0ePnAduhZSOZq3XQZz+TUPNAdwztrMi+/vfzJ5a6LDixz6F6+TIE1da3Y1daJc775j/0++nhzM654eoWXrfXjZQKMzL/jUO5y9OQSX7ddm4f0W6lqbO+rB7CS2vY7ySroWq2ZOu2euemreRvJr/nNK+X4WdFI2/SYLapq7mtNZseqlZmV8U8ud7Xe/I1NWQ3aFYa28AzoWfASpDM51re8VYHXVzmXz1tp74q5CuYVO/bhm1/+HIZYRHQvJ+jFjy/tuxhSh3jNJpgD/k2AEZTV21s9re9XSYPTI7+bIWEBYEFlE655pnTAeqnDNTg92Qj8DXq/mVVu+5k5LV6a2aabUMRp/4IenRQYhEUufkz/luSlFt5pfI5UXbHefmWSuZJsR7vVpjbfPIFEOn4MR/roh9lNzBAVOzyeK5mYs3Zn30xMqVID12MLnNuAv7ayvl+xm1XsE/F3yKEtzf4PgZtNMU8YcuiDKqCvqWtF0R3z+9W8Z8PLsb/7PffjbNS1HMA3bn4fK7a6C6ZJmdaqz1m707KlzsXTluKnj1tXrro118Ps8Q/NK8yp0/yWbTNJAPjj6+k7uF2fprw4tZXLWps22uaixj+8uqZ/kYtF6K5q7MA9Hq6DbHz/vgUZ/V26a9rpiTcE8bzwAnp3LJ5xmWpyEoTUdtSZ+ijL8VPsZHoyLqwMT8uOJK9jd187syxHKSkMXm/iuZA6HEJXzLpCO7UJ6Z9NLaj2d+VmOFu3kvPiepVpBn3XvszHxPFTwQX0m95chx/ev9BTm9CkZDmxX0Nthq2jUL6beeVCvtvrh83PbJqbFpqrn8mu3sSNXPQkzXRS6F7Nb89WOwUX0JcYzfDW1LV67sWVPFTJ331zY7vtQWjZ311w01QtyUFHnkxkcrOlaLHL1fupNsueslYOZDgDl6piS5oWRWZ3lWzMWfAvuICezIVeOWMFznhgoeu/29F6ENVGJYqqYsGmJpz18GK8adNL7+yHF/V1Ay8UYXli+Paf/e2CTZQv7nsi99foYRji6YtrcnbDK7hmi0OtGkybdHTF8NnDB+7W/73no77XcT1UgWhV0bHvQA92d2TWuYKIBp+Lp3krKiuva8Vpx3/B93QUYA69/3tz54g3VtXjW3+ai08cih7iqn0jJIokil7M27GacJiIyC/pOnxlw1UOXUTGAvgLgKEAnlLVe1I+PxzATADfBbAHwKWqWutvUhNSO8Acf2OiV9+Z/+eLfeN3FJfvQGNbJ04eeRT+9tFm3PWfJ/X7m3fX7uyb/eaJRTV4YlENzv3WP6KjK4bpVxRlNYsLEZGTXLVfEMfeXSJDAVQBOAtAPYBSAJep6gbTOr8EcLKqXici4wFcpKqXpttuUVGRlpV5b6J21kOL+s3+QkRUaG445xu4/t+/ltHfishKVS2y+sxNkcupAKpVtUZVuwHMAjAuZZ1xAJ4zXr8G4MeSafsfB1FomkdEg9v9cytzsl03Af1YAOYJHuuNZZbrqGoMwD4AA0r8RWSCiJSJSFlzc2a1yUWjjs7o74iIwuKXZ5yQk+3mtZWLqk4HMB1IFLlkso07LzoJd150kvOKedYV68Xhw9zPSh/rjWPoEOnXkaE7FsfwYUMy2p5ZT28cAmBY0DMEeNAdi6OnN44jjBZKvXGFINGmffiwIX3LgURF+EPzqnDhd76MpVv24IrTvuq6Q0jyN47H1XJAslxSVcTiisMK6LhYUdWMO+CYxXrjfeeo1fE40B3D3w3vH6LW1rdiTV0rrvjXUVl9t1/HX/saV4Sj5MBNQG8AcJzp/UhjmdU69SIyDMDnkKgcHTS8Bl+rYJsM5plsz6wQA8bwYUP67X+yeerRRwwfsO6QIYI/nPMNAMDXvnik5+9JbiPfRASHDQ3HhZ8Nv4KX+RqwOh6pwRwATh55FE4eeVTW3+3X8Q9LIE9yc+WXAhgjIqNFZDiA8QCKU9YpBnCV8foSAB9pGPrBEhENIo45dFWNichEAHORaLY4Q1UrRGQKgDJVLQbwNIDnRaQaQAsSQZ+IiPLIVRm6qpYAKElZdqvpdSeAn/qbNCIi8qLwCluJiMgSAzoRUUQwoBMRRQQDOhFRRDCgExFFhOPgXDn7YpFmANsy/PNjAIRjep7cGgz7yX2MhsGwj0A49vOrqjrC6oPAAno2RKTMbrSxKBkM+8l9jIbBsI9A+PeTRS5ERBHBgE5EFBGFGtCnB52APBkM+8l9jIbBsI9AyPezIMvQiYhooELNoRMRUQoGdCKiiCi4gC4iY0WkUkSqRWRS0OnxSkRqRWSdiJSLSJmx7PMiMk9ENhv/H20sFxF51NjXtSJyimk7VxnrbxaRq+y+Lx9EZIaINInIetMy3/ZJRL5r/GbVxt/mfVYBm328TUQajGNZLiLnmT6bbKS3UkTOMS23PH+N+QY+NZa/Ysw9kFcicpyILBCRDSJSISK/NpZH7Vja7WfhH09VLZh/SIzHvgXA8QCGA1gD4MSg0+VxH2oBHJOy7D4Ak4zXkwDca7w+D8B7AATAaQA+NZZ/HkCN8f/RxuujA9ynHwA4BcD6XOwTgBXGumL87bkh2cfbAPzBYt0TjXPzcACjjXN2aLrzF8BsAOON148D+EUA+/glAKcYr48EUGXsS9SOpd1+FvzxLLQc+qkAqlW1RlW7AcwCMC7gNPlhHIDnjNfPAbjQtHymJiwHcJSIfAnAOQDmqWqLqu4FMA/A2DynuY+qLkZiYhMzX/bJ+OzvVXW5Jq6OmaZt5Y3NPtoZB2CWqnap6lYA1Uicu5bnr5FL/RGA14y/N/9eeaOqO1V1lfG6HcBGJCaAj9qxtNtPOwVzPAstoB8LoM70vh7pD0QYKYAPRGSliEwwlv2Dqu40Xu8C8A/Ga7v9LYTfwa99OtZ4nbo8LCYaxQ0zkkUR8L6PXwDQqqqxlOWBEZFRAL4D4FNE+Fim7CdQ4Mez0AJ6FJyuqqcAOBfA9SLyA/OHRs4lUm1Jo7hPhmkATgDwzwB2Angw0NT4REQ+C+B1AL9R1TbzZ1E6lhb7WfDHs9ACegOA40zvRxrLCoaqNhj/NwF4E4nHtkbjcRTG/03G6nb7Wwi/g1/71GC8Tl0eOFVtVNVeVY0DeBKJYwl438c9SBRXDEtZnncichgSQe5FVX3DWBy5Y2m1n1E4noUW0EsBjDFqkIcjMRl1ccBpck1EjhCRI5OvAZwNYD0S+5BsCXAVgLeN18UArjRaE5wGYJ/x6DsXwNkicrTxWHi2sSxMfNkn47M2ETnNKJu80rStQCWDnOEiJI4lkNjH8SJyuIiMBjAGicpAy/PXyPUuAHCJ8ffm3ytvjN/3aQAbVfUh00eROpZ2+xmJ45mPmlc//yFRs16FRO3yTUGnx2Paj0eiJnwNgIpk+pEoc/sQwGYA8wF83lguAKYa+7oOQJFpWz9HonKmGsA1Ae/Xy0g8ovYgUV74337uE4AiJC6uLQD+BqOHcwj28XljH9YicdF/ybT+TUZ6K2FqyWF3/hrnxgpj318FcHgA+3g6EsUpawGUG//Oi+CxtNvPgj+e7PpPRBQRhVbkQkRENhjQiYgiggGdiCgiGNCJiCKCAZ2IKCIY0ImIIoIBnYgoIv4/Lih+0KeNONEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t_H.cpu().numpy()[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(0.3829)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_W.mm(t_H).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(1.)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_vectors.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyTorch Autograd Introduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[1., 1.],\n",
      "        [1., 1.]], requires_grad=True)\n"
     ]
    }
   ],
   "source": [
    "x = Variable(torch.ones(2, 2), requires_grad=True)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[1., 1.],\n",
      "        [1., 1.]])\n"
     ]
    }
   ],
   "source": [
    "print(x.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "None\n"
     ]
    }
   ],
   "source": [
    "print(x.grad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[3., 3.],\n",
      "        [3., 3.]], grad_fn=<AddBackward0>)\n"
     ]
    }
   ],
   "source": [
    "y = x + 2\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[27., 27.],\n",
      "        [27., 27.]], grad_fn=<MulBackward0>) tensor(108., grad_fn=<SumBackward0>)\n"
     ]
    }
   ],
   "source": [
    "z = y * y * 3\n",
    "out = z.sum()\n",
    "print(z, out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[18., 18.],\n",
      "        [18., 18.]])\n"
     ]
    }
   ],
   "source": [
    "out.backward()\n",
    "print(x.grad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using Autograd for NMF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "lam=1e6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "pW = Variable(tc.FloatTensor(m,d), requires_grad=True)\n",
    "pH = Variable(tc.FloatTensor(d,n), requires_grad=True)\n",
    "pW.data.normal_(std=0.01).abs_()\n",
    "pH.data.normal_(std=0.01).abs_();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "def report():\n",
    "    W,H = pW.data, pH.data\n",
    "    print((M-pW.mm(pH)).norm(2).data[0], W.min(), H.min(), (W<0).sum(), (H<0).sum())\n",
    "\n",
    "def penalty(A):\n",
    "    return torch.pow((A<0).type(tc.FloatTensor)*torch.clamp(A, max=0.), 2)\n",
    "\n",
    "def penalize(): return penalty(pW).mean() + penalty(pH).mean()\n",
    "\n",
    "def loss(): return (M-pW.mm(pH)).norm(2) + penalize()*lam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = Variable(t_vectors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = torch.optim.Adam([pW,pH], lr=1e-3, betas=(0.9,0.9))\n",
    "lr = 0.05\n",
    "#report()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = pH.data.cpu().numpy()\n",
    "show_topics(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(h[0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing Appraoches\n",
    "#### Scikit-Learns's NMF\n",
    "\n",
    "    1.  Fast\n",
    "    2.  No Parameter tuning\n",
    "    3.  Relies on decaded of academic research\n",
    "   \n",
    "#### Using PyTorch and SGD\n",
    "    1. Took an hour to implement, didn't have to be NMF experts\n",
    "    2. Parameters were dfiddy\n",
    "    3. Not as fast\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.0 64-bit",
   "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.0"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "c0f2ef85a013b3e049ebe5ad4c39b42bd34e85c4b0686b7bef6ce5ff62aa91e1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}