Support Mixed BIC & Full-Grid Solves

Unify interface to full-grid and BIC-based selection to allow for "mixed" solves

Allow BIC-based selection to be set separately for each 
parameter (of 4) while other parameters have a full-grid of
solutions. This will be useful for special cases (e.g., two way FPCA
where we want BIC for smoothing parameters but to fix sparsity
parameters to zero) and for visualization.

R/moma_4Dlist_extractor.R

@@ -0,0 +1,25 @@
+get_4Dlist_elem <- function(x, alpha_u_i, lambda_u_i, alpha_v_i, lambda_v_i){
+    if(!inherits(x, "MoMA_4D_list")){
+        moma_error(sQuote("x"), " should be a ", sQuote("MoMA_4D_list"), " object.")
+    }
+    n_alpha_u = dim(x)[1]
+    n_lambda_u = dim(x)[2]
+    n_alpha_v = dim(x)[3]
+    n_lambda_v = dim(x)[4]
+
+    # NOTE: R index starts from 1
+    if(
+        alpha_u_i  <= 0 || alpha_u_i > n_alpha_u ||
+        lambda_u_i <= 0 || lambda_u_i > n_lambda_u ||
+        alpha_v_i  <= 0 || alpha_v_i > n_alpha_v ||
+        lambda_v_i <= 0 || lambda_v_i > n_lambda_v
+    ){
+        moma_error("Invalid index (",alpha_u_i, ",", lambda_u_i, 
+                   ",",alpha_v_i, ",",lambda_v_i,"), dim = ",
+                   dim(x))
+    }
+    return(x[n_lambda_u * n_alpha_v * n_lambda_v * (alpha_u_i-1) + 
+                          n_alpha_v * n_lambda_v * (lambda_u_i-1) + 
+                                      n_lambda_v * (alpha_v_i-1) +
+                                                   lambda_v_i])
+}

R/moma_svd.R

@@ -15,6 +15,18 @@ MOMA_DEFAULT_PROX <- list(
                         prox_eps = 1e-10,
                         # trend filtering
                         l1tf_k = 1)
