updated initialization to have smaller means and updated glmpca scripts

stephenslab · Jan 3, 2024 · 53a9205 · 53a9205
1 parent 6b305fe
commit 53a9205
Show file tree

Hide file tree

Showing 3 changed files with 48 additions and 132 deletions.
diff --git a/R/init.R b/R/init.R
@@ -112,8 +112,8 @@ init_glmpca_pois <- function(
            "of Y.")
 
     # Initialize U and V.
-    U <- matrix(rnorm(n*K,sd = 0.1),n,K)
-    V <- matrix(rnorm(m*K,sd = 0.1),m,K)
+    U <- matrix(rnorm(n*K,sd = 1e-4),n,K)
+    V <- matrix(rnorm(m*K,sd = 1e-4),m,K)
   }
 
   # Check and prepare input arguments X and B.

diff --git a/inst/experiments/fit_droplets_glmpca.R b/inst/experiments/fit_droplets_glmpca.R
@@ -1,80 +1,38 @@
 command_args = commandArgs(trailingOnly = TRUE)
 n_factor = as.integer(command_args[1])
 n_iter = as.integer(command_args[2])
-optimizer = as.character(command_args[3])
-minibatch = as.character(command_args[4])
 
 load("/project2/mstephens/pcarbo/git/fastTopics-experiments/data/droplet.RData")
 
-if (optimizer == "fisher") {
-
-  data <- as.matrix(counts)
-
-}
-
 set.seed(1)
-if (optimizer == "fisher") {
-
-  fit0 <- fastglmpca::init_glmpca(
-    Y = t(data), 
-    K = n_factor, 
-    fit_col_size_factor = TRUE,
-    fit_row_intercept = TRUE
-  )
-
-} else if (optimizer == "avagrad"){
-
-  fit0 <- fastglmpca::init_glmpca_pois(
-    Y = Matrix::t(counts), 
-    K = n_factor, 
-    fit_col_size_factor = TRUE,
-    fit_row_intercept = TRUE
-  )
-
-}
 
+fit0 <- fastglmpca::init_glmpca_pois(
+  Y = Matrix::t(counts),
+  K = n_factor
+)
 
 library(tictoc)
 
-if(optimizer == "fisher") {
-
-  fit <- glmpca::glmpca(
-    Y = t(data),
-    L = n_factor,
-    fam = "poi",
-    optimizer = optimizer,
-    minibatch = minibatch,
-    ctl = list(
-      minIter = n_iter - 1, 
-      maxIter = n_iter, 
-      verbose = TRUE, 
-      tol = .Machine$double.eps,
-      lr = 1e-4
-    ),
-    init = list(factors = fit0$V, loadings = fit0$U)
-  )
-
-} else {
-
-  fit <- glmpca::glmpca(
-    Y = Matrix::t(counts),
-    L = n_factor,
-    fam = "poi",
-    optimizer = optimizer,
-    minibatch = minibatch,
-    ctl = list(
-      minIter = n_iter - 1, 
-      maxIter = n_iter, 
-      verbose = TRUE, 
-      tol = .Machine$double.eps,
-      lr = 1e-4
-    ),
-    init = list(factors = fit0$V, loadings = fit0$U)
+fit <- glmpca::glmpca(
+  Y = Matrix::t(counts),
+  L = n_factor,
+  fam = "poi",
+  optimizer = optimizer,
+  minibatch = minibatch,
+  ctl = list(
+    minIter = n_iter - 1,
+    maxIter = n_iter,
+    verbose = TRUE,
+    tol = .Machine$double.eps,
+    lr = 1e-4
+  ),
+  init = list(
+    factors = fit0$V %*% diag(sqrt(fit0$d)),
+    loadings = fit0$U %*% diag(sqrt(fit0$d))
   )
-
-}
+)
 
 readr::write_rds(
   fit,
-  glue::glue("droplets_glmpca_fit_{n_factor}_factors_{n_iter}_iter_{optimizer}_optimizer_minibatch_{minibatch}.rds")
-)
+  glue::glue("droplets_glmpca_fit_{n_factor}_factors_{n_iter}_iter_avagrad_optimizer_minibatch_stochastic.rds")
+)
diff --git a/inst/experiments/fit_pbmc_glmpca.R b/inst/experiments/fit_pbmc_glmpca.R
@@ -1,80 +1,38 @@
 command_args = commandArgs(trailingOnly = TRUE)
 n_factor = as.integer(command_args[1])
 n_iter = as.integer(command_args[2])
-optimizer = as.character(command_args[3])
-minibatch = as.character(command_args[4])
 
 load("/project2/mstephens/pcarbo/git/fastTopics-experiments/data/pbmc_68k.RData")
 
-if (optimizer == "fisher") {
-
-  data <- as.matrix(counts)
-
-}
-
 set.seed(1)
-if (optimizer == "fisher") {
-
-  fit0 <- fastglmpca::init_glmpca(
-    Y = t(data), 
-    K = n_factor, 
-    fit_col_size_factor = TRUE,
-    fit_row_intercept = TRUE
-  )
-
-} else if (optimizer == "avagrad"){
-
-  fit0 <- fastglmpca::init_glmpca_pois(
-    Y = Matrix::t(counts), 
-    K = n_factor, 
-    fit_col_size_factor = TRUE,
-    fit_row_intercept = TRUE
-  )
-
-}
 
+fit0 <- fastglmpca::init_glmpca_pois(
+  Y = Matrix::t(counts),
+  K = n_factor
+)
 
 library(tictoc)
 
-if(optimizer == "fisher") {
-
-  fit <- glmpca::glmpca(
-    Y = t(data),
-    L = n_factor,
-    fam = "poi",
-    optimizer = optimizer,
-    minibatch = minibatch,
-    ctl = list(
-      minIter = n_iter - 1, 
-      maxIter = n_iter, 
-      verbose = TRUE, 
-      tol = .Machine$double.eps,
-      lr = 1e-4
-    ),
-    init = list(factors = fit0$V, loadings = fit0$U)
-  )
-
-} else {
-
-  fit <- glmpca::glmpca(
-    Y = Matrix::t(counts),
-    L = n_factor,
-    fam = "poi",
-    optimizer = optimizer,
-    minibatch = minibatch,
-    ctl = list(
-      minIter = n_iter - 1, 
-      maxIter = n_iter, 
-      verbose = TRUE, 
-      tol = .Machine$double.eps,
-      lr = 1e-4
-    ),
-    init = list(factors = fit0$V, loadings = fit0$U)
+fit <- glmpca::glmpca(
+  Y = Matrix::t(counts),
+  L = n_factor,
+  fam = "poi",
+  optimizer = optimizer,
+  minibatch = minibatch,
+  ctl = list(
+    minIter = n_iter - 1,
+    maxIter = n_iter,
+    verbose = TRUE,
+    tol = .Machine$double.eps,
+    lr = 1e-4
+  ),
+  init = list(
+    factors = fit0$V %*% diag(sqrt(fit0$d)),
+    loadings = fit0$U %*% diag(sqrt(fit0$d))
   )
-
-}
+)
 
 readr::write_rds(
   fit,
-  glue::glue("pbmc_glmpca_fit_{n_factor}_factors_{n_iter}_iter_{optimizer}_optimizer_minibatch_{minibatch}.rds")
-)
+  glue::glue("pbmc_glmpca_fit_{n_factor}_factors_{n_iter}_iter_avagrad_optimizer_minibatch_stochastic.rds")
+)