Merge pull request #117 from ailurophilia/main

avoid generating large matrix in macro_fisher_null.R

Author: adw96

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: radEmu
 Title: Using Relative Abundance Data to Estimate of Multiplicative Differences in Mean Absolute Abundance 
-Version: 1.2.0.0
+Version: 1.3.0.0
 Authors@R: as.person(c(
     "David Clausen <dsc24@uw.edu> [aut, cre]",
     "Amy Willis [aut]", 

R/fit_null.R

@@ -12,25 +12,25 @@
 #' @param rho_init where to start quadratic penalty parameter
 #' @param tau how much to increment rho by each iteration
 #' @param kappa cutoff above which to increment rho. If distance to feasibility doesn't shrink by at least this factor in an iteration, increment rho by tau.
-#' @param B_tol tolerance for convergence in max_{k,j} |B^t_{kj} - B^{(t - 1)}_{kj}|
+#' @param B_tol tolerance for convergence in $max_{k,j} |B^t_{kj} - B^{(t - 1)}_{kj}|$
 #' @param inner_tol tolerance for inner loop
-#' @param constraint_tol tolerance for |B_kj - g(B_k)|
+#' @param constraint_tol tolerance for $|B_kj - g(B_k)|$
 #' @param max_step maximum step size
 #' @param c1 constant for armijo rule
 #' @param maxit maximum iterations
 #' @param inner_maxit max iterations per inner loop
 #' @param verbose shout at you?
 #' @param trackB track value of beta across iterations and return?
 #' 
-#' @returns A list containing elements 'B', 'k_constr', 'j_constr', 'niter'
-#' 'gap', 'u', 'rho', and 'Bs'. 'B' is a matrix containing parameter estimates 
+#' @return A list containing elements `B`, `k_constr`, `j_constr`, `niter`
+#' `gap`, `u`, `rho`, and `Bs`. `B` is a matrix containing parameter estimates 
 #' under the null (obtained by maximum likelihood on augmented observations Y),
-#' 'k_constr', and 'j_constr' give row and column indexes of the parameter 
-#' fixed to be equal to the constraint function g() under the null. 'niter' is a 
+#' `k_constr`, and `j_constr` give row and column indexes of the parameter 
+#' fixed to be equal to the constraint function $g()$ under the null. `niter` is a 
 #' scalar giving total number of outer iterations used to fit the null model, 
-#' 'gap' gives the final value of g(B_{k_constr}) - B_{k_constr j_constr}, 
-#' 'u' and 'rho' are final values of augmented Lagrangian parameters, and 
-#' 'Bs' is a data frame containing values of B by iteration if trackB was set 
+#' `gap` gives the final value of $g(B_{k constr}) - B_{k constr, j constr}$, 
+#' `u` and `rho` are final values of augmented Lagrangian parameters, and 
+#' `Bs` is a data frame containing values of B by iteration if `trackB` was set 
 #' equal to TRUE (otherwise it contains a NULL value).
 #' 
 fit_null <- function(B,

R/macro_fisher_null.R

@@ -69,10 +69,8 @@ macro_fisher_null <- function(X,
   sm_denom <- as.numeric(as.matrix(1 + Matrix::crossprod(v,info_inverse)%*%v))
   sm_half_num <- info_inverse%*%v
 
-  #strictly speaking, this isn't the inverse info -- its the inverse of 
+  #strictly speaking, info_inverse isn't the inverse info -- its the inverse of 
   # (an approximation to) the hessian of the augmented lagrangian
-  info_inverse <- info_inverse - Matrix::tcrossprod(sm_half_num)/sm_denom
-  # the tcrossprod in the last line is the final step for the sherman-morrison update
 
   #compute derivative of augmented lagrangian
   lag_deriv <- lapply(1:J,
@@ -86,8 +84,8 @@ macro_fisher_null <- function(X,
 
 
   #direction for (approximate) Newton step
-  update_dir <- -info_inverse%*%lag_deriv
-
+  update_dir <- -info_inverse %*%lag_deriv + 
+    sm_half_num%*% Matrix::crossprod(sm_half_num,lag_deriv)/sm_denom
   #shorten step if any of its elements are larger than allowed by max_step
   if(max(abs(update_dir))>max_step){
     update_dir <- update_dir/max(abs(update_dir))
@@ -102,6 +100,12 @@ macro_fisher_null <- function(X,
                                    update_dir[,(1:(J - 1))>=j_ref,drop = FALSE]))
 
 
+  #in debug mode, compute full inverse info for backward compatibility
+  if(debug){
+    info_inverse <- info_inverse - Matrix::tcrossprod(sm_half_num)/sm_denom
+  }
+  
+  
 
 
 

R/score_test.R

@@ -10,7 +10,7 @@
 #' @param k_constr row index of element of B to be tested for equality to row identifiability constraint
 #' @param j_constr column index of element of B to be tested for equality to row identifiability constraint
 #' @param constraint_fn constraint function g to be imposed on each row of B; the null that
-#' is tested is B_{k_constr,j_constr} = g(B_k_constr)
+#' is tested is $B_{k_constr,j_constr} = g(B_k_constr)$
 #' @param constraint_grad_fn derivative of constraint_fn with respect to its
 #' arguments (i.e., elements of a row of B)
 #' @param j_ref column dropped from B in convenience parametrization used to compute 
@@ -70,19 +70,19 @@
 #' be used in computing GEE test statistics. Default is NULL, in which case rows of 
 #' Y are treated as independent.
 #'
-#' @return A list containing elements 'score_stat', 'pval', 'log_pval','niter',
-#' 'convergence', 'gap', 'u', 'rho', 'tau', 'inner_maxit', 'null_B', and 'Bs'. 'score_stat' gives the 
-#' value of the robust score statistic for H_0: B_{k_constr,j_constr} = g(B_{k_constr}).
-#' 'pval' and 'log_pval' are the p-value (on natural and log scales) corresponding to
+#' @return A list containing elements `score_stat`, `pval`, `log_pval`,'niter`,
+#' `convergence`, `gap`, `u`, `rho`, `tau`, `inner_maxit`, `null_B`, and `Bs`. `score_stat` gives the 
+#' value of the robust score statistic for $H_0: B_{k_constr,j_constr} = g(B_{k_constr})$.
+#' `pval` and `log_pval` are the p-value (on natural and log scales) corresponding to
 #' the score statistic (log_pval may be useful when the p-value is very close to zero). 
-#' 'gap' is the final value of g(B_{k_constr}) - B_{k_constr, j_constr} obtained in 
-#' optimization under the null. 'u' and 'rho' are final values of augmented 
-#' Lagrangian parameters returned by null fitting algorithm. 'tau' is the final value of 'tau' that 
-#' is used to update the 'rho' values and 'inner_maxit' is the final maximum number of iterations for 
+#' `gap` is the final value of $g(B_{k_constr}) - B_{k_constr, j_constr}$ obtained in 
+#' optimization under the null. `u` and `rho` are final values of augmented 
+#' Lagrangian parameters returned by null fitting algorithm. `tau` is the final value of `tau` that 
+#' is used to update the `rho` values and `inner_maxit` is the final maximum number of iterations for 
 #' the inner optimization loop in optimization under the null, in which B and z parameter values are
-#' maximized for specific 'u' and 'rho' parameters. 'null_B' is the value of 
-#' B returned but the null fitting algorithm. 'Bs' is by default NULL; if trackB = TRUE,
-#' 'Bs is a data frame containing values of B by outcome category, covariate, and 
+#' maximized for specific `u` and `rho` parameters. `null_B` is the value of 
+#' B returned but the null fitting algorithm. `Bs` is by default `NULL`; if `trackB = TRUE`,
+#' `Bs` is a data frame containing values of B by outcome category, covariate, and 
 #' iteration.
 #' 
 #' 

man/fit_null.Rd

@@ -352,3 +352,101 @@ test_that("We take similar step as we'd take using numerical derivatives in a mo
   expect_equal(update$update,n_update,tolerance = 1e-1)
 
 })
+
+test_that("We take same step as we'd take using numerical derivatives when gap, rho, u are zero", {
+  set.seed(59542234)
+  n <- 10
+  p <- 2
+  X <- cbind(1,rep(c(0,1),each = n/2))
+  J <- 10000
+  z <- rnorm(n) +8
+  b0 <- rnorm(J)
+  b1 <- seq(1,10,length.out = J)
+  b1 <- b1 - mean(b1)
+  b0 <- b0 - mean(b0)
+  b <- rbind(b0,b1)
+  Y <- matrix(NA,ncol = J, nrow = n)
+  
+  k_constr <- 2
+  j_constr <- 1
+  p <- 2
+  
+  constraint_fn <- function(x){ pseudohuber_center(x,0.1)}
+  # constraint_fn <- function(x){mean(x)}
+  
+  ##### Arguments to fix:
+  
+  constraint_grad_fn <- function(x){dpseudohuber_center_dx(x,0.1)}
+  # constraint_grad_fn <- function(x){ rep(1/length(x), length(x))}
+  rho_init = 1
+  tau = 1.2
+  kappa = 0.8
+  obj_tol = 100
+  score_tol <- 1e-3
+  constraint_tol = 1e-5
+  init_tol = 1e6
+  c1 = 1e-4
+  maxit = 1000
+  inner_maxit = 25
+  
+  Y[] <- 0
+  for(i in 1:n){
+    while(sum(Y[i,])==0){
+      for(j in 1:J){
+        temp_mean <- exp(X[i,,drop = FALSE]%*%b[,j,drop = FALSE] + z[i])
+        Y[i,j] <- rpois(1, lambda = temp_mean)
+        # Y[i,j] <- rnbinom(1,mu = temp_mean, size = 3)*rbinom(1,1,0.6)
+      }
+    }
+  }
+  
+  j_ref <- 5
+  
+  full_fit <- #suppressMessages(
+    emuFit_micro_penalized(X = X,
+                           Y = Y,
+                           B = NULL,
+                           constraint_fn = mean,
+                           tolerance = 10,
+                           verbose = FALSE)#)
+  
+  B <- full_fit$B
+  B[1,] <- B[1,] - B[1,j_ref]
+  B[2,] <- B[2,] - B[2,j_ref]
+  
+  Y_aug <- full_fit$Y_augmented
+  
+  B[k_constr,j_constr] <- constraint_fn(B[k_constr,-j_constr])
+  
+  u <- 0
+  rho <- 0
+  
+  z <- update_z(Y,X,B)
+  
+  macro_time <- try(system.time(macro_fisher_null(X = X,
+                              Y = Y_aug,
+                              B = B,
+                              z = z,
+                              J = J,
+                              p = p,
+                              k_constr = k_constr,
+                              j_constr = j_constr,
+                              j_ref = j_ref,
+                              rho = rho,
+                              u = u,
+                              constraint_fn = constraint_fn,
+                              constraint_grad_fn = constraint_grad_fn,
+                              stepsize = 0.5,
+                              c1 = 1e-4,
+                              regularization = 0,
+                              debug= FALSE)))
+
+  if(inherits(macro_time, "try-error")){
+    expect_true(FALSE)
+  }
+  expect_true(macro_time[3]<10)
+  #failing this test indicates that macro_fisher_null is probably directly 
+  #computing the full inverse of the (approximate) hessian matrix of the log likelihood
+  #this includes an outer product that does not need to be computed
+})
+