这是一个重做的sym.nmf
在此过程中进行了一些统计上重要的改进和速度增益。
-
Add a 相对耐受性 (rel.tol
) 参数,当 J[i] 在范围内时中断循环rel.tol
J[i-1] 的百分比。按照您的设置方式,只有当 0 == 1 或机器精度变得比拟合本身更加可变时,您才会停止循环。理论上,你的函数永远不会收敛。
-
Add a seed,因为再现性很重要。沿着这条线,您可能会考虑使用非负双 SVD 进行初始化以获得领先优势。但是,根据您的应用程序,这可能会将您的 NMF 推向局部最小值,而该局部最小值不能代表全局最小值,因此可能很危险。就我而言,我被锁定在类似 SVD 的最小值中,并且 NMF 最终收敛到完全不同于随机初始化的因式分解的状态。
-
Add a 最大迭代次数 (max.iter
),因为有时您不想运行一百万次迭代来达到您的容忍阈值。
-
替代在crossprod
and tcrossprod
基础功能%*%
功能。根据矩阵大小,这可实现约 2 倍的速度增益。
-
减少检查收敛的次数,因为计算残差信号W
减去后HH^T
占用了近一半的计算时间。您可以假设需要数百到数千次迭代才能收敛,因此只需每 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()
and 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
and tcrossprod
)或并行运行多个分解(因为通常需要多次重新启动才能发现可靠的解决方案)。