# R/HALSacc.R In hNMF: Hierarchical Non-Negative Matrix Factorization

```# Project: hNMF_git
#
# Author: nsauwen
# Accelerated hierarchical alternating least squares NMF. This function has been
# written to be compatible with the 'nmf' function from the 'NMF' package.
# It can be added to the list of NMF algorithms available to the 'nmf' function
# using the 'SetNMFMethod' function. For a reference to the method, see N. Gillis,
# Nonnegative matrix factorization: complexity, algorithms and applications
# [Section 4.2, Algo. 6], PhD thesis, Université catholique de Louvain, February 2011.
###############################################################################

#' Accelerated hierarchical alternating least squares NMF. For a reference to the method, see N. Gillis,
#' Nonnegative matrix factorization: complexity, algorithms and applications
#' [Section 4.2, Algo. 6], PhD thesis, Université catholique de Louvain, February 2011.
#' @param X Input data matrix, each column represents one observation
#' and the rows correspond to the different features
#' @param nmfMod Valid NMF model, containing initialized factor matrices
#' (in accordance with the NMF package definition)
#' @param alpha Nonnegative parameter of the accelerated method
#' @param maxiter Maximum number of iterations
#' @param checkDivergence currently not in use, to be implemented
#' @return Resulting NMF model (in accordance with the NMF package definition)
#' @importFrom NMF basis coef
#' @author nsauwen
#' @export
HALSacc <- function (X,nmfMod, alpha = 1, maxiter = 1000, checkDivergence = FALSE) {

# Initialization
nRows <- nrow(X)
nCols <- ncol(X)
iter <- 1
delta <- 0.01
relTol <- 1e-5 # Relative convergence tolerance
relError <- Inf
W <- NMF::basis(nmfMod)
H <- NMF::coef(nmfMod)
err <- rep(0, times = maxiter+1)
err[1] <- norm(X-W%*%H,'f')
frobX <- norm(X,'f')^2

while(iter < maxiter && relError > relTol) {
# Update of W
A <- X%*%t(H)
B <- H%*%t(H)
W <- HALSupdt(t(X),t(W),t(B),t(A),alpha,delta)
W <- t(W)

# Update of H
A <- t(W)%*%X
B <- t(W)%*%W
H <- HALSupdt(X,H,B,A,alpha,delta)

# Evaluation of convergence
err[iter+1] <- sqrt((frobX-2*sum(H*A))+sum(B*(H%*%t(H))))
relError <- abs(err[iter+1]-err[iter])/err[iter]

if(iter%%10 == 0) {
print(paste("Iteration:",as.character(iter),"Relative error:",as.character(relError)))
}
iter <- iter+1
}
NMF::basis(nmfMod) <- W
NMF::coef(nmfMod) <- H
#nmfMod <- nmfModel(W=W,H=H) # nmfModel function from the 'NMF' package
return(nmfMod)
}

## #' Updating step for accelerated HALS NMF
## #' @param M Input data matrix
## #' @param V Factor matrix to be updated
## #' @param UtU Product of the other transposed factor matrix with itself
## #' @param UtM Product of the other transposed factor matrix with the input matrix
## #' @param alpha Nonnegative parameter of the accelerated method
## #' @param delta Convergence parameter
## #' @return Updated factor matrix V
## #' @author Nicolas Sauwen
HALSupdt <- function (M,V,UtU,UtM,alpha,delta) {

nRows <- nrow(V)
nCols <- ncol(V)
cnt <- 1
eps <- 1
eps0 <- 1

while(eps > sqrt(delta)*eps0) {
noDelta <- 0
for(iRow in 1:nRows) {
deltaV <- pmax((UtM[iRow,]-UtU[iRow,]%*%V)/UtU[iRow,iRow],-V[iRow,])
V[iRow,] <- V[iRow,] + deltaV
noDelta <- noDelta + deltaV%*%t(deltaV)
zeroInds <- which(V[iRow,]==0)
V[iRow,zeroInds] <- 1e-16*max(V)
}
if(cnt == 1) {
eps0 <- noDelta
}
eps <- noDelta
cnt <- 0
}
return(V)
}
```

## Try the hNMF package in your browser

Any scripts or data that you put into this service are public.

hNMF documentation built on Jan. 8, 2021, 5:42 p.m.