SVD奇异值分解

来源：互联网发布：mac应用商店英文编辑：程序博客网时间：2024/05/29 15:53

奇异值分解是一个有着很明显的物理意义的一种方法，（马同学一个形象生动的例子）它可以将一个比较复杂的矩阵用更小更简单的几个子矩阵的相乘来表示（矩阵相乘的本质），这些小矩阵描述的是矩阵的重要的特性。就像是描述一个人一样，给别人描述说这个人长得浓眉大眼，方脸，络腮胡，而且带个黑框的眼镜，这样寥寥的几个特征，就让别人脑海里面就有一个较为清楚的认识，实际上，人脸上的特征是有着无数种的，之所以能这么描述，是因为人天生就有着非常好的抽取重要特征的能力，让机器学会抽取重要的特征，SVD是一个重要的方法。

在机器学习领域，有相当多的应用与奇异值都可以扯上关系，比如做feature reduction的PCA，做数据压缩（以图像压缩为代表）的算法，还有做搜索引擎语义层次检索的LSI（Latent Semantic Indexing）。后面会给我出图像压缩的例子。

奇异值分解：来源

下面直观的从一张图片出发，让我们来看看奇异值代表什么意义。这是lena的一张照片，像素为高度512*宽度512,当然不一定就是标准的方型。

我们都知道，图片实际上对应着一个矩阵，矩阵的大小就是像素大小，比如这张图对应的矩阵阶数就是512*512，矩阵上每个元素的数值对应着像素值。我们记这个像素矩阵为 $A$ .

现在我们对矩阵 $A$ 进行奇异值分解。直观上，奇异值分解将矩阵分解成若干个秩一矩阵之和，用公式表示就是：
$(1) \quad\quad \qquad A = \sigma_1 u_1v_1^{\rm T}+\sigma_2 u_2v_2^{\rm T}+...+\sigma_r u_rv_r^{\rm T}$
其中等式右边每一项前的系数 $\sigma$ 就是奇异值， $u$ 和 $v$ 分别表示列向量，秩一矩阵的意思是矩阵秩为1。注意到每一项 $uv^{\rm T}$ 都是秩为1的矩阵。我们假定奇异值满足 $\sigma_1\geq\sigma_2\geq...\geq\sigma_r>0$ （奇异值大于0是个重要的性质，但这里先别在意），如果不满足的话重新排列顺序即可，这无非是编号顺序的问题。

既然奇异值有从大到小排列的顺序，我们自然要问，如果只保留大的奇异值，舍去较小的奇异值，这样(1)式里的等式自然不再成立，那会得到怎样的矩阵——也就是图像？

令 $A_1=\sigma_1 u_1v_1^{\rm T}$ ，这只保留(1)中等式右边第一项，然后作图：

结果就是完全看不清是啥……我们试着多增加几项进来： $A_5=\sigma_1 u_1v_1^{\rm T}+\sigma_2 u_2v_2^{\rm T}+\sigma_3 u_3v_3^{\rm T}+\sigma_4 u_4v_4^{\rm T}+\sigma_5 u_5v_5^{\rm T}$ ，再作图

隐约可以辨别……但还是很模糊，毕竟我们只取了5个奇异值而已。下面我们取20个奇异值试试，也就是(1)式等式右边取前20项构成 $A_{20}$

虽然还有些马赛克般的模糊，但我们总算能辨别出这是lena的脸。当我们取到(1)式等式右边前50项时：

我们得到和原图差别不大的图像。也就是说当 $k$ 从1不断增大时， $A_k$ 不断的逼近 $A$ 。让我们回到公式

$(1) \quad\quad \qquad A = \sigma_1 u_1v_1^{\rm T}+\sigma_2 u_2v_2^{\rm T}+...+\sigma_r u_rv_r^{\rm T}$

矩阵 $A$ 表示一个512*512的矩阵，需要保存512*512=262144个元素的值。等式右边 $u$ 和 $v$ 分别是512*1和512*1的向量，每一项有1+512+512=1025个元素。如果我们要存储很多高清的图片，而又受限于存储空间的限制，在尽可能保证图像可被识别的精度的前提下，我们可以保留奇异值较大的若干项，舍去奇异值较小的项即可。例如在上面的例子中，如果我们只保留奇异值分解的前50项，则需要存储的元素为1025*50=51250，和存储原始矩阵 $A$ 相比，存储量仅为后者的20%。

在图像处理领域，奇异值不仅可以应用在数据压缩上，还可以对图像去噪。如果一副图像包含噪声，我们有理由相信那些较小的奇异值就是由于噪声引起的。当我们强行令这些较小的奇异值为0时，就可以去除图片中的噪声。

部分来源链接：https://www.zhihu.com/question/22237507/answer/53804902

附上代码：

#!/usr/bin/python#  -*- coding:utf-8 -*-import numpy as npimport osfrom PIL import Imageimport matplotlib.pyplot as pltimport matplotlib as mplfrom pprint import pprintdef restore1(sigma, u, v, p):  # 奇异值、左特征向量、右特征向量    m = len(u)    n = len(v)    a = np.zeros((m, n))  # m × n 的0矩阵    for k in range(p):        uk = u[:, k].reshape(m, 1)  # m  行        vk = v[k].reshape(1, n)     # n 列        a += sigma[k] * np.dot(uk, vk)     #矩阵A    a[a < 0] = 0   #去燥    a[a > 255] = 255    # a = a.clip(0, 255)   这个一步到位，去燥    return np.rint(a).astype('uint8')#也可以下面这种定义# def restore2(sigma, u, v, K):  # 奇异值、左特征向量、右特征向量#     m = len(u)#     n = len(v)#     a = np.zeros((m, n))#     for k in range(K+1):#         for i in range(m):#             a[i] += sigma[k] * u[i][k] * v[k]#     a[a < 0] = 0#     a[a > 255] = 255#     return np.rint(a).astype('uint8')if __name__ == "__main__":    A = Image.open("C:\\Anaconda3\\lena.png", 'r')  #打开路径    print(A)    output_path = r'.\SVD_Output'                   #保存到文件夹SVD_Output    # if not os.path.exists(output_path):    #     os.mkdir(output_path)    a = np.array(A)    print(a.shape)    p = 50         #输出前50个奇异值到文件夹SVD_Output    u_r, sigma_r, v_r = np.linalg.svd(a[:,:,0])    u_g, sigma_g, v_g = np.linalg.svd(a[:,:,1])    u_b, sigma_b, v_b = np.linalg.svd(a[:,:,2])    plt.figure(figsize=(11, 9), facecolor='w')    mpl.rcParams['font.sans-serif'] = ['simHei']    mpl.rcParams['axes.unicode_minus'] = False    for k in range(1, p+1):        print(k)        R = restore1(sigma_r, u_r, v_r, k)        G = restore1(sigma_g, u_g, v_g, k)        B = restore1(sigma_b, u_b, v_b, k)        I = np.stack((R, G, B), axis=2)        Image.fromarray(I).save('%s\\svd_%d.png' % (output_path, k)) #将输出的图片保存        #也可以Image.fromarray(I).save("svd_"+str(k)+".png")        if k <= 12:                   #显示前12个            plt.subplot(3, 4, k)            plt.imshow(I)            plt.axis('off')            plt.title('奇异值个数：%d' % k)    plt.suptitle('SVD与图像分解', fontsize=20)    plt.tight_layout(0.3, rect=(0, 0, 1, 0.92))    # plt.subplots_adjust(top=0.9)    plt.show()

阅读全文

1 0