奇异值分解

发表于 | 分类于

奇异值分解在机器学习领域应用广泛，结合MLiA一书中的例子，本文介绍它在图像压缩和推荐系统中的应用。

基本原理

SVD压缩的原理

奇异值分解： \[ X_{m \times n} = U_{m \times m} \Sigma_{m \times n} V^T_{n \times n} \] 选择前\(r\)个奇异值，可以近似原来的数据矩阵 \[ X_{m \times n} \approx U_{m \times r} \Sigma_{r \times r} V^T_{r \times n} \] 由于\(r\)可能远小于\(n\)，也即\(m \times n << m \times r + r \times r + r \times n\)，因此可以对数据进行压缩。

奇异值数量的选取

计算前\(r\)个特征值的能量占总特征值能量的比例: \[ \frac{\sum_{i=1}^{r} \lambda_i^2}{\sum_{i=1}^{n} \lambda_i^2} \] 设定阈值，比如要保留\(90%\)的能量，则使上述比例大于\(0.9\)，最小的\(r\)值作为主成分的数量。

实例1: SVD在推荐系统中的应用

基于协同过滤的推荐系统

协同过滤（collaborative filtering）是通过将用户和其他用户的数据进行对比来实现推荐的。将采集的数据组织成如下表格的形式：

电影1 电影2 电影3
用户1 2 0 3
用户2 0 1 4
用户3 3 4 0

不是根据电影本身的属性，而是通过用户给电影的评分来计算不同电影之间的相似度，如果某部电影和用户看过的电影之间的相似度高，则推荐这部电影给这名用户。

相似度的计算

计算向量\(\mathbf{u}=(u_1, u_2)\)和向量\(\mathbf{v}=(v_1, v_2)\) 1. 欧氏距离： \[ \sqrt{(u_1-v_1)^2+(u_2-v_2)^2} \]

  1. 皮尔逊相关系数（Pearson correlation）： \[ 0.5+0.5*\frac{cov(u,v)}{\sqrt{cov(u,u)*cov(v,v)}} \] 其中协方差\(cov(u,v)=\frac{1}{m} \sum_{i=1}^{m} (u_i-\bar{u})(v_i-\bar{v})\)

  2. 余弦相似度（cosine similarity）： \[ 0.5 + 0.5*\frac{\mathbf{u} \cdot \mathbf{v}}{\|u\| \|v\|} \]

