一、奇异值分解简介
奇异值分解简称SVD（singular value decomposition），可以理解为：将一个比较复杂的矩阵用更小更简单的三个子矩阵的相乘来表示，这三个小矩阵描述了大矩阵重要的特性。SVD的用处有很多，比如：LSA（隐性语义分析）、推荐系统、数据降维、信号处理与统计等。
任何矩阵都可以使用SVD进行分解，对于一个MxN（M>=N）的矩阵M，存在以下的SVD分解：
∑是一个对角矩阵，其中的元素值就是奇异值，并且按照从大到小的顺序排列。
在很多情况下，前10%甚至更少的奇异值的平方就占全部奇异值平方的90%以上了，因此可以用前k个奇异值来近似描述矩阵：
而k的取值是由下面的公式决定：
其中percentage称为“奇异值平方和占比的阈值”，一般取90%，k是一个远小于m和n的值，这样也就达到了降维的目的。
二、MADlib的奇异值分解函数
MADlib的SVD函数可以对稠密矩阵和稀疏矩阵进行奇异值因式分解，并且还提供了一个稀疏矩阵的本地高性能实现函数。
1. 稠密矩阵的SVD函数
（1）语法
svd( source_table,
output_table_prefix,
row_id,
k,
n_iterations,
result_summary_table
);（2）参数
source_table：TEXT类型，源表名（稠密矩阵）。表含有一个row_id列标识每一行，从数字1开始。其它列包含矩阵的数据。可以使用两种稠密格式的任何一个，例如下面示例的2x2矩阵。
格式一：
row_id col1 col2
row1 1 1 0
row2 2 0 1
格式二：
row_id row_vec
row1 1 {1, 0}
row2 2 {0, 1}output_table_prefix：TEXT类型，输出表的前缀。
row_id：TEXT类型，每行的ID。
k：INTEGER类型，计算的奇异值的个数。
n_iterations（可选）：INTEGER类型，运行的迭代次数，必须在[k, 列维度]范围内，k是奇异值的个数。
result_summary_table（可选）：TEXT类型，存储结果摘要表名。
2. 稀疏矩阵的SVD函数
表示为稀疏格式的矩阵使用此函数。为了高效计算，在奇异值分解操作之前，输入矩阵会被转换为稠密矩阵。
（1）语法
svd_sparse( source_table,
output_table_prefix,
row_id,
col_id,
value,
row_dim,
col_dim,
k,
n_iterations,
result_summary_table
);（2）参数
source_table：TEXT类型，源表名（稀疏矩阵）。稀疏矩阵使用行列下标指示矩阵的每个非零条目，非常适合含有很多零元素的矩阵。如下面所示的4x7矩阵，除去零值只有6行。矩阵的维度由行、列的最大值推导出来。需要注意最后一行，即使是0也要包含这一行，因为它标识了矩阵的维度，并暗示了第4行与第7列全是0。
row_id | col_id | value
-------+--------+-------
1 | 1 | 9
1 | 5 | 6
1 | 6 | 6
2 | 1 | 8
3 | 1 | 3
3 | 2 | 9
4 | 7 | 0
(6 rows)
output_table_prefix：TEXT类型，输出表的前缀。
row_id：TEXT类型，包含行下标的列名。
col_id：TEXT类型，包含列下标的列名。
value：TEXT类型，包含值的列名。
row_dim：INTEGER类型，矩阵的行数。
col_dim：INTEGER类型，矩阵的列数。
k：INTEGER类型，计算的奇异值的个数。
n_iterations（可选）：INTEGER类型，运行的迭代次数，必须在[k, 列维度]范围内，k是奇异值的个数。
result_summary_table（可选）：TEXT类型，存储结果摘要表名。
3. 稀疏矩阵的本地实现SVD函数
此函数在计算SVD时使用本地稀疏表示，能够更高效地计算稀疏矩阵，适合高度稀疏的矩阵。
（1）语法
svd_sparse_native( source_table,
output_table_prefix,
row_id,
col_id,
value,
row_dim,
col_dim,
k,
n_iterations,
result_summary_table
);（2）参数
参数含义与svd_sparse函数相同。
4. 输出表
SVD函数的输出是以下三个表：
- 左奇异值矩阵：表名为<output_table_prefix>_u。
- 右奇异值矩阵：表名为<output_table_prefix>_v。
- 奇异值矩阵：表名为<output_table_prefix>_s。
左右奇异值向量表的格式为：
row_idINTEGER类型。每个特征值对应的ID。
row_vecFLOAT8[]类型。该row_id对应的特征向量元素，数组大小为k。
由于只有对角线元素是非零的，奇异值表采用稀疏表格式：
row_idINTEGER类型，第i个奇异值为i。
col_idINTEGER类型，第i个奇异值为i（与row_id相同）。
valueFLOAT8类型，奇异值。
其中的row_id和col_id都是从1开始。
结果摘要表名有以下列：
rows_usedINTEGER类型，SVD计算使用的行数。
exec_timeFLOAT8类型，计算SVD使用的总时间。
iterINTEGER类型，迭代运行次数。
recon_errorFLOAT8类型，这组正交基的质量得分（如近似精度）。计算公式为：
relative_recon_errorFLOAT8类型，相对质量分数。计算公式为：
5. 联机帮助
select madlib.svd();
-- 用法
select madlib.svd('usage');
-- 示例
select madlib.svd('example');三、奇异值分解函数实现推荐算法示例
用稀疏SVD函数解决上一篇介绍的音乐作品推荐问题。（参见“HAWQ + MADlib 玩转数据挖掘之（四）——低秩矩阵分解实现推荐算法”）
1. 建立输入表
（1）建立索引表
从前面的解释可以看到，推荐矩阵的行列下标分别表示用户和音乐作品。然而在业务系统中，userid和musicid很可能不是按从0到N的规则顺序生成的，因此需要建立矩阵下标值与业务表ID之间的映射关系，这里使用HAWQ的BIGSERIAL自增数据类型对应推荐矩阵的索引下标。
-- 用户索引表 drop table if exists tbl_idx_user; create table tbl_idx_user (user_idx bigserial, userid varchar(10)); -- 音乐索引表 drop table if exists tbl_idx_music; create table tbl_idx_music (music_idx bigserial, musicid varchar(10));
（2）建立用户行为业务表
drop table if exists source_data; create table source_data ( userid varchar(10), musicid varchar(10), val float8 );
（3）建立用户行为矩阵表
drop table if exists svd_data; create table svd_data ( row_id int, col_id int, val float8 );
2. 生成输入表数据
（1）生成用户行为业务表数据
insert into source_data values
('u1', 'm1', 5), ('u1', 'm6', -5),
('u2', 'm4', 3),
('u3', 'm3', 1), ('u3', 'm5', 2), ('u3', 'm7', 4),
('u4', 'm2', 4), ('u4', 'm3', 4), ('u4', 'm4', 3), ('u4', 'm7', -2),
('u5', 'm2', 5), ('u5', 'm3', -5), ('u5', 'm5', -5), ('u5', 'm7', 4), ('u5', 'm8', 3),
('u6', 'm3', 4), ('u6', 'm6', 3),
('u7', 'm2', -2), ('u7', 'm6', 5),
('u8', 'm2', -2), ('u8', 'm6', 5), ('u8', 'm8', 5),
('u9', 'm3', 1), ('u9', 'm5', 2), ('u9', 'm7', 4) ;（2）从业务表生成索引表数据
-- 用户表 insert into tbl_idx_user (userid) select distinct userid from source_data order by userid; -- 音乐表 insert into tbl_idx_music (musicid) select distinct musicid from source_data order by musicid;
（3）从业务表生成矩阵表数据
insert into svd_data select t1.user_idx, t2.music_idx, t3.val from tbl_idx_user t1, tbl_idx_music t2, source_data t3 where t1.userid = t3.userid and t2.musicid = t3.musicid;
之所以要用用户业务表作为数据源，是因为矩阵中包含所有有过打分行为的用户和被打过分的作品，但不包括没有任何打分行为相关的用户和作品。如果包含五行为记录的用户或作品，会在计算余弦相似度时出现除零错误或噪声数据。
3. 调用svd_sparse_native函数执行SVD
drop table if exists svd_u, svd_v, svd_s, svd_summary cascade;
select madlib.svd_sparse_native( 'svd_data', -- input table
'svd', -- output table prefix
'row_id', -- column name with row index
'col_id', -- column name with column index
'val', -- matrix cell value
9, -- number of rows in matrix
8, -- number of columns in matrix
7, -- number of singular values to compute
NULL, -- Use default number of iterations
'svd_summary' -- Result summary table
);说明：
- 选择svd_sparse_native函数的原因是测试数据比较稀疏，矩阵数据只有1/3（25/72），该函数效率较高。
- 这里给出的行、列、奇异值个数分别为9、8、7。svd_sparse_native函数要求行数大于等于列数，而奇异值个数小于等于列数，否则会报错。结果U、V矩阵的行数由实际的输入数据做决定，例如测试数据最大的行值为9，最大列值为8，则结果U矩阵的行数为9，V矩阵的行数为8，而不论行、列参数的值是多少。U、V矩阵的列数、S矩阵的行列数均由奇异值个数参数所决定。
4. 查看SVD结果
select array_dims(row_vec) from svd_u; select * from svd_s order by row_id, col_id; select array_dims(row_vec) from svd_v; select * from svd_summary;
输出如下：
dm=# select array_dims(row_vec) from svd_u;
array_dims
------------
[1:7]
[1:7]
[1:7]
[1:7]
[1:7]
[1:7]
[1:7]
[1:7]
[1:7]
(9 rows)
dm=# select * from svd_s order by row_id, col_id;
row_id | col_id | value
--------+--------+------------------
1 | 1 | 10.6650887159422
2 | 2 | 10.0400685494281
3 | 3 | 7.26197376834847
4 | 4 | 6.52278928434469
5 | 5 | 5.11307075598296
6 | 6 | 3.14838515537081
7 | 7 | 2.67251694708376
7 | 7 |
(8 rows)
dm=# select array_dims(row_vec) from svd_v;
array_dims
------------
[1:7]
[1:7]
[1:7]
[1:7]
[1:7]
[1:7]
[1:7]
[1:7]
(8 rows)
dm=# select * from svd_summary;
rows_used | exec_time (ms) | iter | recon_error | relative_recon_error
-----------+----------------+------+----------------+----------------------
9 | 6122.36 | 8 | 0.116171249851 | 0.0523917951113
(1 row)可以看到，结果U、V矩阵的维度分别是9 x 7和8 x 7，奇异值是一个7 x 7的对角矩阵。
5. 对比不同奇异值个数的近似度
让我们按公式3计算一下奇异值的比值，验证k设置为6、7时近似程度。需要将奇异值个数设置成与列数相同，执行SVD，再计算比值。
-- k=8
drop table if exists svd8_u, svd8_v, svd8_s, svd8_summary cascade;
select madlib.svd_sparse_native( 'svd_data', -- input table
'svd8', -- output table prefix
'row_id', -- column name with row index
'col_id', -- column name with column index
'val', -- matrix cell value
9, -- number of rows in matrix
8, -- number of columns in matrix
8, -- number of singular values to compute
NULL, -- Use default number of iterations
'svd8_summary' -- Result summary table
);
-- k=6
drop table if exists svd6_u, svd6_v, svd6_s, svd6_summary cascade;
select madlib.svd_sparse_native( 'svd_data', -- input table
'svd6', -- output table prefix
'row_id', -- column name with row index
'col_id', -- column name with column index
'val', -- matrix cell value
9, -- number of rows in matrix
8, -- number of columns in matrix
6, -- number of singular values to compute
NULL, -- Use default number of iterations
'svd6_summary' -- Result summary table
);对比近似度：
select * from svd6_summary;
select * from svd_summary;
select * from svd8_summary;
select s1/s3, s2/s3
from (select sum(value*value) s1 from svd6_s) t1,
(select sum(value*value) s2 from svd_s) t2,
(select sum(value*value) s3 from svd8_s) t3;结果如下：
dm=# select * from svd6_summary;
t3;
rows_used | exec_time (ms) | iter | recon_error | relative_recon_error
-----------+----------------+------+----------------+----------------------
9 | 3051.51 | 8 | 0.335700790666 | 0.151396899541
(1 row)
dm=# select * from svd_summary;
rows_used | exec_time (ms) | iter | recon_error | relative_recon_error
-----------+----------------+------+----------------+----------------------
9 | 6122.36 | 8 | 0.116171249851 | 0.0523917951113
(1 row)
dm=# select * from svd8_summary;
rows_used | exec_time (ms) | iter | recon_error | relative_recon_error
-----------+----------------+------+-------------------+----------------------
9 | 4182.38 | 8 | 1.52006310774e-15 | 6.85529638348e-16
(1 row)
dm=# select s1/s3, s2/s3
dm-# from (select sum(value*value) s1 from svd6_s) t1,
dm-# (select sum(value*value) s2 from svd_s) t2,
dm-# (select sum(value*value) s3 from svd8_s) t3;
?column? | ?column?
-------------------+-------------------
0.977078978809393 | 0.997255099805013
(1 row)可以看到，随着k值的增加，误差越来越小。在本示例中，奇异值个数为6、7的近似度分别为97.7%和99.7%。后面的计算都使用k=7的结果矩阵。
6. 基于用户的协同过滤算法UserCF生成推荐
（1）定义计算余弦相似度函数
余弦相似度计算公式为：
madlib.cosine_similarity()函数返回两个向量的余弦相似度。
（2）定义基于用户的协同过滤函数
create or replace function fn_user_cf(user_idx int)
returns table(r2 int, s float8, col_id int, val float8, musicid varchar(10)) as
$func$
select r2, s, col_id, val, musicid
from (select r2,s,col_id,val,row_number() over (partition by col_id order by col_id) rn
from (select r2,s,col_id,val
from (select r2,s
from (select r2,s,row_number() over (order by s desc) rn
from (select t1.row_id r1, t2.row_id r2, abs(madlib.cosine_similarity(v1, v2)) s
from (select row_id, row_vec v1 from svd_u where row_id = $1) t1,
(select row_id, row_vec v2 from svd_u) t2
where t1.row_id <> t2.row_id) t) t
where rn <=5 and s < 1) t1, svd_data t2
where t1.r2=t2.row_id and t2.val >=3) t
where col_id not in (select col_id from svd_data where row_id = $1)) t1,
tbl_idx_music t2
where t1.rn = 1 and t1.col_id = t2.music_idx
order by t1.s desc, t1.val desc limit 5;
$func$
language sql;说明：
select t1.row_id r1, t2.row_id r2, abs(madlib.cosine_similarity(v1, v2)) s
from (select row_id, row_vec v1 from svd_u where row_id = $1) t1,
(select row_id, row_vec v2 from svd_u) t2
where t1.row_id <> t2.row_id最内层查询调用madlib.cosine_similarity函数返回指定用户与其他用户的余弦相似度。
select r2,s,row_number() over (order by s desc) rn from ...
外面一层查询按相似度倒序取得排名。
select r2,s from ... where rn <=5 and s < 1
外面一层查询取得最相近的5个用户，同时排除相似度为1的用户，因为相似度为1说明两个用户的作品评分一模一样，而推荐的应该是用户没有打过分的作品。
select r2,s,col_id,val from ... where t1.r2=t2.row_id and t2.val >=3
外面一层查询取得相似用户打分在3及其以上的作品索引ID。
select r2,s,col_id,val,row_number() over (partition by col_id order by col_id) rn from ... where col_id not in (select col_id from svd_data where row_id = $1)
外面一层查询取得作品索引ID的排名，目的是去重，避免相同的作品推荐多次，并且过滤掉被推荐用户已经打过分的作品。
select r2, s, col_id, val, musicid ... where t1.rn = 1 and t1.col_id = t2.music_idx order by t1.s desc, t1.val desc limit 5;
最外层查询关联作品索引表取得作品业务主键，并按相似度和打分推荐前5个作品。
（3）定义接收用户业务ID的函数
create or replace function fn_user_recommendation(i_userid varchar(10))
returns table (r2 int, s float8, col_id int, val float8, musicid varchar(10)) as
$func$
declare
v_rec record;
v_user_idx int:=0;
begin
select user_idx into v_user_idx from tbl_idx_user where userid=i_userid;
for v_rec in (select * from fn_user_cf(v_user_idx)) loop
r2:=v_rec.r2;
s:=v_rec.s;
col_id:=v_rec.col_id;
val :=v_rec.val;
musicid:=v_rec.musicid;
return next;
end loop;
return;
end;
$func$
language plpgsql;通常输入的用户ID是业务系统的ID，而不是索引下标，因此定义一个接收业务系统的ID函数，内部调用fn_user_cf函数生成推荐。
（4）测试推荐结果
select * from fn_user_recommendation('u1');
select * from fn_user_recommendation('u3');
select * from fn_user_recommendation('u9'); 推荐结果如图1所示。可以看到，因为u3和u9的评分完全相同，相似度为1，所以为他们生成的推荐也完全相同。
7. 基于作品的协同过滤算法ItemCF生成推荐
（1）定义基于物品的协同过滤函数
create or replace function fn_item_cf(user_idx int)
returns table(r2 int, s float8, musicid varchar(10)) as
$func$
select t1.r2, t1.s, t2.musicid
from (select t1.r2,t1.s,row_number() over (partition by r2 order by s desc) rn
from (select t1.*, row_number() over (partition by r1 order by s desc) rn
from (select t1.row_id r1, t2.row_id r2, abs(madlib.cosine_similarity(v1, v2)) s
from (select row_id, row_vec v1
from svd_v
where row_id in (select col_id from svd_data where row_id=$1)) t1,
(select row_id, row_vec v2
from svd_v
where row_id not in (select col_id from svd_data where row_id=$1)) t2
where t1.row_id <> t2.row_id) t1) t1
where rn <=3) t1, tbl_idx_music t2
where rn = 1 and t1.r2 = t2.music_idx
order by s desc;
$func$
language sql;说明：
select t1.row_id r1, t2.row_id r2, abs(madlib.cosine_similarity(v1, v2)) s
from (select row_id, row_vec v1
from svd_v
where row_id in (select col_id from svd_data where row_id=$1)) t1,
(select row_id, row_vec v2
from svd_v
where row_id not in (select col_id from svd_data where row_id=$1)) t2
where t1.row_id <> t2.row_id 最内层查询调用madlib.cosine_similarity函数返回指定用户打过分的作品与没打过分的作品的相似度。
select t1.*, row_number() over (partition by r1 order by s desc) rn ...
外面一层查询按相似度倒序取得排名。
select t1.r2,t1.s,row_number() over (partition by r2 order by s desc) rn from ... where rn <=3
外面一层查询取得与每个打分作品相似度排前三的作品，并以作品索引ID分区，按相似度倒序取得排名，目的是去重，避免相同的作品推荐多次。
select t1.r2, t1.s, t2.musicid from ... where rn = 1 and t1.r2 = t2.music_idx order by s desc
最外层查询关联作品索引表取得作品业务主键并推荐。
（2）定义接收用户业务ID的函数
create or replace function fn_item_recommendation(i_userid varchar(10))
returns table (r2 int, s float8, musicid varchar(10)) as
$func$
declare
v_rec record;
v_user_idx int:=0;
begin
select user_idx into v_user_idx from tbl_idx_user where userid=i_userid;
for v_rec in (select * from fn_item_cf(v_user_idx)) loop
r2:=v_rec.r2;
s:=v_rec.s;
musicid:=v_rec.musicid;
return next;
end loop;
return;
end;
$func$
language plpgsql; 通常输入的用户ID是业务系统的ID，而不是索引下标，因此定义一个接收业务系统的ID函数，内部调用fn_item_cf函数生成推荐。
（3）测试推荐结果
select * from fn_item_recommendation('u1');
select * from fn_item_recommendation('u3');
select * from fn_item_recommendation('u9'); 推荐结果如图2所示。可以看到，因为u3和u9的评分作品完全相同，相似度为1，所以按作品相似度为他们生成的推荐也完全相同。
8. 为新用户寻找相似用户
假设一个新用户u10的评分向量为'{0,4,5,3,0,0,-2,0}'，要利用已有的奇异值矩阵找出该用户的相似用户。
（1）添加业务数据
insert into source_data values ('u10', 'm2', 4), ('u10', 'm3', 5), ('u10', 'm4', 3), ('u10', 'm7', -2);
insert into tbl_idx_user (userid)
select distinct userid
from source_data
where userid not in (select userid from tbl_idx_user)
order by userid; （2）确认从评分向量计算svd_u向量的公式
u10[1:8] x svd_v[8:7] x svd_s[7:7]^-1
（3）生成u10用户的向量表和数据
drop table if exists mat_u10;
create table mat_u10(row_id int, row_vec float8[]);
insert into mat_u10 values (1, '{0,4,5,3,0,0,-2,0}'); （4）根据计算公式，先将前两个矩阵相乘
drop table if exists mat_r_10;
select madlib.matrix_mult('mat_u10', 'row=row_id, val=row_vec',
'svd_v', 'row=row_id, val=row_vec',
'mat_r_10'); （5）根据公式，求奇异值矩阵的逆矩阵
drop table if exists svd_s_10; create table svd_s_10 as select row_id, col_id,1/value val from svd_s where value is not null;
（6）根据公式，将（3）、（4）两步的结果矩阵相乘
注意（3）的结果mat_r_10是一个稠密矩阵，（4）的结果svd_s_10是一个稀疏矩阵。
drop table if exists matrix_r_10;
select madlib.matrix_mult('mat_r_10', 'row=row_id, val=row_vec',
'svd_s_10', 'row=row_id, col=col_id, val=val',
'matrix_r_10'); （7）查询与u10相似的用户
select t1.row_id r1, t2.row_id r2, abs(madlib.cosine_similarity(v1, v2)) s
from (select row_id, row_vec v1 from matrix_r_10 where row_id = 1) t1,
(select row_id, row_vec v2 from svd_u) t2
order by s desc; 结果如图3所示。
可以看到，u10与u4的相似度高达99%，从原始的评分向量可以得到验证：
u4：'{0,4,4,3,0,0,-2,0}'
u10：'{0,4,5,3,0,0,-2,0}' （8）将结果向量插入svd_u矩阵
insert into svd_u
select user_idx, row_vec from matrix_r_10, tbl_idx_user where userid = 'u10';
遗留问题：
如果需要给用户行为表中没有的用户做推荐，或者推荐用户行为表没有的物品，如何进行归一化处理？
参考文献：
- Singular Value Decomposition：官方文档对奇异值分解的说明。
- 奇异值分解SVD简介及其在推荐系统中的简单应用：奇异值分解的数学推导，以及一个Python实现的推荐系统示例代码。
- SVD在推荐系统中的应用：比较详细的示例说明。
- 余弦计算相似度度量：说明余弦相似度计算的数学公式。
- 推荐系统_itemCF和userCF：说明基于用户和基于作品的两种协同过滤策略。