Using the main function nmf()


  
# generate a synthetic dataset with known classes: 50 features, 23 samples (10+5+8)
n <- 20; counts <- c(5, 3, 2);
p <- sum(counts)
x <- syntheticNMF(n, counts)
dim(x)
[1] 20 10
# build the true cluster membership
groups <- unlist(mapply(rep, seq(counts), counts))
# run on a data.frame
res <- nmf(data.frame(x), 3)
# missing method: use algorithm suitable for seed
res <- nmf(x, 2, seed=rnmf(2, x))
algorithm(res)
[1] "brunet"
res <- nmf(x, 2, seed=rnmf(2, x, model='NMFns'))
algorithm(res)
[1] "nsNMF"
# compare some NMF algorithms (tracking the approximation error)
res <- nmf(x, 2, list('brunet', 'lee', 'nsNMF'), .options='t')
Compute NMF method 'brunet' [1/3] ... OK Compute NMF method 'lee' [2/3] ... OK Compute NMF method 'nsNMF' [3/3] ... OK
res
Length: 3 Method(s): brunet, lee, nsNMF Total timing: user system elapsed 1.142 0.000 1.142
summary(res, class=groups)
method seed rng metric rank sparseness.basis sparseness.coef purity entropy silhouette.coef silhouette.basis residuals niter brunet brunet random 1 KL 2 0.2932845 0.8254661 0.7 0.4984430 1 1 67.54475 420 lee lee random 1 euclidean 2 0.3725724 0.6479075 0.8 0.3063008 1 1 40.01078 450 nsNMF nsNMF random 1 KL 2 0.3613243 0.9999992 0.7 0.4984430 1 1 71.77176 430 cpu cpu.all nrun brunet 0.109 0.109 1 lee 0.164 0.164 1 nsNMF 0.207 0.207 1
# plot the track of the residual errors
plot(res)

