奇异值分解原理(svd奇异值分解如何计算)

一、特征值分解(EVD)

1.特征值2.特征值分解过程 二、奇异值分解(SVD)

1.奇异值分解的形式2.奇异值分解的方法3.奇异值分解的价值 三、SVD分解在图像处理中的应用

3.1图像处理分析3.2SVD分解python代码3.3输出结果 四、用Eigen库实现SVD分解五、SVD分解的优点

  奇异值分解(Singular Value Decomposition，简称SVD)是在机器学习领域广泛应用的算法，它可以用于降维算法中的特征分解，使用SVD对矩阵进行分解, 能得到代表矩阵最明显变化的矩阵元素，是很多机器学习算法的基石。本文就对SVD的原理进行总结归纳，并举例说明其在图像处理中的应用。

  在介绍SVD分解之前，要先了解一下特征值分解，即EVD分解。关于特征值和特征向量的含义及理解，可以参考这篇文章。

  如果说一个向量 v v v是方阵 A n × n A_{n times n} An×n​的特征向量，将一定可以表示成下面的形式：
A v = λ v Av=lambda v Av=λv
写成矩阵的形式为：
A P = [ λ 1 λ 2 ⋱ λ n ] P AP= begin{bmatrix} lambda _1 \ & lambda _2 \ & & ddots \ & & & lambda _n \ end{bmatrix} P AP=⎣⎢⎢⎡​λ1​​λ2​​⋱​λn​​⎦⎥⎥⎤​P

其中， P n × n P_{ntimes n} Pn×n​为特征向量 v i n × 1 {v_i}_{ntimes 1} vi​n×1​组成的矩阵：
P = [ v 1 v 2 ⋯ v n ] P=begin{bmatrix} v_1 &v_2 & cdots &v_n end{bmatrix} P=[v1​​v2​​⋯​vn​​]

  如果矩阵 A A A是一个 m × m mtimes m m×m的实对称矩阵（即 A = A T A=A^T A=AT ），那么它可以被分解成如下的形式：
A = Q Σ Q T = Q [ λ 1 λ 2 ⋱ λ m ] Q T A=QSigma Q^T=Q begin{bmatrix} lambda _1 \ & lambda _2 \ & & ddots \ & & & lambda _m \ end{bmatrix} Q^T A=QΣQT=Q⎣⎢⎢⎡​λ1​​λ2​​⋱​λm​​⎦⎥⎥⎤​QT
其中 Q m × m Q_{mtimes m} Qm×m​为标准正交阵，即有 Q Q T = I , Q T Q = I QQ^T=I,Q^TQ=I QQT=I,QTQ=I, Σ Sigma Σ为对角矩阵，维度为 m × m mtimes m m×m。 λ i lambda _i λi​为特征值。

  特征值分解（EVD）是一个提取方阵特征很的方法，但是它只适用于方阵。而在机器学习的数据中，绝大部分的数据组成的矩阵都不是方阵，如一组数据有 m m m个样本，每个样本有 n n n个特征，往往这样的数据是样本的 m m m数量非常大，特征的 n n n数量有限，这样的数据组成的 矩阵 m × n mtimes n m×n是无法进行特征值分解的。奇异值分解就是用来分解这种任意矩阵的分解方法。
  奇异值分解是将一个非零的实数矩阵 A m × n A_{m times n} Am×n​分解成由三个是矩阵乘积形式的运算，即进行矩阵的因子分解：
A = U Σ V T A=USigma V^T A=UΣVT
其中， U U U为 m × m mtimes m m×m的单位正交阵， V V V为 n × n ntimes n n×n的单位正交阵，即有 U U T = I , V V T = I UU^T=I,VV^T=I UUT=I,VVT=I。 Σ Sigma Σ是 m × n mtimes n m×n维的对角矩阵其对角线上的数值即为奇异值，并且按照降序排列，如：
Σ = [ σ 1 σ 2 ⋱ ⋱ ] m × n Sigma= begin{bmatrix} sqrt{ sigma _1} \ & sqrt{sigma _2} \ & & ddots \ & & &ddots \ end{bmatrix}_{mtimes n} Σ=⎣⎢⎢⎡​σ1​ ​​σ2​ ​​⋱​⋱​⎦⎥⎥⎤​m×n​
   U Σ V T USigma V^T UΣVT为矩阵 A A A的奇异值分解， σ i sqrt{sigma _i} σi​ ​为矩阵的奇异值， U U U称为左奇异矩阵， V V V称为右奇异矩阵。各矩阵的维度分别为 U ∈ R m × m Uin R^{mtimes m} U∈Rm×m, Σ ∈ R m × n Sigma in R^{mtimes n} Σ∈Rm×n, V ∈ R n × n Vin R^{ntimes n} V∈Rn×n。

  为求解上述 U U U, Σ Sigma Σ, V V V，我们可以利用如下性质：