三种相似度范围均在\([0, 1]\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
from numpy import *
from numpy import linalg as la

def euclidSim(inA, inB):
    return 1.0/(1.0+la.norm(inA - inB))

def pearsSim(inA, inB):
    if len(inA) < 3: return 1.0
    return 0.5+0.5*corrcoef(inA, inB, rowvar=0)[0][1]

def cosSim(inA, inB):
    num = float(inA.T*inB)
    denom = la.norm(inA)*la.norm(inB)
    return 0.5+0.5*(num/denom)
1
2
3
4
5
6
7
8
def loadExData():
    return mat([[0, 0, 0, 2, 2],
           [0, 0, 0, 3, 3],
           [0, 0, 0, 1, 1],
           [1, 1, 1, 0, 0],
           [2, 2, 2, 0, 0],
           [5, 5, 5, 0, 0],
           [1, 1, 1, 0, 0]])
1
2
myMat = loadExData()
print(euclidSim(myMat[:,0],myMat[:,4]),euclidSim(myMat[:,4],myMat[:,4]))
0.129731907557 1.0
1
print(pearsSim(myMat[:,0],myMat[:,4]),pearsSim(myMat[:,4],myMat[:,4]))
0.205965381738 1.0
1
print(cosSim(myMat[:,0],myMat[:,4]),cosSim(myMat[:,4],myMat[:,4]))
0.5 1.0

基于物品的相似度

到底应该采用基于用户的相似度还是基于物品的相似度，取决于用户或者物品的数量。对于大部分产品导向的推荐引擎而言，用户的数量往往大于物品的数量，因此大多采用基于物品的相似度。

推荐系统的评价

将某些已知的评分去掉，然后对他们进行预测，最后计算预测值和真实值之间的差异。通常用于评价推荐系统的指标是最小均方根误差（RMSE），即首先计算均方误差的平均值，然后取其平方根，也即L2模。

推荐系统的工作流程

  1. 寻找用户没有评分的物品
  2. 对每个物品预计一个可能的评分（参考相似度计算结果）。
  3. 对这些评分从高到低排序，返回前\(N\)个物品
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def standEst(dataMat, user, simMeas, item):
    """
    用户对物品的估计评分值
    dataMat: 数据矩阵, (m, n)
    user: 用户编号
    simMeas: 相似度计算方法
    item: 物品编号
    """
    n = shape(dataMat)[1] # 物品的种类数
    simTotal = 0.0; ratSimTotal = 0.0 #初始化
    for j in range(n):
        userRating = dataMat[user, j] # 该用户对其他物品的打分情况
        if userRating == 0: continue
        # 寻找两个用户都评级的物品
        overLap = nonzero(logical_and(dataMat[:,item].A>0, dataMat[:,j].A>0))[0]
        if len(overLap) == 0: similarity = 0 # 如果两个用户没有都评级的物品
        else: similarity = simMeas(dataMat[overLap, item], dataMat[overLap,j]) # 计算两个物品的相似度
        simTotal += similarity # 计算总的相似度
        ratSimTotal += similarity*userRating # 计算相似度加权后的用户评分
    if simTotal == 0: return 0
    else: return ratSimTotal/simTotal #　计算归一化之后的评分结果
    
def recommend(dataMat, user, N=3, simMeas=cosSim, estMethod=standEst):
    """
    dataMat: 数据矩阵
    user: 用户编号
    N: 推荐的前N个结果
    simMeas: 相似度计算结果
    estMethod: 评分估计方法
    """
    unratedItems = nonzero(dataMat[user,:].A==0)[1] # 寻找所有未评分的物品
    if len(unratedItems) == 0: return 'you rated everything'
    itemScores = []
    # 对所有未评分物品循环
    for item in unratedItems:
        estimatedScore = estMethod(dataMat, user, simMeas, item) #估计评分值
        itemScores.append((item, estimatedScore)) # 将物品编号和估计的评分值存储下来
    return sorted(itemScores, key=lambda jj: jj[1], reverse=True)[:N] #　输出评分值排名前N的物品
1
2
3
4
myMat = loadExData()
myMat[0,1] = myMat[0,0] = myMat[1, 0] = myMat[2, 0] = 4
myMat[3,3] = 2
myMat
matrix([[4, 4, 0, 2, 2],
        [4, 0, 0, 3, 3],
        [4, 0, 0, 1, 1],
        [1, 1, 1, 2, 0],
        [2, 2, 2, 0, 0],
        [5, 5, 5, 0, 0],
        [1, 1, 1, 0, 0]])
1
2
# 利用余弦相似度计算未评分物品对用户２的推荐系数
recommend(myMat, 2)
[(2, 2.5), (1, 2.0243290220056256)]
1
2
# 利用欧式距离计算未评分物品对用户２的推荐系数
recommend(myMat, 2, simMeas=euclidSim)
[(2, 3.0), (1, 2.8266504712098603)]
1
2
# 利用Pear相关系数计算未评分物品对用户２的推荐系数
recommend(myMat, 2, simMeas=pearsSim)
[(2, 2.5), (1, 2.0)]

可以看出，三种计算相似度的方法都提示给用户２推荐物品的顺序是２，然后是１

利用SVD提高推荐效果

当数据十分稀疏时，可以采用SVD降维。需要注意的是：用户的个数也就是特征的个数。通过计算前\(N\)个特征值的能量占总能量的比例估计需要保留的特征值。

1
2
3
4
5
6
7
8
9
10
11
12
13
# 一个数据更加稀疏的例子
def loadExData2():
    return mat([[0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 5],
           [0, 0, 0, 3, 0, 4, 0, 0, 0, 0, 3],
           [0, 0, 0, 0, 4, 0, 0, 1, 0, 4, 0],
           [3, 3, 4, 0, 0, 0, 0, 2, 2, 0, 0],
           [5, 4, 5, 0, 0, 0, 0, 5, 5, 0, 0],
           [0, 0, 0, 0, 5, 0, 1, 0, 0, 5, 0],
           [4, 3, 4, 0, 0, 0, 0, 5, 5, 0, 1],
           [0, 0, 0, 4, 0, 4, 0, 0, 0, 0, 4],
           [0, 0, 0, 2, 0, 2, 5, 0, 0, 1, 2],
           [0, 0, 0, 0, 5, 0, 0, 0, 0, 4, 0],
           [1, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0]])
1
2
U, Sigma, VT = la.svd(loadExData2())
Sigma
array([ 15.77075346,  11.40670395,  11.03044558,   4.84639758,
         3.09292055,   2.58097379,   1.00413543,   0.72817072,
         0.43800353,   0.22082113,   0.07367823])
1
2
3
4
5
6
7
def num_sigma(Sigma):
    Sig2 = Sigma**2
    all_sig2 = sum(Sig2)
    sum_sig2 = 0.0
    for i in range(len(Sigma)):
        sum_sig2 += Sig2[i]
        if (sum_sig2 >= sum(Sig2)*0.9): return i+1
1
num_sigma(Sigma)
3
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def svdEst(dataMat, user, simMeas, item):
    """
    用户对物品的估计评分值
    dataMat: 数据矩阵, (m, n)
    user: 用户编号
    simMeas: 相似度计算方法
    item: 物品编号
    """
    n = shape(dataMat)[1] # 物品的种类数
    simTotal = 0.0; ratSimTotal = 0.0 #初始化
    U, Sigma, VT = la.svd(dataMat) # 奇异值分解
    Sig4 = mat(eye(4)*Sigma[:4]) # 转换成矩阵的形式
    xformedItems = dataMat.T * U[:,:4] * Sig4.I # xformedItem:(n, 4)
    for j in range(n):
        userRating = dataMat[user, j] # 用户打分情况
        if userRating == 0: continue
        similarity = simMeas(xformedItems[item,:].T, xformedItems[j,:].T) # 计算两个物品的相似度
        simTotal += similarity # 计算总的相似度
        ratSimTotal += similarity*userRating # 计算相似度加权后的用户评分
    if simTotal == 0: return 0
    else: return ratSimTotal/simTotal #　计算归一化之后的评分结果
1
2
myMat = loadExData2()
recommend(myMat, 1, estMethod=svdEst)
[(4, 3.3447149384692283), (7, 3.3294020724526963), (9, 3.328100876390069)]
1
recommend(myMat, 1, estMethod=svdEst, simMeas=pearsSim)
[(4, 3.3469521867021736), (9, 3.3353796573274699), (6, 3.307193027813037)]

实例2: 基于SVD的图像压缩

利用SVD进行图像压缩的原理十分简单，在本文前面已经介绍过了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def printMat(inMat, thresh=0.8):
    for i in range(32):
        for k in range(32):
            if float(inMat[i,k]) > thresh: print(1,end=' ')
            else: print(0,end=' ')
        print('\n')

def imgCompress(numSV=3, thresh=0.8):
    myl = []
    for line in open('0_5.txt').readlines():
        newRow = []
        for i in range(32):
            newRow.append(int(line[i]))
        myl.append(newRow)
    myMat = mat(myl)
    print("*** Original matrix ****")
    printMat(myMat, thresh)
    U, Sigma, VT = la.svd(myMat)
    SigRecon = mat(zeros((numSV, numSV)))
    for k in range(numSV):
        SigRecon[k, k] = Sigma[k]
    reconMat = U[:, :numSV]*SigRecon*VT[:numSV,:]
    print("*** Reconstructed matrix using %d singular values ****" %numSV)
    printMat(reconMat, thresh)
1
imgCompress(2)
*** Original matrix ****
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 

0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 

0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 

0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 

0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 

0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 

0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 

*** Reconstructed matrix using 2 singular values ****
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

原来的图像大小\(32*32=1024\)，现在变为\(32*2+32*2+2=130\)，获得几乎\(10\)倍的压缩比

结论及讨论

结论

SVD是一种十分强大且实用的工具，在降维、数据压缩、数据可视化等方面都有重要的应用。

讨论

  1. 上面推荐引擎的例子中，SVD可以函数外面执行，从而不需要每次估计一个物品的评分都执行一次分解；另外，为了提高推荐引擎的效率，可以通过离线计算并保存相似度得分；为了节省内存和计算开销，可以只存储采集矩阵的非零元素。
  2. 推荐引擎的另一个问题是冷启动，指的是新物品出现后没有足够的用户评分信息，无法判断每个用户对该物品的喜好。这个问题可以看做搜索问题，可能需要推荐物品的属性。同时，也可以将这些属性作为相似度计算所需要的数据，这常称为基于内容的推荐。

参考资料

  1. Peter Harrington, Machine learning in action, Manning Publications, 2012
  2. 周志华，机器学习，清华大学出版社， 2016