+add_default_prox_args <- function(sparsity_type){
+    # sparsity_type: prox arguments for u and v
+
+    # To call a C function we have to specify
+    # all arguments. However, some arguments
+    # are specific for a particular prox. So
+    # we first assign a default arg list to 
+    # `df_prox_arg_list_u/_v` and
+    # then update them.
+    return(modifyList(MOMA_DEFAULT_PROX, sparsity_type))
+}
+
 
 # This function checks the validity of Omega and alpha
 check_omega <- function(Omega,alpha,n){
@@ -92,16 +104,6 @@ moma_svd <- function(
                         solver = toupper(solver),
                         k = k)
 
-    # Prox arguments for u and v
-    # To call a C function we have to specify
-    # all arguments. However, some arguments
-    # are specific for a particular prox. So
-    # we first assign a default arg list to 
-    # `df_prox_arg_list_u/_v` and
-    # then update them.
-    df_prox_arg_list_u <- MOMA_DEFAULT_PROX
-    df_prox_arg_list_v <- MOMA_DEFAULT_PROX
-
     if (!is.matrix(X)){
         moma_error("X must be a matrix.")
     }
@@ -139,8 +141,8 @@ moma_svd <- function(
                 list(
                     Omega_u = check_omega(Omega_u,alpha_u,n),
                     Omega_v = check_omega(Omega_v,alpha_v,p),
-                    prox_arg_list_u = modifyList(df_prox_arg_list_u,u_sparsity),
-                    prox_arg_list_v = modifyList(df_prox_arg_list_v,v_sparsity)))
+                    prox_arg_list_u = add_default_prox_args(u_sparsity),
+                    prox_arg_list_v = add_default_prox_args(v_sparsity)))
 
     if(is_multiple_para){
         if(select == "gridsearch"){

src/moma.cpp

@@ -47,6 +47,10 @@ MoMA::MoMA(const arma::mat &i_X, // Pass X_ as a reference to avoid copy
             i_EPS_inner,i_MAX_ITER_inner,i_X.n_cols)
      // const reference must be passed to initializer list
 {
+
+    bicsr_u.bind(&solver_u, &PR_solver::bic);
+    bicsr_v.bind(&solver_v, &PR_solver::bic);
+
     MoMALogger::info("Initializing MoMA object:")
     << " lambda_u " << lambda_u
     << " lambda_v " << lambda_v
@@ -214,3 +218,129 @@ int MoMA::reset(double newlambda_u,double newlambda_v,
     solver_v.reset(newlambda_v,newalpha_v);
     return 0;
 }
+
+const arma::vec &set_greedy_grid(const arma::vec &grid, int want_grid){
+    if (want_grid == 1){ 
+        return grid;
+    }
+    else if (want_grid == 0){
+        return MOMA_EMPTY_GRID_OF_LENGTH1;
+    }
+}
+
+arma::vec set_bic_grid(const arma::vec &grid, int want_bic, int i){
+    if(want_bic == 1){
+        return grid;
+    }
+    else if (want_bic == 0){
+        return grid(i) * arma::ones<arma::vec> (1);
+    }
+}
+
+Rcpp::List MoMA::grid_BIC_mix(const arma::vec &alpha_u,
+    const arma::vec &alpha_v,
+    const arma::vec &lambda_u,
+    const arma::vec &lambda_v,
+    int selection_criterion_alpha_u,  // flags; = 0 means grid, = 01 means BIC search
+    int selection_criterion_alpha_v,
+    int selection_criterion_lambda_u,
+    int selection_criterion_lambda_v,
+    int max_bic_iter){
+
+    // If alpha_u is selected via grid search, then the variable
+    // grid_au = alpha_u, bic_au_grid = [-1].
+    // If alpha_u is selected via nested BIC search,
+    // then grid_au = [-1], bic_au_grid = alpha_u
+    const arma::vec &grid_lu = set_greedy_grid(lambda_u, !selection_criterion_lambda_u);
+    const arma::vec &grid_lv = set_greedy_grid(lambda_v, !selection_criterion_lambda_v);
+    const arma::vec &grid_au = set_greedy_grid(alpha_u, !selection_criterion_alpha_u);
+    const arma::vec &grid_av = set_greedy_grid(alpha_v, !selection_criterion_alpha_v);
+
+    // Test that if a grid is set to be BIC-search grid, then
+    // the above code should set grid_xx to the vector [-1]
+    if((selection_criterion_alpha_u == 1 && (grid_au.n_elem != 1 || grid_au(0) != -1))
+        || (selection_criterion_alpha_v == 1 && (grid_av.n_elem != 1 || grid_av(0) != -1))
+        || (selection_criterion_lambda_u == 1 && (grid_lu.n_elem != 1 || grid_lu(0) != -1))
+        || (selection_criterion_lambda_v == 1 && (grid_lv.n_elem != 1 || grid_lv(0) != -1)) )
+    {
+        MoMALogger::error("Wrong grid-search grid!")
+        << "grid_lu.n_elem=" << grid_lu.n_elem
+        << ", grid_av.n_elem=" << grid_av.n_elem
+        << ", grid_lu.n_elem" << grid_lu.n_elem
+        << ", grid_lv.n_elem" << grid_lv.n_elem;
+    }
+
+    int n_lambda_u = grid_lu.n_elem;
+    int n_lambda_v = grid_lv.n_elem;
+    int n_alpha_u = grid_au.n_elem;
+    int n_alpha_v = grid_av.n_elem;
+
+
+    RcppFourDList four_d_list(n_alpha_u, n_lambda_u, n_alpha_v, n_lambda_v);
+
+    // nested-BIC search returns a list that
+    // contains (lambda, alpha, bic, selected vector)
+    Rcpp::List u_result;
+    Rcpp::List v_result;
+
+    // to facilitate warm-start
+    arma::vec oldu;
+    arma::vec oldv;
+
+    for(int i = 0; i < n_alpha_u; i++){
+        for(int j = 0; j < n_lambda_u; j++){
+            for(int k = 0; k < n_alpha_v; k++){
+                for(int m = 0; m < n_lambda_v; m++){
+
+                    arma::vec bic_au_grid = set_bic_grid(alpha_u, selection_criterion_alpha_u, i);
+                    arma::vec bic_lu_grid = set_bic_grid(lambda_u, selection_criterion_lambda_u, j);
+                    arma::vec bic_av_grid = set_bic_grid(alpha_v, selection_criterion_alpha_v, k);
+                    arma::vec bic_lv_grid = set_bic_grid(lambda_v, selection_criterion_lambda_v, m);
+
+                    if((selection_criterion_alpha_u == 0 && bic_au_grid.n_elem != 1)
+                        || (selection_criterion_alpha_v == 0 && bic_av_grid.n_elem != 1)
+                        || (selection_criterion_lambda_u == 0 && bic_lu_grid.n_elem != 1)
+                        || (selection_criterion_lambda_v == 0 && bic_lv_grid.n_elem != 1) )
+                    {
+                            MoMALogger::error("Wrong BIC search grid!");
+                    }
+
+                    tol = 1;
+                    iter = 0;
+
+                    // We conduct 2 BIC searches over 2D grids here instead 
+                    // of 4 searches over 1D grids. It's consistent with 
+                    // Genevera's code
+                    while(tol > EPS && iter < MAX_ITER && iter < max_bic_iter){
+                        iter++;
+                        oldu = u;
+                        oldv = v;
+
+                        // choose lambda/alpha_u
+                        MoMALogger::debug("Start u search.");
+                        u_result = bicsr_u.search(X*v, u, bic_au_grid, bic_lu_grid);
+                        u = Rcpp::as<Rcpp::NumericVector>(u_result["vector"]);
+
+                        MoMALogger::debug("Start v search.");
+                        v_result = bicsr_v.search(X.t()*u, v, bic_av_grid, bic_lv_grid);
+                        v = Rcpp::as<Rcpp::NumericVector>(v_result["vector"]);
+
+                        tol = norm(oldu - u) / norm(oldu) + norm(oldv - v) / norm(oldv);
+                        MoMALogger::message("Finish BIC search outer loop. (iter, tol) = (") 
+                                    << iter << "," << tol << "), "
+                                    << "(bic_u, bic_v) = (" 
+                                    << (double)u_result["bic"] << "," 
+                                    << (double)v_result["bic"] << ")";
+                    }
+
+                    four_d_list.insert(Rcpp::List::create(
+                                Rcpp::Named("u") = u_result,
+                                Rcpp::Named("v") = v_result), i, j, k, m);
+
+                }
+            }
+        }
+    }
+
+    return four_d_list.get_list();
+}

src/moma.h

@@ -10,6 +10,12 @@
 // Prox operators
 #include "moma_solver.h"
 
+// BIC searcher
+#include "moma_solver_BICsearch.h"
+
+// 4-D list
+#include "moma_fourdlist.h"
+
 // Prototypes
 // moma_logging.cpp
 void moma_set_logger_level_cpp(int);
@@ -28,6 +34,13 @@ class MoMA{
     double lambda_v;
 
 public:
+
+    // Receiver a grid of parameters
+    // and perform greedy BIC search
+    BIC_searcher bicsr_u;
+    BIC_searcher bicsr_v;
+
+
     // Our own copy of the data matrix
     // Modified when finding rank-k svd
     arma::mat X;
@@ -96,6 +109,17 @@ class MoMA{
     // change penalty level
     int reset(double newlambda_u,double newlambda_v,
                 double newalpha_u,double newalpha_v);
+
+    Rcpp::List grid_BIC_mix( 
+        const arma::vec &alpha_u,
+        const arma::vec &alpha_v,
+        const arma::vec &lambda_u,
+        const arma::vec &lambda_v,
+        int selection_criterion_alpha_u,  // flags; = 0 means grid, = 01 means BIC search
+        int selection_criterion_alpha_v,
+        int selection_criterion_lambda_u,
+        int selection_criterion_lambda_v,
+        int max_bic_iter=5);
 };
 
 #endif

src/moma_R_function.cpp

@@ -40,12 +40,12 @@ Rcpp::List cpp_sfpca(
               MAX_ITER_inner,
               solver);
 
-    int n_lu = lambda_u.n_elem;
-    int n_lv = lambda_v.n_elem;
-    int n_au = alpha_u.n_elem;
-    int n_av = alpha_v.n_elem;
+    int n_lambda_u = lambda_u.n_elem;
+    int n_lambda_v = lambda_v.n_elem;
+    int n_alpha_u = alpha_u.n_elem;
+    int n_alpha_v = alpha_v.n_elem;
 
-    int n_more_than_one = int(n_lv > 1) + int(n_lu > 1) + int(n_au > 1) + int(n_av > 1);
+    int n_more_than_one = int(n_lambda_v > 1) + int(n_lambda_u > 1) + int(n_alpha_u > 1) + int(n_alpha_v > 1);
     if(n_more_than_one > 0){
         MoMALogger::error("We don't allow a range of parameters in finding a rank-k svd.");
     }
@@ -95,21 +95,21 @@ Rcpp::List cpp_sfpca_grid(
     int k = 1){
 
     // We only allow changing two parameters
-    int n_lu = lambda_u.n_elem;
-    int n_lv = lambda_v.n_elem;
-    int n_au = alpha_u.n_elem;
-    int n_av = alpha_v.n_elem;
+    int n_lambda_u = lambda_u.n_elem;
+    int n_lambda_v = lambda_v.n_elem;
+    int n_alpha_u = alpha_u.n_elem;
+    int n_alpha_v = alpha_v.n_elem;
 
-    int n_more_than_one = int(n_lv > 1) + int(n_lu > 1) + int(n_au > 1) + int(n_av > 1);
+    int n_more_than_one = int(n_lambda_v > 1) + int(n_lambda_u > 1) + int(n_alpha_u > 1) + int(n_alpha_v > 1);
     if(n_more_than_one > 2){
         MoMALogger::error("We only allow changing two parameters.");
     }
 
-    if(n_lv == 0 || n_lu == 0 || n_au == 0 || n_av == 0){
+    if(n_lambda_v == 0 || n_lambda_u == 0 || n_alpha_u == 0 || n_alpha_v == 0){
         MoMALogger::error("Please specify all four parameters.");
     }
 
-    int n_total = n_lv * n_lu * n_au * n_av;
+    int n_total = n_lambda_v * n_lambda_u * n_alpha_u * n_alpha_v;
 
     // NOTE: arguments should be listed
     // in the exact order of MoMA constructor
@@ -137,10 +137,10 @@ Rcpp::List cpp_sfpca_grid(
     arma::vec d(n_total);
 
     int problem_id = 0;
-    for(int i = 0; i < n_lu; i++){
-        for(int j = 0; j < n_lv; j++){
-            for(int k = 0; k < n_au; k++){
-                for(int m = 0; m < n_av; m++){
+    for(int i = 0; i < n_lambda_u; i++){
+        for(int j = 0; j < n_lambda_v; j++){
+            for(int k = 0; k < n_alpha_u; k++){
+                for(int m = 0; m < n_alpha_v; m++){
                     MoMALogger::info("Setting up model:")
                                 << " lambda_u " << lambda_u(i)
                                 << " lambda_v " << lambda_v(j)
@@ -194,21 +194,21 @@ Rcpp::List cpp_sfpca_nestedBIC(
     int k = 1){
 
     // We only allow changing two parameters
-    int n_lu = lambda_u.n_elem;
-    int n_lv = lambda_v.n_elem;
-    int n_au = alpha_u.n_elem;
-    int n_av = alpha_v.n_elem;
+    int n_lambda_u = lambda_u.n_elem;
+    int n_lambda_v = lambda_v.n_elem;
+    int n_alpha_u = alpha_u.n_elem;
+    int n_alpha_v = alpha_v.n_elem;
 
-    int n_more_than_one = int(n_lv > 1) + int(n_lu > 1) + int(n_au > 1) + int(n_av > 1);
+    int n_more_than_one = int(n_lambda_v > 1) + int(n_lambda_u > 1) + int(n_alpha_u > 1) + int(n_alpha_v > 1);
     if(n_more_than_one > 2){
         MoMALogger::error("We only allow changing two parameters.");
     }
 
-    if(n_lv == 0 || n_lu == 0 || n_au == 0 || n_av == 0){
+    if(n_lambda_v == 0 || n_lambda_u == 0 || n_alpha_u == 0 || n_alpha_v == 0){
         MoMALogger::error("Please specify all four parameters.");
     }
 
-    int n_total = n_lv * n_lu * n_au * n_av;
+    int n_total = n_lambda_v * n_lambda_u * n_alpha_u * n_alpha_v;
 
     // NOTE: arguments should be listed
     // in the exact order of MoMA constructor

src/moma_base.h

@@ -17,5 +17,6 @@
 // we add a small "nugget" here to regularize the computations
 #define MOMA_EIGENVALUE_REGULARIZATION 0.01
 static constexpr double MOMA_INFTY = std::numeric_limits<double>::infinity();
+static const arma::vec MOMA_EMPTY_GRID_OF_LENGTH1 = - arma::ones<arma::vec> (1);
 #define MOMA_FUSEDLASSODP_BUFFERSIZE 5000
 #endif