奇异值分解是一种矩阵分解的方法,大学线性代数里面也讲过奇异值分解的方法,因此这是一个为大家所熟知的算法。
1 SVD 时间复杂度分析
给定一个
m
×
n
m \times n
m×n 的矩阵
a
\boldsymbol{a}
a,按照下面公式做分解,其中
Σ
\Sigma
Σ
是一个
m
×
n
m \times n
m×n 的矩阵,除对角线以为其他元素都是 0,对角线上的值就是所谓的奇异值。
U
\boldsymbol{U}
U 是
m
×
n
m \times n
m×n 的酉矩阵,
V
\boldsymbol{V}
V 是 KaTeX parse error: Undefined control sequence: \n at position 1: \̲n̲ ̲\times n 的酉矩阵,即满足
U
T
∗
U
=
I
\boldsymbol{U} ^T * \boldsymbol{U} = \boldsymbol{I}
UT∗U=I
A
=
U
Σ
V
T
\boldsymbol{A} = \boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{V} ^T
A=UΣVT
- 将 A \boldsymbol{A} A 的转置和 A \boldsymbol{A} A 做矩阵乘法,这样就会得到一个 n × n n \times n n×n 的方阵 A T ∗ A \boldsymbol{A} ^T ∗ \boldsymbol{A} AT∗A,然后运用方阵特征分解,得到 KaTeX parse error: Undefined control sequence: \fbf at position 64: …= \lambda _i ∗ \̲f̲b̲f̲{v} _i
- 得到矩阵 A T ∗ A \boldsymbol{A} ^T ∗ \boldsymbol{A} AT∗A 的 n 个特征值和对应的 n 个特征向量 v \bf{v} v,将所有特征向量 v \bf{v} v 张成一个 KaTeX parse error: Undefined control sequence: \n at position 1: \̲n̲ ̲\times n 的矩阵 V \boldsymbol{V} V,就是前面公式里的 V \boldsymbol{V} V 矩阵,一般叫其中 V \boldsymbol{V} V 中的每个特征向量是 A \boldsymbol{A} A 的右奇异向量。
- 同样的对 A \boldsymbol{A} A 和 A \boldsymbol{A} A 的转置做乘法,就得到 KaTeX parse error: Undefined control sequence: \m at position 1: \̲m̲ ̲\times m 的方阵 A ∗ A T \boldsymbol{A} ∗ \boldsymbol{A} ^T A∗AT,运用方阵特征分解,得到 A ∗ A T ∗ u i = λ i ∗ u i \boldsymbol{A} ∗ \boldsymbol{A} ^T * \bf{u} _i = \lambda _i * \bf{u} _i A∗AT∗ui=λi∗ui
- 得到矩阵 A ∗ A T \boldsymbol{A} ∗ \boldsymbol{A} ^T A∗AT 的 m 个特征值和对应的 m 个特征向量 u \bf{u} u,将所有特征向量 u \bf{u} u 张成一个 KaTeX parse error: Undefined control sequence: \m at position 1: \̲m̲ ̲\times m 的矩阵 U \boldsymbol{U} U,就是前面公式里的 U \boldsymbol{U} U 矩阵,一般叫 U \boldsymbol{U} U 中的每个特征向量是 A \boldsymbol{A} A 的左奇异向量。
- 计算奇异矩阵
Σ
\Sigma
Σ
A = U ∗ Σ ∗ V T ⇒ A ∗ U = U ∗ Σ ∗ V T ∗ V ⇒ A ∗ V = U ∗ Σ ⇒ A ∗ v i = σ i ∗ u i ⇒ σ i = A ∗ v i / u i \boldsymbol{A} = \boldsymbol{U} ∗ \boldsymbol{\Sigma} ∗ \boldsymbol{V} ^T \Rightarrow \boldsymbol{A} ∗ \boldsymbol{U} = \boldsymbol{U} ∗ \boldsymbol{\Sigma} ∗ \boldsymbol{V} ^T ∗ \boldsymbol{V} \Rightarrow \boldsymbol{A} ∗ \boldsymbol{V} = \boldsymbol{U} ∗ \boldsymbol{\Sigma} \Rightarrow \boldsymbol{A} ∗ \bf{v} _i = \bf{\sigma} _i ∗ \bf{u} _i \Rightarrow \bf{\sigma} _i = \boldsymbol{A} ∗ \bf{v} _i / \bf{u} _i A=U∗Σ∗VT⇒A∗U=U∗Σ∗VT∗V⇒A∗V=U∗Σ⇒A∗vi=σi∗ui⇒σi=A∗vi/ui
这样得到奇异值 σ \sigma σ, 组成奇异值矩阵 Σ \Sigma Σ。
这是标准的按照线性代数理论进行分解的方法。由于矩阵乘法的实时间复杂度是 O ( n 3 ) O(n^3) O(n3),那么不难得知,按照线性代数理论进行奇异值分解的时间复杂度是 O ( max ( m , n ) 3 ) O(\max(m, n)^3) O(max(m,n)3)。
SVD 在工业上通常用于推荐系统,物料数和用户数通常最少也是千万级别,三次方下来就10,0000亿,一个很恐怖的时间复杂度,接下来介绍一个 SVD 的性质,然后继续分析怎么优化这个算法。
2 SVD 分解算法优化
2.1 奇异值分解性质
对于奇异值分解结果的中间矩阵,它对角线上的值称为奇异值,按照从大到小排序,并且奇异值下降特别快,正常情况下前10%的奇异值占有了全部奇异值的99%,那么可以用前 k 个奇异值近似表示原始矩阵。
A
m
×
N
=
U
m
×
m
∗
Σ
m
×
n
∗
V
n
×
n
T
≈
U
m
×
k
∗
Σ
k
×
k
∗
V
k
×
n
T
\boldsymbol{A} _{m \times N} = \boldsymbol{U} _{m \times m} * \Sigma _{m \times n} * \boldsymbol{V} _{n \times n} ^T \approx \boldsymbol{U} _{m \times k} * \Sigma _{k \times k} * \boldsymbol{V} _{k \times n} ^T
Am×N=Um×m∗Σm×n∗Vn×nT≈Um×k∗Σk×k∗Vk×nT
一般来说,这里
k
<
<
min
(
m
,
n
)
k << \min(m,n)
k<<min(m,n),可以达到数据压缩或者降维的作用。
2.2 截断奇异值分解
利用上面提到的奇异值分解的性质,可以借助截断奇异值分解来优化奇异值分解。
首先不直接分解矩阵,而是用机器学习的方法,直接去求第二个矩阵。定义损失函数
C
=
∑
(
i
,
j
)
∈
R
[
(
a
(
i
j
)
−
u
i
∗
v
j
T
)
2
+
λ
(
u
i
2
+
v
j
2
)
]
\boldsymbol{C} = \sum _{(i, j) \in \boldsymbol{R}} [(\bf{a}_(ij) - \bf{u} _i * \bf{v} _j ^T)^2 + \lambda (\bf{u} _i ^2 + \bf{v}_j ^2)]
C=(i,j)∈R∑[(a(ij)−ui∗vjT)2+λ(ui2+vj2)]
第一项使用平方差定义分解后和原始矩阵的 RMSE,其中 b f a i j bf{a}_{ij} bfaij 是原始矩阵的第 i 行第 j 列, u i \bf{u} _i ui 是推荐系统场景中的用户特征向量, v j \bf{v} _j vj 是物品特征向量,后面一项是正则化项。
有了损失函数就可以用 ALS 或者梯度下降法求解了。这里的时间复杂度是 O ( m n k ) O(mnk) O(mnk),而 k 是个很小的数,那么复杂度就是相当于是 O ( m ∗ n ) O(m*n) O(m∗n),这样在千万的数据上计算复杂度瞬间变成100亿了。
另外,在算法实现上,还可以使用并行计算的方式来进行 SVD 分解计算,限于篇幅这里暂时不做展开,Mars算法实践
3 变种 SVD 算法
3.1 改进的 SVD 算法
前面说到,奇异值分解常用在推荐系统中,最简单的方法就是直接把评分矩阵分解去做推荐,但是实际中发现不管是用户还是物品都存在一定偏差,比如有些用户给的评分就是比其他人高0.5,或者有些电影用户倾向于给高分,但是建模的时候没有考虑进去这种情况。那么把用户和物品的偏差考虑进来,那么就可以对 SVD 算法进行改进,评分 = 兴趣 + 偏见。
a u , i = u + b u + b i + U u ∗ V i T \bf{a} _{u, i} = \bf{u} + \bf{b}_u + \bf{b}_i + \boldsymbol{U}_u * \boldsymbol{V}_i^T au,i=u+bu+bi+Uu∗ViT
其中 b u \bf{b}_u bu 表示用户偏见, b i \bf{b}_i bi 表示物品偏见, u \bf{u} u 表示全局均值,损失函数为:
C = ∑ ( u , i ) ∈ R ( a u , i − b i − b u − U u ∗ V i T ) 2 + λ ( ∣ ∣ U u ∣ ∣ 2 + ∣ ∣ V i ∣ ∣ 2 + b u 2 b i 2 ) \boldsymbol{C} = \sum _{(u, i) \in \boldsymbol{R}} \left( \bf{a}_{u, i} - \bf{b}_i - \bf{b}_u - \boldsymbol{U}_u * \boldsymbol{V}_i^T \right) ^2 + \lambda \left( ||\boldsymbol{U}_u|| ^2 + ||\boldsymbol{V}_i|| ^2 + \bf{b}_u ^2 \bf{b}_i ^2 \right) C=(u,i)∈R∑(au,i−bi−bu−Uu∗ViT)2+λ(∣∣Uu∣∣2+∣∣Vi∣∣2+bu2bi2)
3.2 SVD++ 算法
实际使用中,除了用户或者物品偏见,还有一个问题就是行为数据中的评分数据很少,但是隐式行为数据有很多,把隐式行为数据建模进去从而缓解评分稀疏提高推荐效果,这就是 SVD++ 算法。SVD++ 在 SVD 中加入了用户对物品的隐式行为,SVD++ 的假设条件是评分 = 显式兴趣 + 隐式兴趣 + 偏见
a u , i = u + b u + b i + V i T ∗ ( U u + ∣ I u ∣ − 1 2 ∗ ∑ j ∈ I u y j ) \bf{a} _{u, i} = \bf{u} + \bf{b}_u + \bf{b}_i + \boldsymbol{V}_i^T * \left( \boldsymbol{U}_u + |\boldsymbol{I}_u| ^{-\frac{1}{2}} * \sum _{j \in \boldsymbol{I} _u} y_j \right) au,i=u+bu+bi+ViT∗ Uu+∣Iu∣−21∗j∈Iu∑yj
其中 I u \boldsymbol{I} _u Iu 是该用户有隐式行为的所有物品集合,而 y j y_j yj 是隐式评分电影 j 反应出的喜好值,其中对 I u \boldsymbol{I} _u Iu 取根号是一个经验值,这样就把系统中大量的隐式行为也建模到 SVD 算法中,虽然这么做对计算复杂度提高了不少,但是能够缓解评分稀疏提高推荐效果。同样的,损失函数为:
C = ∑ ( u , i ) ∈ R ( a u , i − b i − b u − V i T ∗ ( U u + min ( ∣ I u ∣ − 1 2 ∗ ∑ j ∈ I u y j ) ) ) 2 + λ ( ∑ u ( b i , u 2 + ∣ ∣ U u ∣ ∣ 2 ) + ∑ i ( b i 2 + ∣ ∣ V i ∣ ∣ 2 + ∣ ∣ y i ∣ ∣ 2 ) ) \boldsymbol{C} = \sum _{(u, i) \in \boldsymbol{R}} \left( \bf{a}_{u, i} - \bf{b}_i - \bf{b}_u - \boldsymbol{V}_i^T * (\boldsymbol{U}_u + \min(|\boldsymbol{I}_u| ^{-\frac{1}{2}} * \sum _{j \in \boldsymbol{I} _u} y_j)) \right) ^2 + \lambda \left( \sum _u (\bf{b}_{i, u}^2 + ||\boldsymbol{U} _u||^2) + \sum _i (\bf{b}_i^2 + ||\boldsymbol{V}_i|| ^2 + ||\bf{y}_i||^2) \right) C=(u,i)∈R∑ au,i−bi−bu−ViT∗(Uu+min(∣Iu∣−21∗j∈Iu∑yj)) 2+λ(u∑(bi,u2+∣∣Uu∣∣2)+i∑(bi2+∣∣Vi∣∣2+∣∣yi∣∣2))
3.3 timeSVD 算法
2010 年,Koren 发现在 Netflix 的数据中,个体用户的兴趣会随着时间转移,论文中称作 concepte drift (观念转移)。比如大家都知道的《大话西游》,刚上映票房很低,大家都给出的评分也不高,但是随着时间流逝,大家给出的评分越来越高。另外可能有些电影在某些特定的节日或者某天会获得高分,作者希望建立模型能捕捉到这些。
timeSVD 在 SVD 中加入了时间因素,也可以叫做有时序的SVD。timeSVD 的假设条件是评分 = 显式兴趣 + 时序兴趣 + 热点 + 偏见
按照原始论文的介绍来说,这个模型是为了对观念转移建模。从三个角度出发,首先是时间窗口概念,另外是实例权重,采用时间衰减,最后是集成学习的概念,引入多个基础预估方法。
static predictor 是最基础的预估方法:
b
u
,
i
(
t
)
=
u
+
b
u
+
b
i
\bf{b}_{u,i}(t) = \bf{u} + \bf{b}_u + \bf{b}_i
bu,i(t)=u+bu+bi
mov方案 mov predictor:
b
u
,
i
(
t
)
=
u
+
b
u
+
b
i
+
b
i
,
B
i
n
(
t
)
\bf{b}_{u,i}(t) = \bf{u} + \bf{b}_u + \bf{b}_i + \bf{b}_{i, Bin(t)}
bu,i(t)=u+bu+bi+bi,Bin(t)
这里对 b i \bf{b}_i bi 是 b i ( t ) = b f b i + b i , B i n ( t ) \bf{b}_i(t) = bf{b}_i + \bf{b}_{i, Bin(t)} bi(t)=bfbi+bi,Bin(t),分为Static和time changing两部分,引入物品的时间变化因素,对时间片进行划分。论文中是以10周为一片,共划分了30片,赋予一个 B i n ( t ) Bin(t) Bin(t),即1-30之间。时间片划分也可以是其他的,小的时间片可以有更好的性能,大的时间片能够让每片有更多的数据。
linear predictor,引入用户兴趣的变化,首先定义
d
e
v
u
(
t
)
=
s
i
g
n
(
t
−
t
u
)
∗
∣
t
−
t
u
∣
β
dev_u(t) = sign(t - t_u)*|t - t_u|^{\beta}
devu(t)=sign(t−tu)∗∣t−tu∣β
表示当前用户当前评分时间和平均评分时间的距离,论文中
β
=
0.4
β=0.4
β=0.4
然后对
b
u
\bf{b}_u
bu 引入时间变化到线性模型,多了一个需要训练的参数
α
\alpha
α
b
u
(
1
)
=
b
f
b
u
+
α
u
∗
d
e
v
u
(
t
)
\bf{b}_u^{(1)} = bf{b}_u + \alpha _u * dev_u(t)
bu(1)=bfbu+αu∗devu(t)
引入到公式中得到
b
u
,
i
(
t
)
=
u
+
b
u
+
α
u
∗
d
e
v
u
(
t
)
+
b
i
+
b
i
,
B
i
n
(
t
)
\bf{b}_{u, i}(t) = \bf{u} + \bf{b}_u + \alpha _u * dev_u(t) + \bf{b}_i + \bf{b}_{i, Bin(t)}
bu,i(t)=u+bu+αu∗devu(t)+bi+bi,Bin(t)
spline predictor 通过另一种方法引入用户兴趣变化的模型,区别于前面的线性模型,这是一个曲线式模型
b
u
,
i
(
t
)
=
u
+
b
u
+
∑
l
=
1
k
u
e
−
r
∣
t
−
t
u
∣
∗
b
t
,
l
u
∑
l
=
1
k
u
e
−
r
∣
t
−
t
u
∣
+
b
i
+
b
i
,
B
i
n
(
t
)
\bf{b}_{u, i}(t) = \bf{u} + \bf{b}_u + \frac{\sum _{l=1} ^{k_u} e^{-r|t-t^u|} * \bf{b}_{t, l}^{u}}{\sum _{l=1} ^{k_u} e^{-r|t-t^u|}} + \bf{b}_i + \bf{b}_{i, Bin(t)}
bu,i(t)=u+bu+∑l=1kue−r∣t−tu∣∑l=1kue−r∣t−tu∣∗bt,lu+bi+bi,Bin(t)
linear+ predictor 引入实时特定,比如用户某天心情,或者生日等,引入一个参数
b
u
,
t
\bf{b}_{u, t}
bu,t,表示每日的特定偏差
b
u
,
i
(
t
)
=
u
+
b
u
+
α
u
∗
d
e
v
u
(
t
)
+
b
i
+
b
u
,
t
+
b
i
,
B
i
n
(
t
)
\bf{b}_{u, i}(t) = \bf{u} + \bf{b}_u + \alpha _u * dev_u(t) + \bf{b}_i + \bf{b}_{u, t} + \bf{b}_{i, Bin(t)}
bu,i(t)=u+bu+αu∗devu(t)+bi+bu,t+bi,Bin(t)
spline+ predictor 曲线模型也引入,
b
u
,
i
(
t
)
=
u
+
b
u
+
∑
l
=
1
k
u
e
−
r
∣
t
−
t
u
∣
∗
b
t
,
l
u
∑
l
=
1
k
u
e
−
r
∣
t
−
t
u
∣
+
b
i
+
b
u
,
t
+
b
i
,
B
i
n
(
t
)
\bf{b}_{u, i}(t) = \bf{u} + \bf{b}_u + \frac{\sum _{l=1} ^{k_u} e^{-r|t-t^u|} * \bf{b}_{t, l}^{u}}{\sum _{l=1} ^{k_u} e^{-r|t-t^u|}} + \bf{b}_i + \bf{b}_{u, t} + \bf{b}_{i, Bin(t)}
bu,i(t)=u+bu+∑l=1kue−r∣t−tu∣∑l=1kue−r∣t−tu∣∗bt,lu+bi+bu,t+bi,Bin(t)
通过梯度下降法进行训练,损失函数为
C
=
∑
u
,
i
,
t
∈
K
(
r
u
,
i
(
t
)
−
u
−
b
u
−
α
u
∗
d
e
v
u
(
t
)
−
b
i
−
b
u
,
t
−
b
i
,
B
i
n
(
t
)
)
2
+
λ
∗
(
b
u
2
+
α
u
2
+
b
u
,
t
2
+
b
i
2
+
b
i
,
B
i
n
t
2
)
\boldsymbol{C} = \sum _{u, i, t \in \boldsymbol{K}} \left( \bf{r}_{u, i}(t) - \bf{u} - \bf{b}_u - \alpha _u * dev_u(t) - \bf{b}_i - \bf{b}_{u, t} - \bf{b}_{i, Bin(t)} \right) ^2 + \lambda * \left(\bf{b}_u^2 + \bf{\alpha}_u^2 + \bf{b}_{u, t}^2 + \bf{b}_i^2 + \bf{b}_{i, Bin{t}}^2 \right)
C=u,i,t∈K∑(ru,i(t)−u−bu−αu∗devu(t)−bi−bu,t−bi,Bin(t))2+λ∗(bu2+αu2+bu,t2+bi2+bi,Bint2)
原论文中对各个 predictor 用 RMSE 作为指标的结果
这里通过时序建模,用户矩阵的参数表示为:
u
u
(
t
)
[
k
]
=
p
u
,
k
+
α
u
,
k
∗
d
e
v
u
(
t
)
+
p
u
,
k
∗
k
=
1
,
⋯
,
f
\bf{u}_u(t)[k] = \bf{p}_{u, k} + \bf{\alpha}_{u, k} * dev_u(t) + \bf{p}_{u, k} * \bf{k} = 1, \cdots, f
uu(t)[k]=pu,k+αu,k∗devu(t)+pu,k∗k=1,⋯,f
这里 k 代表第 k 个 predictor,每一个 predictor 单独训练,最终得到如下公式,
a
u
,
i
=
u
+
b
i
(
t
)
+
b
u
(
t
)
+
V
i
T
∗
(
u
u
(
t
)
+
∣
I
u
∣
−
1
2
∗
∑
j
∈
I
u
y
j
)
\bf{a}_{u, i} = \bf{u} + \bf{b}_i(t) + \bf{b}_u(t) + \boldsymbol{V}_i^T * (\bf{u}_u(t) + |\boldsymbol{I}_u|^{-\frac{1}{2}} * \sum _{j \in \boldsymbol{I}_u} y_j)
au,i=u+bi(t)+bu(t)+ViT∗(uu(t)+∣Iu∣−21∗j∈Iu∑yj)
通过在 Netflix 数据集测试,对比三种算法,可以得到