我的问题涉及从算法输出的奇异值.如果你做完整的SVD,为什么这些值不等于第一个k-奇异值?
下面我在R中有一个简单的实现.任何关于提高性能的建议都将受到赞赏.
rsvd = function(A, k=10, p=5) { n = nrow(A) y = A %*% matrix(rnorm(n * (k+p)), nrow=n) q = qr.Q(qr(y)) b = t(q) %*% A svd = svd(b) list(u=q %*% svd$u, d=svd$d, v=svd$v) } > set.seed(10) > A <- matrix(rnorm(500*500),500,500) > svd(A)$d[1:15] [1] 44.94307 44.48235 43.78984 43.44626 43.27146 43.15066 42.79720 42.54440 42.27439 42.21873 41.79763 41.51349 41.48338 41.35024 41.18068 > rsvd.o(A,10,5)$d [1] 34.83741 33.83411 33.09522 32.65761 32.34326 31.80868 31.38253 30.96395 30.79063 30.34387 30.04538 29.56061 29.24128 29.12612 27.61804
我认为你的算法是对Martinsson et al.算法的修改.如果我理解正确,这特别适用于低秩矩阵的近似.我可能错了.
差异很容易通过A(500)的实际等级与k(10)和p(5)的值之间的巨大差异来解释.此外,Martinsson等人提到p的值实际上应该大于k的选择值.
因此,如果我们将您的考虑因素考虑在内,请使用:
set.seed(10) A <- matrix(rnorm(500*500),500,500) # rank 500 B <- matrix(rnorm(500*50),500,500) # rank 50
我们发现,与原始svd算法相比,使用更大p值仍会导致巨大的加速时间.
> system.time(t1 <- svd(A)$d[1:5]) user system elapsed 0.8 0.0 0.8 > system.time(t2 <- rsvd(A,10,5)$d[1:5]) user system elapsed 0.01 0.00 0.02 > system.time(t3 <- rsvd(A,10,30)$d[1:5]) user system elapsed 0.04 0.00 0.03 > system.time(t4 <- svd(B)$d[1:5] ) user system elapsed 0.55 0.00 0.55 > system.time(t5 <-rsvd(B,10,5)$d[1:5] ) user system elapsed 0.02 0.00 0.02 > system.time(t6 <-rsvd(B,10,30)$d[1:5] ) user system elapsed 0.05 0.00 0.05 > system.time(t7 <-rsvd(B,25,30)$d[1:5] ) user system elapsed 0.06 0.00 0.06
但我们看到,对较低秩矩阵使用较高的p确实给出了更好的近似.如果我们让k也接近等级,那么真实解与近似之间的差异就变成了appx. 0,而速度增益仍然很大.
> round(mean(t2/t1),2) [1] 0.77 > round(mean(t3/t1),2) [1] 0.82 > round(mean(t5/t4),2) [1] 0.92 > round(mean(t6/t4),2) [1] 0.97 > round(mean(t7/t4),2) [1] 1
总的来说,我相信可以得出结论:
应该选择> p,所以p> k(如果我是对的话,Martinsson称之为l)
> k不应与等级(A)有太大差别
>对于低秩矩阵,结果通常更好.
Optimalization:
就我而言,这是一种巧妙的方式.实际上,我真的找不到更优化的方式.我唯一可以说的是建议t(q)%*%A被建议不要.应该使用crossprod(q,A),这应该更快一点.但在你的例子中,差异是不存在的.