R上的对称非负矩阵分解

Question

我正在尝试基于以下公式在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
        }

}

Answer 1

这是对sym.nmf函数的修改，在此过程中进行了一些统计上重要的改进和速度提升。

当Ji在Ji-1的rel.tol百分比内时，

1. 会添加一个relative the (rel.tol)参数来中断循环。按照您设置的方式，只有当0 == 1或机器精度变得比拟合本身更可变时，才会停止循环。理论上，你的函数永远不会收敛。

1. 添加种子，因为可重复性很重要。沿着这条线，你可能会考虑用非负的双SVD来初始化，以获得领先的开始。但是，根据您的应用程序，这可能会使您的NMF陷入不能代表全局最小值的局部最小值，因此它可能是危险的。在我的例子中，我被锁定在一个类似奇异值分解的最小值中，NMF最终以一种完全不同于随机initialization.

的因式分解的状态收敛。

1. 添加最大迭代次数( number of iterations，max.iter)，因为有时您不希望运行一百万次迭代来达到您的容差阈值。

crossprod tcrossprod和

1. 函数中的 %*%替代基本函数。根据矩阵大小的不同，这会带来大约2倍的速度增益。

1. Reduce在中检查收敛的次数，因为在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()和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)或通过并行运行多个因子分解(因为通常需要多次重新启动才能发现健壮的解决方案)来实现进一步的收益。