A A T = U Σ V T ( U Σ V T ) T = U Σ V T V Σ T U T = U Σ Σ T U T A T A = ( U Σ V T ) T U Σ V T = V Σ T U T U Σ V T = V Σ T Σ V T AA^T=USigma V^T(USigma V^T)^T=USigma V^TVSigma ^TU^T=USigma Sigma^TU^T\ A^TA=(USigma V^T)^TUSigma V^T=VSigma ^TU^TUSigma V^T=VSigma ^TSigma V^T AAT=UΣVT(UΣVT)T=UΣVTVΣTUT=UΣΣTUTATA=(UΣVT)TUΣVT=VΣTUTUΣVT=VΣTΣVT
其中 Σ Σ T ∈ R m × m Sigma Sigma^Tin R^{mtimes m} ΣΣT∈Rm×m, Σ T Σ ∈ R n × n Sigma^T Sigma in R^{ntimes n} ΣTΣ∈Rn×n,即：

Σ Σ T = [ σ 1 σ 2 ⋱ ⋱ ] m × m Sigma Sigma^T= begin{bmatrix} sigma _1 \ & sigma _2 \ & & ddots \ & & &ddots \ end{bmatrix}_{mtimes m} ΣΣT=⎣⎢⎢⎡​σ1​​σ2​​⋱​⋱​⎦⎥⎥⎤​m×m​

Σ T Σ = [ σ 1 σ 2 ⋱ ⋱ ] n × n Sigma^TSigma= begin{bmatrix} sigma _1 \ & sigma _2 \ & & ddots \ & & &ddots \ end{bmatrix}_{ntimes n} ΣTΣ=⎣⎢⎢⎡​σ1​​σ2​​⋱​⋱​⎦⎥⎥⎤​n×n​
可以看到 Σ Σ T Sigma Sigma^T ΣΣT和 Σ T Σ Sigma^T Sigma ΣTΣ的形式非常接近，进一步分析，我们可以发现 A A T AA^T AAT和 A T A A^TA ATA也是对称矩阵，那么就可以做特征值分解(EVD)。对 A A T AA^T AAT特征值分解，得到的特征矩阵即为 U U U 。对 A T A A^TA ATA特征值分解，得到的特征矩阵即为 V V V 。对 A A T AA^T AAT或 A T A A^TA ATA中的特征值开方，可以得到所有的奇异值。
由此可求得特征值为 σ 1 sqrt{ sigma _1} σ1​ ​ 、 σ 2 sqrt{ sigma _2} σ2​ ​ 、…、 σ k sqrt{ sigma _k} σk​ ​ 。所以矩阵 A A A：
A = U Σ V T = [ u 1 u 2 ⋯ u m   ] [ σ 1 σ 2 ⋱ ⋱ ] m × n [ v 1 v 2 ⋮ v n ] A=USigma V^T= begin{bmatrix} u_1 &u_2&cdots&u_m end{bmatrix}begin{bmatrix} sqrt{sigma _1} \ &sqrt{ sigma _2} \ & & ddots \ & & &ddots \ end{bmatrix}_{mtimes n} begin{bmatrix} v_1 \ v_2 \ vdots \ v_n \ end{bmatrix} A=UΣVT=[u1​​u2​​⋯​um​ ​]⎣⎢⎢⎡​σ1​ ​​σ2​ ​​⋱​⋱​⎦⎥⎥⎤​m×n​⎣⎢⎢⎢⎡​v1​v2​⋮vn​​⎦⎥⎥⎥⎤​
进一步化简得到
A = σ 1 u 1 v 1 T + σ 2 u 2 v 2 T + ⋯ + σ k u k v k T A=sqrt{sigma _1} u_1v_1^T+sqrt{sigma _2} u_2v_2^T+cdots+sqrt{sigma _k} u_kv_k^T A=σ1​ ​u1​v1T​+σ2​ ​u2​v2T​+⋯+σk​ ​uk​vkT​
即：
A = λ 1 u 1 v 1 T + λ 2 u 2 v 2 T + ⋯ + λ k u k v k T A=lambda _1 u_1v_1^T+lambda_2 u_2v_2^T+cdots+lambda _k u_kv_k^T A=λ1​u1​v1T​+λ2​u2​v2T​+⋯+λk​uk​vkT​
矩阵 A A A被分解为 k k k个小矩阵，每个矩阵都是 λ lambda λ乘以 u i v i T u_iv_i^T ui​viT​。此时可以看到，若 λ lambda λ取值较大，其对应的矩阵在矩阵 A A A的占比也大，所以取以去前面主要的 λ lambda λ来近似代表系统，近似的系统可以表示为：
A m × n = U m × m Σ m × n V n × n T ≈ U m × k Σ k × k V k × n T A_{mtimes n}=U_{mtimes m}Sigma _{mtimes n}V^T_{ntimes n} approx U_{mtimes k}Sigma _{ktimes k}V^T_{ktimes n} Am×n​=Um×m​Σm×n​Vn×nT​≈Um×k​Σk×k​Vk×nT​

  这种近似的应用就是SVD分解的主要价值，可以用来数据的降维或者减少数据的存储空间。例如原先矩阵 A ∈ R m × n Ain R^{mtimes n} A∈Rm×n需要存储的数据为 m × n mtimes n m×n个，而将其按照特征值分解后，假如只取第一个分量近似其矩阵的话，则存储的数据为 m + n + 1 m+n+1 m+n+1个，为保证其尽可能的近似于原矩阵，特征值可以选取前 k k k个，则其存储的数据为 k × ( m + n + 1 ) ktimes (m+n+1) k×(m+n+1)个，显然大大减少了存储空间。

  在图像处理中，往往要面对数据量比较大的图像。然而在实际的问题处理中，只需要提取图像的主要特征即可，不需要对其所有的数据进行处理，此时便需要提取主要的特征。例如下方的一张 700 × 500 × 3 700{times}500{times}3 700×500×3的图片，其数据量为 700 × 500 × 3 = 1050000 700{times}500{times}3=1050000 700×500×3=1050000个。对于这样的一张图片，其需要的存储空间为1050000个字节。运用SVD分解可以合理的减少其存储空间，从而大大简化后续的图像的其他操作的运算空间。下面介绍如何对一幅图像进行SVD分解，并得到其主要特征。

