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 package NMF 可以做到这一点,如果你不想尝试重新发明*** 如果你对方法感兴趣,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),因为有时您不想运行一百万次迭代来达到容差阈值。

    crossprodtcrossprod 函数替换基本 %*% 函数。根据矩阵大小,这将获得大约 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 以进一步提高速度。并行化还可以在目标函数内(并行化crossprodtcrossprod)或通过并行运行多个因式分解(因为通常需要多次重新启动来发现稳健的解决方案)来获得进一步的收益。

【讨论】:

以上是关于R中的对称非负矩阵分解的主要内容,如果未能解决你的问题,请参考以下文章

推荐系统中矩阵分解算法-funkSVD和ALS

推荐系统笔记:基于非负矩阵分解的协同过滤

NMF非负矩阵分解初探

基于R语言利用NMF(非负矩阵分解)替代层次聚类进行肿瘤分型

推荐算法——非负矩阵分解(NMF)

SVD(奇异值分解)+NMF(非负矩阵分解)