奇异值分解 SVD

奇异值分解（Singular Value Decomposition）是矩阵分解的一种，是一个应用非常广泛的无监督算法，常用于数据降维、压缩。

在讲奇异值分解，首先讲讲奇异值分解的特例特征值分解

特征值分解

给定实对称方阵 A\in R^n，A 一定可以分解为如下形式：

大学线代学过对于实对称方阵 A 一定可以进行特征值分解表示为A=P^{-1}\Sigma P,P 为 A 的特征向量组成的矩阵，\Sigma为 A 的特征值组成的对角矩阵,\lambda_i对应的特征向量为q_i，矩阵的逆运算不好求对 P 正交化单位化得到正交阵 Q，则可以分解为A=Q^T\Sigma Q。

正交矩阵：若AA^T=E，A 为正交矩阵。正交阵的各行各列是单位向量且两两正交，A^T=A^{-1}。

正交：q_i^Tq_j=0，即q_i,q_j内积为 0 则q_i,q_j正交。

奇异值分解

特征值分解要求 A 是方阵，而在实际中 A\in R^{m*n}，m 为样本个数，n 为特征个数，m 通常都是远大于 n 的，这时候需要一般化的矩阵分解即奇异值分解：

其中 U，V 均为正交阵，U 中的向量称为左奇异向量，V 中的向量称为右奇异向量，\Sigma仅在主对角线有值，这些值称为奇异值，其他元素均为 0。

奇异值求解

AA^T\in R^{m*m}，A^TA\in R^{n*n}为对称半正定矩阵，证明：

对称性：(AA^T)^T=(A^T)^TA^T=AA^T

半正定性：对于任意x\in R^m,x^TAA^Tx=(A^Tx)^T(A^Tx)\ge0

因为对称矩阵一定可以对角化，所以我们可以通过对角化AA^T,A^TA来间接求奇异值：

AA^T=U\Sigma V^TV\Sigma^TU^T=U\Sigma \Sigma^TU^T \\ A^TA=V\Sigma^TU^TU\Sigma V^T=V\Sigma^T\Sigma V^T

$\Sigma \Sigma^T=\begin{bmatrix}
\sigma_1^2 &0 &0 &0 \
0 &\sigma_2^2 &0 &0 \
\vdots &\vdots &\ddots &\vdots \
\cdots &\cdots &\cdots &\sigma_m^2
\end{bmatrix}{m*m}
$ $\Sigma^T \Sigma=\begin{bmatrix}
\sigma_1^2 &0 &0 &0 \
0 &\sigma_2^2 &0 &0 \
\vdots &\vdots &\ddots &\vdots \
\cdots &\cdots &\cdots &\sigma_n^2
\end{bmatrix}{n*n}$

可以看出AA^T,A^TA的特征值对角矩阵是 A 的奇异值矩阵的平方，即AA^T,A^TA特征值是 A 的奇异值的平方。

奇异值分解分析

由上面可知对于任意A\in R^{m*n}（n<m）可以分解为 n 个矩阵之和:

A=\sigma_1u_1v_1^T +\sigma_2u_2v_2^T + \cdots + \sigma_nu_nv_n^T

\Sigma矩阵中的奇异值从大到小排列，实际中奇异值从大到小排列的衰减速度非常快，那么我们可以选取前 r 大个奇异值和奇异向量外积相加来近似 A，即 A 的特征值分解可表示为：

A \approx \sigma_1u_1v_1^T +\sigma_2u_2v_2^T + \cdots + \sigma_ru_rv_r^T=U_{m*r}\Sigma_{r*r}V_{r*n}^T

内存角度

分解之前 A 需要存储 m*n 个数，分解之后需要存储 r(m+1+n)，而 r 是远小于 min(m,n)的。

运算次数角度

假定A\in R^{100000*1000},x\in^{1000*1},考虑计算 Ax 需要计算的乘法次数：

分解之前需要计算 100000*1000={10}^8次,分解后通过矩阵乘法的结合律优化计算次数，假定 r 取 100,从右往左依次相乘需要计算的乘法次数为{}10^5+10^4+10^6,可以看出分解后乘法次数大约为之前的百分之一。

SVD 在图像压缩中的应用

输入图像 960x960x3，选择前 r 大个奇异值重构图像，这里选择 r 为 200，100，50，10，下面我们看看图片展示：

可以看出，当 r 取 50 已经能保持不错的精度了。下面看看奇异值从大到小排列的的衰减变化以及前 n 个奇异值占所有奇异值的比例：

可以看出奇异值呈指数衰减，前 50 个奇异值就占了奇异值总和的 80%，前 200 个奇异值就占了奇异值总和的 90%。

代码

import cv2
import matplotlib.pyplot as plt
from PIL import Image
import numpy as np

#读取图片
img =Image.open(r'D:\pictures\lisa.jpg')
img_array = np.asarray(img).reshape(960, -1)
#奇异值分解
U,sigma,V = np.linalg.svd(img_array)
#选取前sval_nums个奇异值重构图像
#sval_nums为保留最大奇异值个数
sval_nums = 50
img_restruct1 = (U[:,0:sval_nums]).dot(np.diag(sigma[0:sval_nums])).dot(V[0:sval_nums,:])
img_restruct1 = img_restruct1.reshape(960,960,3)
sval_nums = 10
img_restruct2 = (U[:,0:sval_nums]).dot(np.diag(sigma[0:sval_nums])).dot(V[0:sval_nums,:])
img_restruct2 = img_restruct2.reshape(960,960,3)
#图示
fig, ax = plt.subplots(1, 3, figsize=(24,32))
ax[0].imshow(img)
ax[0].set(title='src')
ax[1].imshow(img_restruct1.astype(np.uint8))
ax[1].set(title='nums of sigma = 50')
ax[2].imshow(img_restruct2.astype(np.uint8))
ax[2].set(title='nums of sigma = 10')
#观察奇异值变化
x = np.arange(960)
sum_ = sum(sigma)
tmp = sigma.copy()
for i in range(len(sigma)):
    tmp[i] = sum(sigma[:i]) / sum_
fig,ax = plt.subplots(1,2,figsize=(20,10))
ax[0].plot(x, sigma)
ax[0].set_xlabel('index')
ax[0].set_ylabel('singular value')
ax[1].plot(x, tmp)
ax[1].set_xlabel('index')
ax[1].set_ylabel('rate')

总结

特征值分解和奇异值分解都是常见的矩阵分解算法，奇异值分解更为一般，通过奇异值分解可以进行数据降维提取矩阵的重要信息。奇异值从大到小排列通常呈指数衰减，所以我们选取前 r 个奇异值（r 远小于矩阵维度）及对应的左右奇异向量外积相加基本就能近似原矩阵了。

Reference

SVD（奇异值分解）小结 - EndlessCoding - 博客园