《线性代数及其应用》笔记18.矩阵分解-奇异值分解.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 奇异值分解\n",
    "分解条件：条件最低，任何矩阵都有可能。\n",
    "## 1.正交对角化\n",
    "$\\quad$**对称矩阵**是一个满足$A^{T}=A$的矩阵$A$。自然地，对称矩阵是一个方阵。  \n",
    "$\\quad$**定理**：如果$A$是对称矩阵，那么不同特征空间的任意两个特征向量是正交的。（补充：无论是不是对称矩阵，相异的特征值对应的特征向量组成的集合都是线性无关的）证：\n",
    "\n",
    "<img src=\"./_image/18_1.png\" width=\"400\" height=\"350\" />  \n",
    "\n",
    "$\\quad$我们知道，[对角化](《线性代数及其应用》笔记11.矩阵相似与对角化.ipynb)是对于一个方阵$A$，能找到$n$个线性无关的特征向量组成的矩阵$P$、特征值组成的对角矩阵$D$使得$A=PDP^{-1}$。而当$P$为正交矩阵，则称上式为**正交对角化**。我们再引出**定理**：一个$n\\times n$矩阵$A$可正交对角化的充要条件是$A$是对称矩阵。  \n",
    "$\\quad$证明如下：\n",
    "\n",
    "<img src=\"./_image/18_2.png\" width=\"300\" height=\"250\" />  \n",
    "\n",
    "$\\quad$给一个例子：\n",
    "\n",
    "<img src=\"./_image/18_3.png\" width=\"400\" height=\"350\" />  \n",
    "\n",
    "## 2.奇异值分解介绍\n",
    "$\\quad$对角化在许多应用中都很重要，但是要对一个矩阵进行对角化，该矩阵必须是方阵。下面将介绍另一种矩阵分解：**奇异值分解**，该分解对任意$m\\times n$矩阵$A$都有可能，是线性代数应用中最有用的矩阵分解，被用于求矩阵的秩、降维、压缩图像等。  \n",
    "$\\quad$奇异值分解的形式如下：\n",
    "\n",
    "$$A=U\\varSigma V^{T}$$\n",
    "\n",
    "其中$U$为$m\\times m$正交矩阵，$V$为$n\\times n$正交矩阵，$U$的列称为$A$的左奇异向量，$V$的列称为$A$的右奇异向量；$\\varSigma$为$m\\times n$”对角矩阵“，$D$为一个$r\\times r$对角矩阵，$\\sigma_{1},\\cdots,\\sigma_{r}$为矩阵的$r$个奇异值，形式如下：\n",
    "\n",
    "<img src=\"./_image/18_4.png\" width=\"300\" height=\"250\" />  \n",
    "\n",
    "## 3.奇异值分解算法\n",
    "$\\quad$令矩阵$A$为$m\\times n$矩阵，我们希望得到$A$的奇异值分解$A=U\\varSigma V^{T}$，接下来会从这个式子出发，逐步推出$U,\\varSigma,V$的求解方法以及之间的关系。  \n",
    "### 3.1结论一\n",
    "$\\quad$我们计算$A^{T}A,AA^{T}$有：\n",
    "\n",
    "<img src=\"./_image/18_5.png\" width=\"300\" height=\"250\" />  \n",
    "\n",
    "形式与对角化是一致，其中而对称矩阵是可以正交对角化的，为此可以下一个结论：$A^{T}A$以及$AA^{T}$的特征向量规范化（指化为单位正交基）后可构成$V$和$U$，特征值的非零项开方称之为**奇异值**，记作$\\sigma_{i}$，奇异值构成$\\varSigma$。\n",
    "### 3.2结论二\n",
    "$\\quad$令$\\{\\mathbf{v}_{1},\\cdots,\\mathbf{v}_{n}\\}$是$\\Bbb{R}^n$的单位正交基且构成$A^{T}A$的特征向量，$\\lambda_{1},\\cdots,\\lambda_{n}$是$A^{T}A$对应的特征值，那么我们有:\n",
    "\n",
    "<img src=\"./_image/18_6.png\" width=\"300\" height=\"250\" />  \n",
    "\n",
    "这里得到两个结论：  \n",
    "* 奇异值非负；\n",
    "* $\\{\\mathbf{v}_{1},\\cdots,\\mathbf{v}_{n}\\}$为$A^{T}A$的特征向量，$A$的奇异值是向量$A\\mathbf{v}_{1},\\cdots,A\\mathbf{v}_{n}$的长度. \n",
    "\n",
    "### 3.3结论三\n",
    "$\\quad$我们对奇异值分解的式子变换如下：\n",
    "$$AV=U\\varSigma$$\n",
    "我们将其展开：\n",
    "\n",
    "<img src=\"./_image/18_7.png\" width=\"400\" height=\"350\" />  \n",
    "\n",
    "得出结论：\n",
    "\n",
    "<img src=\"./_image/18_8.png\" width=\"100\" height=\"50\" />  \n",
    "\n",
    "是奇异值、原矩阵、左右奇异向量的关系，在我们通过结论一求出了奇异值、其中一个奇异向量时，可以直接求出另一个奇异向量。通过这些知识我们已经可以对A进行奇异值分解了。  \n",
    "### 3.4实战\n",
    "$\\quad$例子：\n",
    "\n",
    "<img src=\"./_image/18_9.png\" width=\"400\" height=\"350\" />  \n",
    "\n",
    "## 4.奇异值分解应用\n",
    "### 4.1矩阵求秩\n",
    "$\\quad$计算机估计大矩阵$A$的[秩](《线性代数及其应用》笔记07.维数和秩.ipynb)时，最可靠的方法就是计算非零奇异值的个数（特别小的非零奇异值在实际计算中假定为0），可参考MATLAB的rank源码。  \n",
    "$\\quad$**定理一**：假若$\\{\\mathbf{v}_{1},\\cdots,\\mathbf{v}_{n}\\}$是包含$A^{T}A$的特征向量的$\\Bbb{R}^{n}$上的单位正交基，重新整理使得对应的$A^{T}A$的特征值满足$\\lambda_{1}\\geq\\cdots\\geq\\lambda_{n}$。假如$A$有$r$个非零奇异值，那么$\\{A\\mathbf{v}_{1},\\cdots,A\\mathbf{v}_{r}\\}$是$ColA$的一个正交基，且$rankA=r$。  \n",
    "$\\quad$定理一说明了，矩阵的奇异值的个数即为矩阵的秩。我们结合结论三得出定理二：  \n",
    "$\\quad$**定理二**：给定$m\\times n$矩阵$A$的一个奇异值分解，取$\\mathbf{u}_{1},\\cdots,\\mathbf{u}_{m}$是左奇异向量，$\\mathbf{v}_{1},\\cdots,\\mathbf{v}_{n}$是右奇异向量，且$\\sigma_{1},\\cdots,\\sigma_{n}$是奇异值，$r$为$A$的秩；$\\{\\mathbf{u}_{1},\\cdots,\\mathbf{u}_{r}\\}$是$ColA$的一个单位正交基，$\\{\\mathbf{v}_{r+1},\\cdots,\\mathbf{v}_{n}\\}$是$NulA$的一个单位正交基，$\\{\\mathbf{v}_{1},\\cdots,\\mathbf{v}_{r}\\}$是$RowA$的一个单位正交基。  \n",
    "$\\quad$上述的两条定理表明了奇异值个数、列空间、行空间、零空间之间的关系。举一个例子：\n",
    "\n",
    "<img src=\"./_image/18_10.png\" width=\"400\" height=\"350\" />  \n",
    "\n",
    "### 4.2图像压缩\n",
    "$\\quad$我们知道一个图像其实就是一个矩阵，矩阵上每个元素都代表着像素。  \n",
    "$\\quad$假设$A$为一个$m\\times n$矩阵，对$A$做奇异值分解有：$A=U\\varSigma V^{T}$，我们展开计算有：\n",
    "\n",
    "<img src=\"./_image/18_11.png\" width=\"400\" height=\"350\" />  \n",
    "\n",
    "其中$\\sigma_{i}$按照从大到小排列，$\\mathbf{u}_{i}\\mathbf{v}_{i}^T$为一个秩一矩阵。我们知道，$r$个不同的秩一矩阵相加得到的矩阵秩为$r$，可得到两个结论：  \n",
    "1. $A$的奇异值个数为$A$的秩；\n",
    "2. 每个奇异值可看作是一个秩一矩阵的权重。  \n",
    "\n",
    "$\\quad$奇异值一般数值间差异很大，我们可以取前$k$项来近似等于$A$从而压缩：\n",
    "\n",
    "<img src=\"./_image/18_12.png\" width=\"200\" height=\"150\" />  \n",
    "\n",
    "我们读取图片，再对矩阵做奇异值分解，代码如下:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np  \n",
    "import cv2  \n",
    "import matplotlib.pyplot as plt  \n",
    "from scipy.linalg import svd  \n",
    "  \n",
    "def recover_from_k(U, S, V, k):  \n",
    "    return np.dot(np.dot(U[:, :k], np.diag(S[:k])), V[:k, :])  \n",
    "  \n",
    "if __name__ == \"__main__\":  \n",
    "  \n",
    "    pic = cv2.imread('/Users/lin/study/linear/test.jpg', cv2.IMREAD_GRAYSCALE)  \n",
    "    plt.imshow(pic, cmap=plt.cm.gray)  \n",
    "    plt.show()  \n",
    "    U, S, V = svd(pic)  \n",
    "    print(sorted(S, reverse=True))  \n",
    "  \n",
    "    pic = recover_from_k(U, S, V, 10)  \n",
    "    plt.imshow(pic, cmap=plt.cm.gray)  \n",
    "    plt.show()  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图片大小为250x250，共有250个奇异值，前5个奇异值为：\n",
    "\n",
    "$$28703.2,5937.4,4087.2,3820.3,3148.4$$\n",
    "\n",
    "$\\quad$原图为：\n",
    "\n",
    "<img src=\"./_image/18_13.png\" width=\"200\" height=\"150\" />  \n",
    "\n",
    "我们取前10个奇异值还原，如下:\n",
    "\n",
    "<img src=\"./_image/18_14.png\" width=\"200\" height=\"150\" />  \n",
    "\n",
    "取前20个还原，如下：\n",
    "\n",
    "<img src=\"./_image/18_15.png\" width=\"200\" height=\"150\" />  \n",
    "\n",
    "取前50个还原，可以看出已经非常接近原图了：\n",
    "\n",
    "<img src=\"./_image/18_16.png\" width=\"200\" height=\"150\" />  \n",
    "\n",
    "$\\quad$原图为250x250，总共有62500个元素。我们取前50个奇异值还原，$\\mathbf{u}_{1}$和$\\mathbf{v}_{1}^T$均包含250个元素，总共需要50*(250+250+1)=25050个元素，成功达到压缩目地。  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}