5. Discretized Joint Non-negative Matrix
Factrozation (djNMF)
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.
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 Xk (N × Mk) as the matrix product of W (N × J) and Vk (J Mk):
Xk ≈ (W + Vk)Hk s.t. W, Vk, Hk ∈ {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 ( ≤ 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:
## 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.04e-56 4.12e-56 2.27e-54 2.49e-55 7.58e-56 ...
## ..$ : num [1:100, 1:4] 1.65e-63 2.34e-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.17e-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.33e-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.96e-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 Hks 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 4.4.3 (2025-02-28)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [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
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] nnTensor_1.3.0 fields_16.3.1 viridisLite_0.4.2 spam_2.11-1
## [5] dcTensor_1.3.0 rmarkdown_2.29
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 jsonlite_1.9.1 compiler_4.4.3 maps_3.4.2.1
## [5] Rcpp_1.0.14 plot3D_1.4.1 tagcloud_0.6 jquerylib_0.1.4
## [9] scales_1.3.0 yaml_2.3.10 fastmap_1.2.0 ggplot2_3.5.1
## [13] R6_2.6.1 tcltk_4.4.3 knitr_1.49 MASS_7.3-65
## [17] dotCall64_1.2 misc3d_0.9-1 tibble_3.2.1 maketools_1.3.2
## [21] munsell_0.5.1 pillar_1.10.1 bslib_0.9.0 RColorBrewer_1.1-3
## [25] rlang_1.1.5 cachem_1.1.0 xfun_0.51 sass_0.4.9
## [29] sys_3.4.3 cli_3.6.4 magrittr_2.0.3 digest_0.6.37
## [33] grid_4.4.3 rTensor_1.4.8 lifecycle_1.0.4 vctrs_0.6.5
## [37] evaluate_1.0.3 glue_1.8.0 buildtools_1.0.0 colorspace_2.1-1
## [41] pkgconfig_2.0.3 tools_4.4.3 htmltools_0.5.8.1