Hierarchical Graph Factorization (Soft) Clustering. See: Yu, Kai, Shipeng Yu, and Volker Tresp. "Soft clustering on graphs." Advances in Neural Information Processing Systems 18 (2006): 1553.

hgfc.R

# Ozan Irsoy, 2013

# Graph factorization.
# W: adjacency matrix
# m: # clusters (# vertices in U)
gfact = function(W, m) {
  n = dim(W)[1]
  H = matrix(runif(n*m,1,2),n,m)
  L = diag(runif(m,1,2))
  
  for (iter in 1:100) {
    denom = H %*% L %*% t(H)
    W_ = W/denom
    W_[is.nan(W_)] = 1  # Let 0/0 = 1 because of divergence loss

    H_ = H*(W_ %*% H %*% L)
    H_ = scale(H_, center=F, scale=colSums(H_))
    L_ = L*(t(H)%*%W_%*%H)
    L_ = L_ / (sum(L_)/sum(W))
    H = H_
    L = L_
 
    #diver = sum(W*log(W_)-W+denom, na.rm=T) # let 0*-Inf = 0
    #print(diver)
  }
  return( list(B = H%*%L, L=L) )
}

# Hierarchical graph factorization clustering.
# W: adjacency matrix
# m: either the final # clusters, or an array of decreasing
#    # clusters for each level
# list.Bs: boolean. if true, will list B's for each level and 
#          return.
hgfc = function(W, m=NULL, list.Bs=FALSE) {
  n = dim(W)[1]
  if (is.null(m)) {m=(n-1):1}
  else if (length(m)==1) {m=(n-1):m} # this line needs to be modified

  Bs = list()
  
  for (l in 1:length(m)) {
    B = gfact(W, m[l])$B
    D = diag(rowSums(B))
    # TODO: sometimes D has a 0 entry in the diagonal, which means
    # there is an unused cluster. Fix!
    W = t(B)%*%solve(D)%*%B
    if (l==1) {B_ = B}
    else {B_ = B_%*%solve(D)%*%B}
    
    if (list.Bs) {Bs[[l]] = B}
  }
  
  if (list.Bs) {return (list(B=B_,Bs=Bs))}
  else {return(B_)}
}