奇异值分解

SVD（Singular Value Decomposition，奇异值分解）
numpy.linalg模块中的svd函数可以对矩阵进行奇异值分解。

分解的目标：

• 是一种因子分解运算，将一个矩阵分解为3个矩阵的乘积
• 3个矩阵: U, $\Sigma$ 和 V，其中U和V是正交矩阵，分别称为左奇异值、右奇异值，$\Sigma$为奇异值。

numpy的实现

numpy.linalg.svd(A,full_matrices=1,compute_uv=1)

参数：

• A是一个形如(m,n)矩阵
• full_matrices的取值是为0或者1，默认值为1,表示输出的U和V是否为全矩阵。
• compute_uv的取值是为0或者1，默认值为1，表示输出3个矩阵。为0的时候只输出$\Sigma$。

返回值：

返回3个矩阵U、$\Sigma$、V

• full_matrices=1时：U大小为(m,m)，$\Sigma$大小为(k)，k=min(m,n)，V大小为(n,n)。$A_{(m,n)}=U_{(m,m)}{\Sigma }_k{V}_{(n,n)}$
• full_matrices=0时：U大小为(m,k)，$\Sigma$大小为(k)，k=min(m,n)，V大小为(k,n)。$A_{(m,n)}=U_{(m,k)}{\Sigma }_k{V}_{(k,n)}$
• 其中**$\Sigma$**是对矩阵S的奇异值分解。$\Sigma$对角元素不为0，其他元素都为0，并且对角元素从大到小排列。$\Sigma$中的奇异值，一般排在后面的比较接近0，所以前r个奇异值可保留主要信息。

主要应用：压缩

$\Sigma$中有k个奇异值，但是由于排在后面的很多接近0，所以我们可以仅保留比较大的r个奇异值：
$A_{(m,n)}\approx U_{(m,r)}{\Sigma }_{(r,r)}{V}_{(r,n)}$
实际应用中，仅需保留着3个比较小的矩阵，就可近似表示A，不仅节省存储量，在计算的时候更是减少了计算量。SVD在信息检索（隐性语义索引）、图像压缩、推荐系统、金融等领域都有应用。

示例

import numpy as np

# 分解矩阵
A = np.mat("4 11 14 3;8 7 -2 0;5 2 8 9")
#D = np.mat(gray)
print('A:'), print(A)

# 使用svd函数分解矩阵
U,Sigma,V = np.linalg.svd(A,full_matrices=False)
print('U:'), print(U)
print('Sigma:'), print(Sigma)
print('V:'), print(V)
  
# 使用diag函数生成完整的奇异值矩阵。将分解出的3个矩阵相乘
print('np.diag(Sigma):')
print(np.diag(Sigma))
print('U * np.diag(Sigma) * V:')
print (U * np.diag(Sigma) * V)

A:
[[ 4 11 14  3]
 [ 8  7 -2  0]
 [ 5  2  8  9]]
U:
[[ 0.822 -0.049  0.567]
 [ 0.238  0.935 -0.264]
 [ 0.517 -0.352 -0.78 ]]
Sigma:
[21.893  9.921  7.436]
V:
[[ 0.355  0.536  0.693  0.325]
 [ 0.556  0.534 -0.542 -0.334]
 [-0.504  0.381  0.299 -0.716]]
np.diag(Sigma):
[[21.893  0.     0.   ]
 [ 0.     9.921  0.   ]
 [ 0.     0.     7.436]]
U * np.diag(Sigma) * V:
[[ 4. 11. 14.  3.]
 [ 8.  7. -2.  0.]
 [ 5.  2.  8.  9.]]

参数 full_matrices（2维数据的计算）：

a1 = np.random.randint(0,9,24).reshape(6,-1) # 随机生成0到9之间24个整数，并变换为3X4的数组
a2 = np.random.randint(0,9,24).reshape(4,-1) # 随机生成0到9之间24个整数，并变换为4X3的数组
print(a1),print(a2)

[[7 5 6 7]
 [5 5 6 3]
 [0 6 5 0]
 [2 1 1 7]
 [7 3 5 3]
 [4 7 2 4]]
[[0 2 5 0 3 5]
 [0 8 7 0 6 5]
 [0 2 4 3 8 3]
 [1 6 7 0 0 8]]

u1, s1, v1 = np.linalg.svd(a1, full_matrices=True)
u1.shape, s1.shape, v1.shape

((6, 6), (4,), (4, 4))

k = s1.shape[0]
np.dot(u1[:,:k] * s1, v1[:k,:])

[[ 7.,  5.,  6.,  7.],
   [ 5.,  5.,  6.,  3.],
   [ 0.,  6.,  5., -0.],
   [ 2.,  1.,  1.,  7.],
   [ 7.,  3.,  5.,  3.],
   [ 4.,  7.,  2.,  4.]]

u2, s2, v2 = np.linalg.svd(a2, full_matrices=False)
u2.shape, s2.shape, v2.shape

((4, 4), (4,), (4, 6))

k = s2.shape[0]
np.dot(u2[:,:k] * s2, v2[:k,:])

  [[ 0.,  2.,  5.,  0.,  3.,  5.],
   [-0.,  8.,  7., -0.,  6.,  5.],
   [-0.,  2.,  4.,  3.,  8.,  3.],
   [ 1.,  6.,  7., -0.,  0.,  8.]]

# 测试压缩结果
r = k-1
np.dot(u2[:,:r] * s2[:r], v2[:r,:]).round()

  [[ 0.,  2.,  5.,  1.,  3.,  5.],
   [ 0.,  8.,  7.,  0.,  6.,  5.],
   [-0.,  2.,  4.,  3.,  8.,  3.],
   [ 1.,  6.,  7., -1.,  0.,  8.]])

3维数据的计算：

a3 = np.random.randint(0,9,24).reshape(2,3,-1) # 随机生成0到9之间24个整数，并变换为3X3X4的数组
a3

   [[[6, 8, 6, 1],
     [0, 6, 7, 6],
     [3, 4, 3, 5]],
     
	[[8, 1, 7, 5],
    [2, 1, 6, 6],
    [7, 1, 2, 3]]]

u3, s3, v3 = np.linalg.svd(a3, full_matrices=False)
u3.shape, s3.shape, v3.shape

((2, 3, 3), (2, 3), (2, 3, 4))

np.matmul(u3*s3[:,None,:], v3)

  [[[6., 8., 6., 1.],
    [0., 6., 7., 6.],
    [3., 4., 3., 5.]],

   [[8., 1., 7., 5.],
    [2., 1., 6., 6.],
    [7., 1., 2., 3.]]]

# 测试压缩结果
r = 2
u3_ = u3[:,:,:r]
s3_ = s3[:,None,:r]
v3_ = v3[:,:r,:]
u3_.shape, s3_.shape, v3_.shape

((2, 3, 2), (2, 1, 2), (2, 2, 4))

np.matmul(u3_ * s3_, v3_).round()

  [[[6., 8., 6., 1.],
    [1., 6., 6., 7.],
    [2., 4., 4., 4.]],

   [[8., 1., 6., 6.],
    [2., 1., 6., 6.],
    [7., 1., 3., 2.]]]