# specify algorithm by its name
res <- nmf(x, 3, 'nsNMF', seed=123) # nonsmooth NMF
# names are partially matched so this also works
identical(res, nmf(x, 3, 'ns', seed=123))
[1] FALSE
res <- nmf(x, 3, 'offset') # NMF with offset
# run a custom algorithm defined as a standard function
myfun <- function(x, start, alpha){ # update starting point # ... basis(start) <- 3 * basis(start) # return updated point start }
res <- nmf(x, 2, myfun, alpha=3)
algorithm(res)
[1] "nmf_1eff752fc498"
# error: alpha missing
try( nmf(x, 2, myfun) )
# Model: features: 20 basis/rank: 2 samples: 10 # Details: algorithm: nmf_1eff6a26d7ae seed: random RNG: 403L, 111L, ..., 905346294L [32127b1f5d128e86611ca13100c4e072] distance metric: 'euclidean' residuals: 93757.77 Timing: user system elapsed 0.004 0.000 0.003
# possibly the algorithm fits a non-standard NMF model, e.g. NMFns model
res <- nmf(x, 2, myfun, alpha=3, model='NMFns')
modelname(res)
[1] "NMFns"
# assume a known NMF model compatible with the matrix `x`
y <- rnmf(3, x)
# fits an NMF model (with default method) on some data using y as a starting point
res <- nmf(x, y)
# the fit can be reproduced using the same starting point
nmf.equal(nmf(x, y), res)
[1] TRUE
# missing method: use default algorithm
res <- nmf(x, 3)
# Fit a 3-rank model providing an initial value for the basis matrix
nmf(x, rmatrix(nrow(x), 3), 'snmf/r')
# Model: features: 20 basis/rank: 3 samples: 10 # Details: algorithm: snmf/r seed: NMF RNG: 403L, 471L, ..., 905346294L [5b3bc938a22917243a2ad8c9fc79d83d] distance metric: residuals: 32.96566 Iterations: 70 Timing: user system elapsed 0.160 0.000 0.164
# Fit a 3-rank model providing an initial value for the mixture coefficient matrix
nmf(x, rmatrix(3, ncol(x)), 'snmf/l')
# Model: features: 20 basis/rank: 3 samples: 10 # Details: algorithm: snmf/l seed: NMF RNG: 403L, 515L, ..., 905346294L [a4cdbeabe17bc6669787b00b08588da7] distance metric: residuals: 32.93398 Iterations: 65 Timing: user system elapsed 0.151 0.000 0.151
# default fit
res <- nmf(x, 2)
summary(res, class=groups)
rank sparseness.basis sparseness.coef purity entropy silhouette.coef silhouette.basis residuals 2.0000000 0.3836735 0.6690929 0.8000000 0.3811979 1.0000000 1.0000000 66.2048250 niter cpu cpu.all nrun 450.0000000 0.1010000 0.1010000 1.0000000
# run default algorithm multiple times (only keep the best fit)
res <- nmf(x, 3, nrun=10)
# NOTE - CRAN check detected: limiting maximum number of cores [2/4] # NOTE - CRAN check detected: limiting maximum number of cores [2/4]
res
Method: brunet Runs: 10 RNG: 407L, 735887289L, 356183078L, -363657105L, -463838236L, -682061419L, -2098846382L Total timing: user system elapsed 2.854 0.228 3.019
summary(res, class=groups)
rank sparseness.basis sparseness.coef purity entropy silhouette.coef 3.0000000 0.4680864 0.7657260 0.9000000 0.1738140 0.9266244 silhouette.basis residuals niter cpu cpu.all nrun 0.8109707 44.6115320 420.0000000 0.1200000 2.8540000 10.0000000 cophenetic dispersion silhouette.consensus 0.9802005 0.7080000 0.7582481
# run default algorithm multiple times keeping all the fits
res <- nmf(x, 3, nrun=10, .options='k')
# NOTE - CRAN check detected: limiting maximum number of cores [2/4] # NOTE - CRAN check detected: limiting maximum number of cores [2/4]
res
Method: brunet Runs: 10 RNG: 407L, 1017738740L, -1532700571L, 258173410L, -924957701L, 1735077632L, 1899471617L Total timing: user system elapsed 4.230 0.372 2.651 Sequential timing: user system elapsed 1.323 0.001 1.324
summary(res, class=groups)
rank sparseness.basis sparseness.coef purity entropy silhouette.coef 3.0000000 0.4681643 0.7694107 0.9000000 0.1738140 0.9071356 silhouette.basis residuals niter cpu cpu.all nrun 0.7994680 44.6115286 430.0000000 0.1150000 4.2300000 10.0000000 cophenetic dispersion silhouette.consensus 0.9847183 0.7840000 0.8386278
## Note: one could have equivalently done
# res <- nmf(V, 3, nrun=10, .options=list(keep.all=TRUE))
# use a method that fit different model
res <- nmf(x, 2, 'nsNMF')
fit(res)
features: 20 basis/rank: 2 samples: 10 theta: 0.5
# pass parameter theta to the model via `...`
res <- nmf(x, 2, 'nsNMF', theta=0.2)
fit(res)
features: 20 basis/rank: 2 samples: 10 theta: 0.2
## handling arguments in `...` and model parameters
myfun <- function(x, start, theta=100){ cat("theta in myfun=", theta, "\n\n"); start }
# no conflict: default theta
fit( nmf(x, 2, myfun) )
theta in myfun= 100
features: 20 basis/rank: 2 samples: 10
# no conlfict: theta is passed to the algorithm
fit( nmf(x, 2, myfun, theta=1) )
theta in myfun= 1
features: 20 basis/rank: 2 samples: 10
# conflict: theta is used as model parameter
fit( nmf(x, 2, myfun, model='NMFns', theta=0.1) )
theta in myfun= 100
features: 20 basis/rank: 2 samples: 10 theta: 0.1
# conflict solved: can pass different theta to model and algorithm
fit( nmf(x, 2, myfun, model=list('NMFns', theta=0.1), theta=5) )
theta in myfun= 5
features: 20 basis/rank: 2 samples: 10 theta: 0.1
## USING SEEDING METHODS
# run default algorithm with the Non-negative Double SVD seeding method ('nndsvd')
res <- nmf(x, 3, seed='nndsvd')
## Note: partial match also works
identical(res, nmf(x, 3, seed='nn'))
[1] FALSE
# run nsNMF algorithm, fixing the seed of the random number generator
res <- nmf(x, 3, 'nsNMF', seed=123456)
nmf.equal(nmf(x, 3, 'nsNMF', seed=123456), res)
[1] TRUE
# run default algorithm specifying the starting point following the NMF standard model
start.std <- nmfModel(W=matrix(0.5, n, 3), H=matrix(0.2, 3, p))
nmf(x, start.std)
# Model: features: 20 basis/rank: 3 samples: 10 # Details: algorithm: brunet seed: NMF RNG: 403L, 342L, ..., 1047588911L [8ed89da992168bfe13a600ef97943baa] distance metric: 'KL' residuals: 89.23127 Iterations: 420 Timing: user system elapsed 0.092 0.000 0.092
# to run nsNMF algorithm with an explicit starting point, this one
# needs to follow the 'NMFns' model:
start.ns <- nmfModel(model='NMFns', W=matrix(0.5, n, 3), H=matrix(0.2, 3, p))
nmf(x, start.ns)
# Model: features: 20 basis/rank: 3 samples: 10 theta: 0.5 # Details: algorithm: nsNMF seed: NMF RNG: 403L, 342L, ..., 1047588911L [8ed89da992168bfe13a600ef97943baa] distance metric: 'KL' residuals: 70.70244 Iterations: 540 Timing: user system elapsed 0.245 0.000 0.244
# Note: the method name does not need to be specified as it is infered from the
# when there is only one algorithm defined for the model.
# if the model is not appropriate (as defined by the algorihtm) an error is thrown
# [cf. the standard model doesn't include a smoothing parameter used in nsNMF]
try( nmf(x, start.std, method='nsNMF') )
## Callback functions
# Pass a callback function to only save summary measure of each run
res <- nmf(x, 3, nrun=3, .callback=summary)
# NOTE - CRAN check detected: limiting maximum number of cores [2/4] # NOTE - CRAN check detected: limiting maximum number of cores [2/4]
# the callback results are simplified into a matrix
res$.callback
[,1] [,2] [,3] rank 3.0000000 3.0000000 3.0000000 sparseness.basis 0.4690227 0.4646504 0.4580022 sparseness.coef 0.8045101 0.7769446 0.8068001 silhouette.coef 0.9546248 0.8902991 0.8762803 silhouette.basis 0.8394310 0.7919504 0.7687868 residuals 44.6117396 44.6136628 45.3878071 niter 450.0000000 430.0000000 860.0000000 cpu 0.1150000 0.1080000 0.1870000 cpu.all 0.1150000 0.1080000 0.1870000 nrun 1.0000000 1.0000000 1.0000000
res <- nmf(x, 3, nrun=3, .callback=summary, .opt='-S')
# NOTE - CRAN check detected: limiting maximum number of cores [2/4] # NOTE - CRAN check detected: limiting maximum number of cores [2/4]
# the callback results are simplified into a matrix
res$.callback
[[1]] rank sparseness.basis sparseness.coef silhouette.coef silhouette.basis residuals niter cpu 3.0000000 0.3608160 0.8514485 0.9150132 0.7120422 50.4900252 480.0000000 0.1220000 cpu.all nrun 0.1220000 1.0000000 [[2]] rank sparseness.basis sparseness.coef silhouette.coef silhouette.basis residuals niter cpu 3.0000000 0.4684599 0.7880077 0.9424514 0.8147630 44.6115438 420.0000000 0.1080000 cpu.all nrun 0.1080000 1.0000000 [[3]] rank sparseness.basis sparseness.coef silhouette.coef silhouette.basis residuals niter cpu 3.0000000 0.4649104 0.8032040 0.9475480 0.8329278 44.6134353 440.0000000 0.1050000 cpu.all nrun 0.1050000 1.0000000
# Pass a custom callback function
cb <- function(obj, i){ if( i %% 2 ) sparseness(obj) >= 0.5 }
res <- nmf(x, 3, nrun=3, .callback=cb)
# NOTE - CRAN check detected: limiting maximum number of cores [2/4] # NOTE - CRAN check detected: limiting maximum number of cores [2/4]
res$.callback
[[1]] basis coef FALSE TRUE [[2]] NULL [[3]] basis coef FALSE TRUE
# Passs a callback function which throws an error
cb <- function(){ i<-0; function(object){ i <<- i+1; if( i == 1 ) stop('SOME BIG ERROR'); summary(object) }}
res <- nmf(x, 3, nrun=3, .callback=cb())
# NOTE - CRAN check detected: limiting maximum number of cores [2/4] # NOTE - CRAN check detected: limiting maximum number of cores [2/4] Warning message: NMF::nmf - All NMF fits were successful but 2/3 callback call(s) threw an error. # Callback error(s) thrown: - run #1: SOME BIG ERROR
## PARALLEL COMPUTATIONS
# try using 3 cores, but use sequential if not possible
res <- nmf(x, 3, nrun=3, .options='p3')
# NOTE - CRAN check detected: limiting maximum number of cores [2/4]
# force using 3 cores, error if not possible
res <- nmf(x, 3, nrun=3, .options='P3')
# NOTE - CRAN check detected: limiting maximum number of cores [2/4] Error: Parallel computation aborted: only 2 core(s) available [requested 3 core(s)]
# use externally defined cluster
library(parallel)
cl <- makeCluster(6)
res <- nmf(x, 3, nrun=3, .pbackend=cl)
# NOTE - CRAN check detected: limiting maximum number of cores [2/4] # NOTE - CRAN check detected: limiting maximum number of cores [2/4]
# use externally registered backend
registerDoParallel(cl)
res <- nmf(x, 3, nrun=3, .pbackend=NULL)