{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Multiomer Prediction with SVD + Regression","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle"}},{"cell_type":"markdown","source":"In this notebook, we present a solution to Multiomer prediction. More specifically, we will use \n\n1. We will use the sparse matrix to minimize RAM footprint from [this notebook](https://www.kaggle.com/code/stautxie/reduce-memory-footprint-of-competition-data) and [this dataset](https://www.kaggle.com/datasets/stautxie/sparse-measurement-data-open-problems-multimodal).\n\n2. Truncated Singular Value Decomposition to reduce the size of input matrix.\n\n3. We will explore lasso and ridge regression to predict the result.\n\n4. We develop a novel vectorization implementation of the competition metric.","metadata":{}},{"cell_type":"markdown","source":"While the overall idea is similar to [this wonderful notebook](https://www.kaggle.com/code/fabiencrom/msci-multiome-quickstart-w-sparse-matrices), the new contributions of this work are\n\n1. We also include result in both lasso and ridge regression while the previous work is only on ridge regression.  Our OOF study shows that lasso regression slightly outperforms ridge regression.\n2. We develop a new (and much straightforward) implementation of the scoring function without explicit looping.\n3. This is an R-implementation while all previous work is in Python.  While Scikit Learn has incorporated many relevant machine learning tools, an R solution typically needs to assemble various libraries. Hopefully, this work makes it easier for subsequent R-speaking Kagglers in this competition.","metadata":{}},{"cell_type":"code","source":"library(data.table)\nlibrary(Matrix)\nlibrary(magrittr)\nlibrary(tictoc)\nlibrary(ggplot2)\nlibrary(magrittr)\nlibrary(glmnet)\nlibrary(RSpectra)\nlibrary(caret)\nlibrary(ggplot2)","metadata":{"execution":{"iopub.status.busy":"2022-09-03T03:55:49.958835Z","iopub.execute_input":"2022-09-03T03:55:49.960815Z","iopub.status.idle":"2022-09-03T03:55:52.484149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DAT_DIR <- \"../input/open-problems-multimodal\"\nMAT_DIR <- \"../input/sparse-measurement-data-open-problems-multimodal\"","metadata":{"execution":{"iopub.status.busy":"2022-09-03T03:55:52.486423Z","iopub.execute_input":"2022-09-03T03:55:52.524870Z","iopub.status.idle":"2022-09-03T03:55:52.540286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let us first load the input matrix and use truncated SVD to reduce the dimension of input matrix so as to save both CPU and memory footprint.","metadata":{}},{"cell_type":"code","source":"tr_mu_inputs <- readRDS(file.path(MAT_DIR, \"sp_train_multi_inputs.rds\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-03T03:55:52.543034Z","iopub.execute_input":"2022-09-03T03:55:52.544661Z","iopub.status.idle":"2022-09-03T03:57:38.257038Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_comp <- 16\n\ntic()\nsvd_input <- svds(tr_mu_inputs,\n                  k = n_comp)\ntoc()\nsaveRDS(svd_input, \"mu_svd_input.rds\")","metadata":{"execution":{"iopub.status.busy":"2022-09-03T03:57:38.260000Z","iopub.execute_input":"2022-09-03T03:57:38.261576Z","iopub.status.idle":"2022-09-03T04:01:17.346229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data.table(d2 = svd_input$d[2:n_comp]^2 / max(svd_input$d[2:n_comp]^2),\n           i = seq_along(svd_input$d[2:n_comp])) %>%\n  ggplot(aes(x = i, y = d2)) +\n  geom_point() +\n  geom_line() +\n  labs(x = NULL, y = \"d squared / d[2] squared\") ","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:01:17.348770Z","iopub.execute_input":"2022-09-03T04:01:17.350125Z","iopub.status.idle":"2022-09-03T04:01:17.850636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see above, the first four singular values decrease abruptly and the trend becomes much slower afterwards.","metadata":{}},{"cell_type":"markdown","source":"To save memory footprint, we will randomly select 10% rows of the input and response matrices to save memory footprint.","metadata":{}},{"cell_type":"code","source":"set.seed(1234)\n\n#~ Take a sample of training data\nselected_pct <- 0.1\nselected_rows <- sample(1:nrow(tr_mu_inputs), floor(selected_pct * nrow(tr_mu_inputs)), replace=FALSE)\n\ntr_mu_input_svd <- svd_input$u %*% diag(svd_input$d, nrow=length(svd_input$d), ncol=length(svd_input$d))\ncolnames(tr_mu_input_svd) <- paste(\"SVD_\", 1:n_comp, sep=\"\")\n\ntr_mu_input_svd_sampled <- tr_mu_input_svd[selected_rows, ]\n","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:01:17.854244Z","iopub.execute_input":"2022-09-03T04:01:17.855709Z","iopub.status.idle":"2022-09-03T04:01:17.907470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(tr_mu_input_svd); gc()","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:01:17.913041Z","iopub.execute_input":"2022-09-03T04:01:17.916456Z","iopub.status.idle":"2022-09-03T04:01:18.293214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(tr_mu_inputs); gc()","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:01:18.295387Z","iopub.execute_input":"2022-09-03T04:01:18.296648Z","iopub.status.idle":"2022-09-03T04:01:19.063231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr_mu_targets <- readRDS(file.path(MAT_DIR, \"sp_train_multi_targets.rds\"))","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:01:19.066220Z","iopub.execute_input":"2022-09-03T04:01:19.067545Z","iopub.status.idle":"2022-09-03T04:02:20.784592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr_mu_targets_sampled <- tr_mu_targets[selected_rows,]\nrm(tr_mu_targets); gc()","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:02:20.788478Z","iopub.execute_input":"2022-09-03T04:02:20.790177Z","iopub.status.idle":"2022-09-03T04:02:30.279990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr_mu_targets_sampled <- as.matrix(tr_mu_targets_sampled)","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:02:30.282178Z","iopub.execute_input":"2022-09-03T04:02:30.283413Z","iopub.status.idle":"2022-09-03T04:02:33.855945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Scoring Function","metadata":{}},{"cell_type":"markdown","source":"The official scoring function can be found [here](https://www.kaggle.com/competitions/open-problems-multimodal/overview/evaluation) In summary, the predicted target is compared with the group solution for each sample (i.e., per row) and is measured by Pearson correlation coefficient. The final score is the mean of the row-wise correlation cofficient. ","metadata":{}},{"cell_type":"markdown","source":"A popular Python implementation of the scoring function can be found [here](https://www.kaggle.com/code/ambrosm/msci-citeseq-quickstart) where explicit loop is required to calculate the mean. However, we found a much more straightforward solution that uses vectorization instead of explicit looping.  This vectorization solution not only leads to more succinct representation but also potentially more efficient as well.","metadata":{}},{"cell_type":"markdown","source":"Our implementation of the scoring function is shown below in row_mean_corr(). This implementation takes advantage of\n\n1. Base R function cor() returns a matrix $A$ where each cell $a_{ij} := correlation(column_{i}, column_{j})$. By transposing the true target matrix $Y_{true}$ and prediction matrix $Y_{pred}$, the matrix $t(A)$ contains row-wise correlation.\n\n2. Base R function diag() returns diagonal values of a matrix.\n\n3. mean(diag(.)) returns the mean of row-wise correlation coefficient, which is the competition metric.","metadata":{}},{"cell_type":"code","source":"row_mean_corr <- function(Y_true, Y_pred) {\n  cor_mat <- cor(t(Y_true), t(Y_pred), method=\"pearson\")\n  return (mean(diag(cor_mat)))\n}","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:02:33.858771Z","iopub.execute_input":"2022-09-03T04:02:33.860167Z","iopub.status.idle":"2022-09-03T04:02:33.871301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will explore ridge and lasso regression below. To save time, we will use $\\alpha$=1 for both regressions.","metadata":{}},{"cell_type":"markdown","source":"### Ridge regression","metadata":{}},{"cell_type":"code","source":"set.seed(1234)\ntest_folds <- createFolds(1:nrow(tr_mu_input_svd_sampled), k=5, list=TRUE, returnTrain=FALSE)\n\ntest_fit_scores <- list()\ncor_fit_scores <- list()\n\nfor (fold in seq_along(test_folds)) {\n  tr_X <- tr_mu_input_svd_sampled[-test_folds[[fold]], ]\n  tr_Y <- tr_mu_targets_sampled[-test_folds[[fold]], ] %>% as.matrix()\n  te_X <- tr_mu_input_svd_sampled[test_folds[[fold]], ]\n  te_Y <- tr_mu_targets_sampled[test_folds[[fold]], ] %>% as.matrix()\n  fit <- glmnet(tr_X,\n                tr_Y,\n                family = \"mgaussian\",\n                alpha = 0,\n                nlambda = 1,\n                trace.it = TRUE)\n  \n  pred <- predict(fit, newx=te_X, s=1)\n  dim(pred) <- dim(te_Y)\n  mse <- mean((pred - te_Y)^2)\n  test_fit_scores[[fold]] <- mse\n  mean_cor_val <- row_mean_corr(te_Y, pred)\n  cor_fit_scores[[fold]] <- mean_cor_val\n  cat(paste(\"fold = \", fold, \", mse = \", mse, \", cor = \", mean_cor_val, \"\\n\"))\n  rm(tr_X, tr_Y, te_X, te_Y); gc()\n}","metadata":{"execution":{"iopub.status.busy":"2022-09-03T04:02:33.873534Z","iopub.execute_input":"2022-09-03T04:02:33.874760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(paste(\"OOF mean MSE: \", mean(as.numeric(test_fit_scores)), \"\\n\", sep=\"\"))\ncat(paste(\"OOF mean Cor: \", mean(as.numeric(cor_fit_scores)), \"\\n\", sep=\"\"))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Lasso Regression","metadata":{}},{"cell_type":"code","source":"test_folds <- createFolds(1:nrow(tr_mu_input_svd_sampled), k=5, list=TRUE, returnTrain=FALSE)\n\ntest_fit_scores <- list()\ncor_fit_scores <- list()\n\nfor (fold in seq_along(test_folds)) {\n  tr_X <- tr_mu_input_svd_sampled[-test_folds[[fold]], ]\n  tr_Y <- tr_mu_targets_sampled[-test_folds[[fold]], ] %>% as.matrix()\n  te_X <- tr_mu_input_svd_sampled[test_folds[[fold]], ]\n  te_Y <- tr_mu_targets_sampled[test_folds[[fold]], ] %>% as.matrix()\n  fit <- glmnet(tr_X,\n                tr_Y,\n                family = \"mgaussian\",\n                alpha = 1,\n                lambda = 1,\n                trace.it = TRUE)\n  \n  pred <- predict(fit, newx=te_X, s=1)\n  dim(pred) <- dim(te_Y)\n  mse <- mean((pred - te_Y)^2)\n  test_fit_scores[[fold]] <- mse\n  mean_cor_val <- row_mean_corr(te_Y, pred)\n  cor_fit_scores[[fold]] <- mean_cor_val\n  cat(paste(\"fold = \", fold, \", mse = \", mse, \", cor = \", mean_cor_val, \"\\n\"))\n  rm(tr_X, tr_Y, te_X, te_Y); gc()\n}","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(paste(\"OOF mean MSE: \", mean(as.numeric(test_fit_scores)), \"\\n\", sep=\"\"))\ncat(paste(\"OOF mean Cor: \", mean(as.numeric(cor_fit_scores)), \"\\n\", sep=\"\"))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Fit Prediction Model","metadata":{}},{"cell_type":"markdown","source":"Since lasso regression produces slightly better model on OOF, we will fit a lasso regression model with all the sampled data.","metadata":{}},{"cell_type":"code","source":"fit <- glmnet(tr_mu_input_svd_sampled,\n              tr_mu_targets_sampled %>% as.matrix(),\n              family = \"mgaussian\",\n              alpha = 1,\n              lambda = 1,\n              trace.it = TRUE)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"saveRDS(fit, \"fit_mu_lasso.rds\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(tr_mu_input_svd_sampled, tr_mu_targets_sampled); gc()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Predicting on Test data","metadata":{}},{"cell_type":"markdown","source":"We load the testing data set, predict with the fitted lasso model, and save the results for subsequent use.","metadata":{}},{"cell_type":"code","source":"te_mu_inputs <- readRDS(file.path(MAT_DIR, \"sp_test_multi_inputs.rds\"))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will use the fitted SVD transformation matrix to convert testing input as well.","metadata":{}},{"cell_type":"code","source":"svd_input <- readRDS(\"mu_svd_input.rds\")\nte_mu_inputs_svd <- te_mu_inputs %*% svd_input$v\nrm(te_mu_inputs); gc()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We load the saved model to predict on the testing sample.","metadata":{}},{"cell_type":"code","source":"fit <- readRDS(\"fit_mu_lasso.rds\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Predicting on the complete test data set will exceed 16GB memory limit, we will just predict on the first 5000 rows. It is possible to write an iterative routine to circumvent the memory limit.","metadata":{}},{"cell_type":"code","source":"te_pred <- predict(fit, newx=te_mu_inputs_svd[1:5000,], s=1)\ndim(te_pred) <- dim(te_pred)[-3]\nsaveRDS(te_pred, \"mu_te_pred.rds\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}