djNMF)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)
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.
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:
## 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 \(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)
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 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.20.so; LAPACK version 3.10.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.2.0 fields_15.2 viridisLite_0.4.2 spam_2.9-1
## [5] dcTensor_1.3.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.4 jsonlite_1.8.7 highr_0.10 compiler_4.3.1
## [5] maps_3.4.1 Rcpp_1.0.11 plot3D_1.4 tagcloud_0.6
## [9] jquerylib_0.1.4 scales_1.2.1 yaml_2.3.7 fastmap_1.1.1
## [13] ggplot2_3.4.3 R6_2.5.1 tcltk_4.3.1 knitr_1.43
## [17] MASS_7.3-60 dotCall64_1.0-2 misc3d_0.9-1 tibble_3.2.1
## [21] munsell_0.5.0 pillar_1.9.0 bslib_0.5.1 RColorBrewer_1.1-3
## [25] rlang_1.1.1 utf8_1.2.3 cachem_1.0.8 xfun_0.40
## [29] sass_0.4.7 cli_3.6.1 magrittr_2.0.3 digest_0.6.33
## [33] grid_4.3.1 rTensor_1.4.8 lifecycle_1.0.3 vctrs_0.6.3
## [37] evaluate_0.21 glue_1.6.2 fansi_1.0.4 colorspace_2.1-0
## [41] rmarkdown_2.24 pkgconfig_2.0.3 tools_4.3.1 htmltools_0.5.6