4. Discretized Joint Non-negative Matrix Factrozation (djNMF)

Koki Tsuyuzaki

Laboratory for Bioinformatics Research, RIKEN Center for Biosystems Dynamics Research
k.t.the-answer@hotmail.co.jp

2023-03-17

Introduction

In this vignette, we consider approximating non-negative multiple matrices as a product of binary (or non-negative) low-rank matrices (a.k.a., factor matrices).

Test data is available from toyModel.

library("dcTensor")
suppressMessages(library("nnTensor"))
X <- nnTensor::toyModel("siNMF_Hard")

You will see that there are some blocks in the data matrices as follows.

suppressMessages(library("fields"))
layout(t(1:3))
image.plot(X[[1]], main="X1", legend.mar=8)
image.plot(X[[2]], main="X2", legend.mar=8)
image.plot(X[[3]], main="X3", legend.mar=8)

Semi-Binary Simultaneous Matrix Factorization (SBSMF)

Here, we consider the approximation of \(K\) binary data matrices \(X_{k}\) (\(N \times M_{k}\)) as the matrix product of \(W\) (\(N \times J\)) and \(V_{k}\) (J \(M_{k}\)):

\[ X_{k} \approx (W + V_{k}) H_{k} \ \mathrm{s.t.}\ W,V_{k},H_{k} \in \{0,1\} \]

This is the combination of binary matrix factorization (BMF (Zhang 2007)) and joint non-negative matrix decomposition (jNMF (Zi 2016; CICHOCK 2009)), which is implemented by adding binary regularization against jNMF. See also jNMF function of nnTensor package.

Basic Usage

In SBSMF, a rank parameter \(J\) (\(\leq \min(N, M)\)) is needed to be set in advance. Other settings such as the number of iterations (num.iter) or factorization algorithm (algorithm) are also available. For the details of arguments of djNMF, see ?djNMF. After the calculation, various objects are returned by djNMF. SBSMF is achieved by specifying the binary regularization parameter as a large value like the below:

set.seed(123456)
out_djNMF <- djNMF(X, Bin_W=1E-1, J=4)
str(out_djNMF, 2)
## List of 7
##  $ W            : num [1:100, 1:4] 0.343 0.338 0.346 0.344 0.342 ...
##  $ V            :List of 3
##   ..$ : num [1:100, 1:4] 2.02e-56 4.10e-56 2.26e-54 2.47e-55 7.54e-56 ...
##   ..$ : num [1:100, 1:4] 1.65e-63 2.33e-64 2.07e-60 2.49e-62 6.55e-61 ...
##   ..$ : num [1:100, 1:4] 0.156 0.143 0.157 0.155 0.15 ...
##  $ H            :List of 3
##   ..$ : num [1:300, 1:4] 4.16e-06 3.30e-06 3.38e-06 3.85e-06 7.51e-07 ...
##   ..$ : num [1:200, 1:4] 7.05e-20 7.47e-20 2.01e-20 4.32e-19 4.83e-20 ...
##   ..$ : num [1:150, 1:4] 95.3 95.9 96.4 94.1 94.9 ...
##  $ RecError     : Named num [1:101] 1.00e-09 1.14e+04 1.03e+04 9.94e+03 9.98e+03 ...
##   ..- attr(*, "names")= chr [1:101] "offset" "1" "2" "3" ...
##  $ TrainRecError: Named num [1:101] 1.00e-09 1.14e+04 1.03e+04 9.94e+03 9.98e+03 ...
##   ..- attr(*, "names")= chr [1:101] "offset" "1" "2" "3" ...
##  $ TestRecError : Named num [1:101] 1e-09 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 ...
##   ..- attr(*, "names")= chr [1:101] "offset" "1" "2" "3" ...
##  $ RelChange    : Named num [1:101] 1.00e-09 1.95e-01 1.12e-01 3.46e-02 3.95e-03 ...
##   ..- attr(*, "names")= chr [1:101] "offset" "1" "2" "3" ...

The reconstruction error (RecError) and relative error (RelChange, the amount of change from the reconstruction error in the previous step) can be used to diagnose whether the calculation is converged or not.