参考

import numpy as np
import matplotlib.image as mping
import matplotlib.pyplot as plt
import matplotlib as mpl


def image_svd(n, pic):
    a, b, c = np.linalg.svd(pic)#a,b,c= U,sigma,VT
    svd = np.zeros((a.shape[0], c.shape[1]))
    for i in range(0, n):
        svd[i, i] = b[i]
    img = np.matmul(a, svd)
    img = np.matmul(img, c)
    img[img >= 255] = 255
    img[0 >= img] = 0
    img = img.astype(np.uint8)
    return img


if __name__ == '__main__':
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    mpl.rcParams['axes.unicode_minus'] = False

    path = './cat.jpg'
    img = mping.imread(path)
    print(img.shape)

    r = img[:, :, 0]
    g = img[:, :, 1]
    b = img[:, :, 2]
    plt.figure(figsize=(50, 100))
    for i in range(1, 31):
        r_img = image_svd(i, r)
        g_img = image_svd(i, g)
        b_img = image_svd(i, b)
        pic = np.stack([r_img, g_img, b_img], axis=2)
        print(i)
        plt.subplot(5, 6, i)
        plt.title("图像的SVD分解，使用前 %d 个特征值" % (i))
        plt.axis('off')
        plt.imshow(pic)
    plt.suptitle("图像的SVD分解")
    plt.subplots_adjust()
    plt.show()

  通过分析结果图可以发现，不同数量的奇异值返回的图像效果不同。如果只取前5个奇异值，返回的图像较模糊，返回前15个奇异值的图像，则能较好的表现原图像，而返回30个奇异值生成的图像，则几乎与原图像误差，根据不同的应用场合，可以合理的选取奇异值，这样的大大简化了原先700个奇异值的运算。

用Eigen库实现SVD分解。

1.内存少
  奇异值分解矩阵中 Σ Sigma Σ里面的奇异值按从大到小的顺序排列，奇异值从大到小的顺序减小的特别快。在大多数情况下，前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上。也就是说，剩下的90%甚至99%的奇异值几乎没有什么作用。因此，可以用前面 个大的奇异值来近似描述矩阵，于是奇异值分解公式可以写成如下：
A m × n = U m × m Σ m × n V n × n T ≈ U m × k Σ k × k V k × n T A_{mtimes n}=U_{mtimes m}Sigma _{mtimes n}V^T_{ntimes n} approx U_{mtimes k}Sigma _{ktimes k}V^T_{ktimes n} Am×n​=Um×m​Σm×n​Vn×nT​≈Um×k​Σk×k​Vk×nT​
其中 k k k是一个远远小于 m m m 和 n n n 的数，右边的三个矩阵相乘的结果将会使一个接近 A A A的矩阵。如果 k k k越接近于 n n n ，则相乘的结果越接近于 A A A 。如果 k k k 的取值远远小于 n n n ，从计算机内存的角度来说，右边三个矩阵的存储内存要远远小于矩阵 A A A 。
2.SVD可以获取两个方向的主成分
  左奇异矩阵 U U U可以用于对行数的压缩；右奇异矩阵 V V V可以用于对列(即特征维度)的压缩。
3.滤波作用
  由于SVD分解最终是取得了多个权重较大的奇异值，因此其余较小的奇异值会被舍去。往往这些小的奇异值表征的是一些杂波，即不相关的干扰数据，运用SVD分解取得重要特征之后，恰好可以滤掉这些干扰数据，从而起到滤波的作用。