fastai_NLA_nbs/NMF_and_SVD.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "# 多行输出\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NMF and SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![NMF and SVD](../imgs/NMF1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 传统的奇异值分解方法只能对稠密矩阵进行分解，即不允许所分解矩阵有空值。因而，若采用奇异值分解，需要首先填充数值，这样造成了两个问题：\n",
    "    - 其一，填充大大增加了数据量，增加了算法复杂度\n",
    "    - 其二，简单粗暴的数据填充很容易造成数据失真\n",
    "\n",
    "- 使用矩阵分解具有以下优点：\n",
    "    - 比较容易编程实现，随机梯度下降方法依次迭代即可训练出模型\n",
    "    - 比较低的时间和空间复杂度，高维矩阵映射为两个低维矩阵节省了存储空间，训练过程比较费时，但是可以离线完成\n",
    "    - 评分预测一般在线计算，直接使用离线训练得到的参数，可以实时推荐\n",
    "    - 预测的精度比较高，预测准确率要高于基于领域的协同过滤以及内容过滤等方法\n",
    "    - 非常好的扩展性，很方便在用户特征向量和物品特征向量中添加其它因素\n",
    "\n",
    "- 矩阵分解的不足主要有：\n",
    "    - 模型训练比较费时。\n",
    "    - 推荐结果不具有很好的可解释性，分解出来的用户和物品矩阵的每个维度无法和现实生活中的概念来解释，无法用现实概念给每个维度命名，只能理解为潜在语义空间。\n",
    "\n",
    "- 非负矩阵分解(non-negative matrix factorization)1，或非负矩阵近似(non-negative matrix approximation)，是多变量分析和线性代数的算法。\n",
    "    - 输入：非负矩阵V\n",
    "    - 输出：两个非负矩阵W和H\n",
    "    - 目标：V=WH\n",
    "- 优点\n",
    "    - 处理大规模数据更快更便捷；\n",
    "    - 实现简便性、分解形式和分解结果上的可解释性，占用存储空间少。\n",
    "- 缺点\n",
    "    - NMF中只用一层表示隐变量，无法处理复杂学习问题；\n",
    "    - NMF只约束了W和H的非负性（这是唯一先验，只要求满足这个），而没有考虑到对于该先验，H内部元素间的相关性。\n",
    "- 应用领域\n",
    "　　- 计算机视觉，文档聚类，化学统计学，音频信号处理，推荐系统"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![NMF and SVD](../imgs/NMF2.png)\n",
    "![NMF and SVD](../imgs/NMF3.png)"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "## 数据集\n",
    "newsgroups数据集包含20个话题的18000篇新闻。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Downloading 20news dataset. This may take a few minutes.\n",
      "Downloading dataset from https://ndownloader.figshare.com/files/5975967 (14 MB)\n"
     ]
    }
   ],
   "source": [
    "categories = ['alt.atheism', 'talk.religion.misc', 'comp.graphics', 'sci.space']\n",
    "remove = ('headers', 'footers', 'quotes')\n",
    "# 只下载 上述的4个话题数据并删除 '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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2034,), (2034,))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 把文件名理解为 x 话题理解为 y\n",
    "newsgroups_train.filenames.shape, newsgroups_train.target.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "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"
     ]
    }
   ],
   "source": [
    "# 查看数据信息\n",
    "len(newsgroups_train.data[:3])\n",
    "print('\\n'.join(newsgroups_train.data[:3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['alt.atheism', 'comp.graphics', 'sci.space', 'talk.religion.misc']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array(['/home/lab/scikit_learn_data/20news_home/20news-bydate-train/comp.graphics/38816',\n",
       "       '/home/lab/scikit_learn_data/20news_home/20news-bydate-train/talk.religion.misc/83741',\n",
       "       '/home/lab/scikit_learn_data/20news_home/20news-bydate-train/sci.space/61092',\n",
       "       '/home/lab/scikit_learn_data/20news_home/20news-bydate-train/alt.atheism/51306',\n",
       "       '/home/lab/scikit_learn_data/20news_home/20news-bydate-train/sci.space/60903'],\n",
       "      dtype='<U90')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([1, 3, 2, 0, 2])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.target_names  # 数据集的类别\n",
    "newsgroups_train.filenames[:5] # 前三个样本的x和y\n",
    "newsgroups_train.target[:5] # 类别是类名的索引"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 生成词典"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        ...,\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0]], dtype=int64)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(2034, 26576)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectorizer = CountVectorizer(stop_words='english') # sparse matrix\n",
    "vectors = vectorizer.fit_transform(newsgroups_train.data).todense() # (documents, vocab)\n",
    "vectors  # 词汇统计\n",
    "vectors.shape #, vectors.nnz / vectors.shape[0], row_means.shape\n",
    "# 总共出现了 26576 个词汇，一共有2034个训练样本"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26576,)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array(['factors', 'factory', 'facts', 'factsnet', 'factual', 'factually',\n",
       "       'faculty', 'fade', 'fades', 'fading', 'fag', 'faget', 'faggots',\n",
       "       'fahrenheit', 'fai', 'fail', 'failed', 'failing', 'fails',\n",
       "       'failsafe', 'failure', 'failures', 'faily', 'fain', 'faint',\n",
       "       'fainter', 'fair', 'fairchild', 'fairfield', 'fairgrove',\n",
       "       'fairing', 'fairly', 'fairness', 'fairy', 'fait', 'faith',\n",
       "       'faithful', 'faithfully', 'faiths', 'fake', 'faking', 'falguiere',\n",
       "       'falkow', 'fall', 'fallacies', 'fallacious', 'fallaciously',\n",
       "       'fallacy', 'fallacys', 'fallback', 'fallen', 'fallible', 'falling',\n",
       "       'falls', 'false', 'falsehood', 'falsely', 'falseness',\n",
       "       'falsifiable', 'falsification', 'falsified', 'falsify',\n",
       "       'faltering', 'fama', 'fame', 'familar', 'familia', 'familial',\n",
       "       'familiar', 'familiarity', 'familiarize', 'familiarized',\n",
       "       'families', 'familiy', 'family', 'famine', 'famous', 'fan',\n",
       "       'fanatic', 'fanatical', 'fanatically', 'fanatics', 'fanciful',\n",
       "       'fancy', 'fanning', 'fans', 'fantasies', 'fantastical', 'fantasy',\n",
       "       'fantasyland', 'faq', 'faq1', 'faqs', 'far', 'faraway', 'farce',\n",
       "       'faries', 'farin', 'farm', 'farming'], dtype='<U80')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 查看词汇\n",
    "vocab = np.array(vectorizer.get_feature_names())\n",
    "vocab.shape\n",
    "vocab[10000:10100]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVD 解决方案"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SVD算法将一个矩阵分解为一个正交列矩阵和一个正交行矩阵\n",
    "![SVD](../imgs/svd_fb.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4min 30s, sys: 34.4 s, total: 5min 5s\n",
      "Wall time: 26.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "U, s, Vh = linalg.svd(vectors, full_matrices=False) # 这里full_matrices=False保证，我们得到的是 (M, K) 和 (K, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'U' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-12-92c6c110279b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# 上述运行太慢\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mU\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mVh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'U' is not defined"
     ]
    }
   ],
   "source": [
    "# 上述运行太慢\n",
    "print(U.shape, s.shape, Vh.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, s, Vh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 奇异值\n",
    "plt.plot(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(s[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_topics, num_top_words=8, 10  # 从图像可以看到，刚开始奇异值较大但是在后面的部分奇异值趋于平缓，\n",
    "# 所以，我们选取前8个作为最重要的特征"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_topics(a):\n",
    "#      排序获取索引，升序， num_top_words 个\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vh[5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_topics(Vh[:num_topics])\n",
    "newsgroups_train.target_names\n",
    "c = newsgroups_train.target[:10]\n",
    "[newsgroups_train.target_names[i] for i in c]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NMF 方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- NMF是限制元素是非负的，SVD是限制矩阵是正交的\n",
    "- 正值更容易解释\n",
    "- 非精确分解\n",
    "- $V=WH$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m, n = vectors.shape\n",
    "m, n\n",
    "d=5  # num topics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = decomposition.NMF(n_components=d, random_state=1)\n",
    "W1 = clf.fit_transform(vectors)\n",
    "H1 = clf.components_\n",
    "W1.shape, H1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_topics(H1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 使用 SGD 学习NMF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "损失函数看做 MSE 的回归问题"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 非负分解,当元素 < 0 时,惩罚\n",
    "def penalty(A):\n",
    "    return torch.pow((A<0) * torch.clamp(A, max=0.), 2)\n",
    "\n",
    "def penalize():\n",
    "    return penalty(pW).mean() + penalty(pH).mean()\n",
    "\n",
    "def loss():\n",
    "    return (M-pW.mm(pH)).norm(2) + penalize()*lam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* NMF 计算快的原因是,其只计算关心的那些列,(原来是26576列,现在只保留5列)而SVD要全部计算,可以通过截断SVD来实现.(选取奇异值最大的几列)\n",
    "- 计算完整的SVD非常缓慢\n",
    "- 随机SVD只计算我们需要的列.比如我们计算的是 (2034,5),(5,)(5,26576)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LU 分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在线性代数与数值分析中，LU分解是矩阵分解的一种，将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积，有时需要再乘上一个置换矩阵。LU分解可以被视为高斯消元法的矩阵形式。在数值计算上，LU分解经常被用来解线性方程组、且在求反矩阵和计算行列式中都是一个关键的步骤\n",
    "\n",
    "[LU分解](https://zh.wikipedia.org/wiki/LU%E5%88%86%E8%A7%A3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(linalg.lu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl, u = linalg.lu(A, permute_l=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.allclose(A - pl @ u, np.zeros((4, 4)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QR 分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QR分解法是三种将矩阵分解的方式之一。这种方式，把矩阵分解成一个半正交矩阵与一个上三角矩阵的积。QR分解经常用来解线性最小二乘法问题。QR分解也是特定特征值算法即QR算法的基础\n",
    "\n",
    "[QR分解](https://zh.wikipedia.org/wiki/QR%E5%88%86%E8%A7%A3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}