layout(t(1:2))
plot(log10(out_djNMF$RecError[-1]), type="b", main="Reconstruction Error")
plot(log10(out_djNMF$RelChange[-1]), type="b", main="Relative Change")

The products of \(W\) and \(H_{k}\)s show whether the original data matrices are well-recovered by djNMF.

recX1 <- lapply(seq_along(X), function(x){
  out_djNMF$W %*% t(out_djNMF$H[[x]])
})
recX2 <- lapply(seq_along(X), function(x){
  out_djNMF$V[[x]] %*% t(out_djNMF$H[[x]])
})
layout(rbind(1:3, 4:6, 7:9))
image.plot(X[[1]], legend.mar=8, main="X1")
image.plot(X[[2]], legend.mar=8, main="X2")
image.plot(X[[3]], legend.mar=8, main="X3")
image.plot(recX1[[1]], legend.mar=8, main="Reconstructed X1 (Common Factor)")
image.plot(recX1[[2]], legend.mar=8, main="Reconstructed X2 (Common Factor)")
image.plot(recX1[[3]], legend.mar=8, main="Reconstructed X3 (Common Factor)")
image.plot(recX2[[1]], legend.mar=8, main="Reconstructed X1 (Specific Factor)")
image.plot(recX2[[2]], legend.mar=8, main="Reconstructed X2 (Specific Factor)")
image.plot(recX2[[3]], legend.mar=8, main="Reconstructed X3 (Specific Factor)")

The histogram of \(W\) shows that the factor matrix \(W\) looks binary.

layout(rbind(1:4, 5:8))
hist(out_djNMF$W, main="W", breaks=100)
hist(out_djNMF$H[[1]], main="H1", breaks=100)
hist(out_djNMF$H[[2]], main="H2", breaks=100)
hist(out_djNMF$H[[3]], main="H3", breaks=100)
hist(out_djNMF$V[[1]], main="V1", breaks=100)
hist(out_djNMF$V[[2]], main="V2", breaks=100)
hist(out_djNMF$V[[3]], main="V3", breaks=100)

Session Information

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS/LAPACK: /home/koki/miniconda3/lib/libopenblasp-r0.3.17.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] nnTensor_1.1.12   fields_13.3       viridis_0.6.2     viridisLite_0.4.0
## [5] spam_2.8-0        dcTensor_1.0.1   
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.8         highr_0.9          RColorBrewer_1.1-2 rTensor_1.4.8     
##  [5] bslib_0.3.1        compiler_3.6.3     pillar_1.7.0       jquerylib_0.1.4   
##  [9] tools_3.6.3        dotCall64_1.0-1    digest_0.6.29      jsonlite_1.8.0    
## [13] evaluate_0.15      lifecycle_1.0.1    tibble_3.1.2       gtable_0.3.0      
## [17] pkgconfig_2.0.3    rlang_0.4.11       DBI_1.1.2          yaml_2.3.5        
## [21] xfun_0.29          fastmap_1.1.0      gridExtra_2.3      stringr_1.4.0     
## [25] dplyr_1.0.6        knitr_1.37         generics_0.1.2     sass_0.4.0        
## [29] vctrs_0.3.8        maps_3.4.0         plot3D_1.4         tidyselect_1.1.1  
## [33] grid_3.6.3         glue_1.4.2         R6_2.5.1           fansi_1.0.2       
## [37] tcltk_3.6.3        rmarkdown_2.11     purrr_0.3.4        ggplot2_3.3.5     
## [41] magrittr_2.0.2     scales_1.1.1       htmltools_0.5.2    ellipsis_0.3.2    
## [45] MASS_7.3-55        tagcloud_0.6       misc3d_0.9-1       assertthat_0.2.1  
## [49] colorspace_2.0-3   utf8_1.2.2         stringi_1.7.6      munsell_0.5.0     
## [53] crayon_1.5.0

References

CICHOCK, A. et al. 2009. Nonnegative Matrix and Tensor Factorizations. Wiley.
Zhang, Z. et al. 2007. “Binary Matrix Factorization with Applications.” ICDM 2007, 391–400.
Zi, et al., Yang. 2016. “A Non-Negative Matrix Factorization Method for Detecting Modules in Heterogeneous Omics Multi-Modal Data.” Bioinformatics 32(1): 1–8.