R中的对称非负矩阵分解
Posted
技术标签:
【中文标题】R中的对称非负矩阵分解【英文标题】:Symmetric Non negative Matrix Factorization in R 【发布时间】:2016-02-11 18:41:19 【问题描述】:我正在尝试根据以下公式在 R 中实现 NMF: H 最初是猜测,然后根据这个公式迭代更新。我写了这段代码,但它需要像以往一样执行。我怎样才能重写这段代码? W是相似度矩阵。
sym.nmf <- function ( W )
N <- ncol(W)
set.seed(1234)
H <- matrix(runif(N * k, 0, 1),N,k)
J1 <- 0
while (0 < 1)
HT <- t(H)
A <- W %*% H
B <- H %*% HT %*% H
H <- 0.5 * ( H * ( 1 + ( A / B )))
J = W - (H %*% t(H))
J = sum (J^2)
if ( (J1 != 0 ) && (J > J1) )
return (H1)
H1 <- H
J1 <- J
【问题讨论】:
有一个 R packageNMF
可以做到这一点,如果你不想尝试重新发明***
如果你对方法感兴趣,you could check out their implementation。
@TomNash 不幸的是,这是一种新方法,没有任何实现
我没有看到任何明显慢的东西,所以如果你想要更快的速度,可能是时候看看rcpp
或类似的东西了。
@Gregor 矩阵 W 的维数较大时速度较慢。就我而言,它是 1500* 1500。
【参考方案1】:
这是对 sym.nmf
函数的修改,在此过程中进行了一些重要的统计改进和速度提升。
添加相对容差 (rel.tol
) 参数以在 J[i] 位于 J[i-1] 的rel.tol
百分比范围内时中断循环。按照您的设置方式,您只会在 0 == 1 或机器精度变得比拟合本身更具可变性时停止循环。理论上,你的函数永远不会收敛。
添加种子,因为可重复性很重要。沿着这条线,您可能会考虑使用非负双 SVD 进行初始化以抢占先机。但是,根据您的应用程序,这可能会使您的 NMF 进入不代表全局最小值的局部最小值,因此可能很危险。在我的例子中,我被锁定在一个类似于 SVD 的最小值中,并且 NMF 最终收敛到一个完全不同于随机初始化分解的状态。
添加最大迭代次数 (max.iter
),因为有时您不想运行一百万次迭代来达到容差阈值。
用 crossprod
和 tcrossprod
函数替换基本 %*%
函数。根据矩阵大小,这将获得大约 2 倍的速度增益。
减少收敛检查的次数,因为在减去HH^T
之后计算W
中的残差信号需要将近一半的计算时间。您可以假设需要数百到数千次迭代才能收敛,因此只需每 100 个周期检查一次收敛。
更新功能:
sym.nmf <- function (W, k, seed = 123, max.iter = 10000, rel.tol = 1e-10)
set.seed(seed)
H <- matrix(runif(ncol(W) * k, 0, 1),ncol(W),k)
J <- c()
for(i in 1:max.iter)
H <- 0.5*(H*(1+(crossprod(W,H)/tcrossprod(H,crossprod(H)))))
# check for convergence every 100 iterations
if(i %% 100 == 0)
J <- c(J,sum((W - tcrossprod(H))^2))
plot(J, xlab = "iteration", ylab = "total residual signal", log = 'y')
cat("Iteration ",i,": J =",tail(J)[1],"\n")
if(length(J) > 3 && (1 - tail(J, 1)/tail(J, 2)[1]) < rel.tol)
return(H)
if(i == max.iter)
warning("Max.iter was reached before convergence\n")
return(H)
目标函数也可以隔离,Rfast也可以用于Rfast::Crossprod()
和Rfast::Tcrossprod()
的并行计算。
sym.nmf <- function (W, k, seed = 123, max.iter = 100, rel.tol = 1e-10)
set.seed(seed)
require(Rfast)
H <- matrix(runif(ncol(W) * k, 0, 1),ncol(W),k)
J <- c()
for(i in 1:max.iter)
H <- 0.5 * fit_H(W,H, num.iter = 100)
J <- c(J,sum((W - tcrossprod(H))^2))
plot(J, xlab = "iteration", ylab = "total residual signal", log = 'y')
cat("Iteration ",i,": J =",tail(J, n = 1),"\n")
if(length(J) > 3 && (1 - tail(J, 1)/tail(J, 2)[1]) < rel.tol)
return(H)
if(i == max.iter)
warning("Max.iter was reached before convergence\n")
return(H)
fit_H <- function(W,H, num.iter)
for(i in 1:num.iter)
H <- 0.5*(H*(1+(Rfast::Crossprod(W,H)/Rfast::Tcrossprod(H,Rfast::Crossprod(H,H)))))
H
现在可以将此目标函数转换为 Rcpp 以进一步提高速度。并行化还可以在目标函数内(并行化crossprod
和tcrossprod
)或通过并行运行多个因式分解(因为通常需要多次重新启动来发现稳健的解决方案)来获得进一步的收益。
【讨论】:
以上是关于R中的对称非负矩阵分解的主要内容,如果未能解决你的问题,请参考以下文章