Merge pull request #138 from cvxgrp/better_memory_mgmt

Better memory mgmt

Author: bnaras

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: CVXR
 Type: Package
 Title: Disciplined Convex Optimization
-Version: 1.0-12
+Version: 1.0-13
 VignetteBuilder: knitr
 Authors@R: c(
            person("Anqi", "Fu",
@@ -63,7 +63,7 @@ LinkingTo: Rcpp, RcppEigen
 License: Apache License 2.0 | file LICENSE
 LazyData: true
 Collate:
-    'CVXR.R'
+    'CVXR-package.R'
     'data.R'
     'globals.R'
     'generics.R'
@@ -101,7 +101,7 @@ Collate:
     'CVXcanon-R6.R'
     'Deque.R'
     'canonInterface.R'
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 Encoding: UTF-8
 Enhances: 
     Rcplex,

NEWS.md

@@ -1,3 +1,10 @@
+# CVXR 1.0-13
+
+* Address inefficient processing of cones for MOSEK (Issue 137
+  reported by aszekMosek)
+* Fix `extract_quadratic_coeffs` to use sparse matrix and sweep in
+  place for better memory use (reported by Marissa Reitsma) 
+
 # CVXR 1.0-12
 
 * `Rmosek` to be removed from CRAN, so moved to drat repo

R/CVXR.R renamed to R/CVXR-package.R

@@ -20,12 +20,10 @@
 #' @import Matrix
 #' @importFrom stats rnorm runif
 #' @aliases CVXR-package CVXR
-#' @docType package
 #' @author Anqi Fu, Balasubramanian Narasimhan, John Miller, Steven Diamond, Stephen Boyd
-#'
-#' Maintainer: Anqi Fu<[email protected]>
-#' @keywords package
-NULL
+#' Maintainer: Anqi Fu<[email protected]>
+#' @keywords internal
+"_PACKAGE"
 
 
 

R/RcppExports.R

@@ -24,6 +24,14 @@
     .Call('_CVXR_cpp_convolve', PACKAGE = 'CVXR', xa, xb)
 }
 
+.sweep_dgCmat_in_place <- function(A, c_part) {
+    invisible(.Call('_CVXR_multiply_dgCMatrix_vector', PACKAGE = 'CVXR', A, c_part))
+}
+
+.sweep_in_place <- function(P, c_part) {
+    invisible(.Call('_CVXR_sweep_in_place', PACKAGE = 'CVXR', P, c_part))
+}
+
 #' Create a new LinOp object.
 #'
 #' @return an external ptr (Rcpp::XPtr) to a LinOp object instance.

R/coeff_extractor.R

@@ -126,14 +126,30 @@ setMethod("extract_quadratic_coeffs", "CoeffExtractor", function(object, affine_
       orig_id <- as.character(id(quad_forms[[var_id]][[3]]@args[[1]]))
       var_offset <- affine_id_map[[var_id]][[1]]
       var_size <- affine_id_map[[var_id]][[2]]
+      ind_seq <- seq.int(from = var_offset + 1, length.out = var_size)
       if(!any(is.na(value(quad_forms[[var_id]][[3]]@P)))) {
         P <- value(quad_forms[[var_id]][[3]]@P)
-        if(is(P, "sparseMatrix"))
-          P <- as.matrix(P)
-        c_part <- c[1, (var_offset + 1):(var_offset + var_size)]
-        P <- sweep(P, MARGIN = 2, FUN = "*", c_part)
-      } else
-        P <- sparseMatrix(i = 1:var_size, j = 1:var_size, x = c[1, (var_offset + 1):(var_offset + var_size)])
+        ## Assuming P is a matrix or a sparsematrix
+        ## Force P to be sparse
+        if (is.matrix(P)) {
+          ind <- which(P != 0, arr.ind = TRUE)
+          P <- sparseMatrix(p = c(0L, cumsum(tabulate(ind[, 2L], ncol(P)))),
+                            i = ind[, 1L] - 1L,
+                            x = P[ind], dims = dim(P), index1 = FALSE)
+        }
+        ## if(is(P, "sparseMatrix")) {          
+        ##  P <- as.matrix(P)
+        ##c_part <- c[1, (var_offset + 1):(var_offset + var_size)]
+        c_part <- c[1L, ind_seq]
+        ## P <- sweep(P, MARGIN = 2, FUN = "*", c_part)
+        ##if(is(P, "sparseMatrix"))
+        .Call('_CVXR_multiply_dgCMatrix_vector', PACKAGE = 'CVXR', P, c_part)
+        ##else
+        ##  .Call('_CVXR_sweep_in_place', PACKAGE = 'CVXR', P, c_part)          
+      } else {
+        ##P <- sparseMatrix(i = 1:var_size, j = 1:var_size, x = c[1, (var_offset + 1):(var_offset + var_size)])
+        P <- sparseMatrix(i = seq_len(var_size), j = seq_len(var_size), x = c[1L, ind_seq])
+      }
       if(orig_id %in% names(coeffs)) {
         coeffs[[orig_id]]$P <- coeffs[[orig_id]]$P + P
         coeffs[[orig_id]]$q <- coeffs[[orig_id]]$q + rep(0, nrow(P))
@@ -145,13 +161,16 @@ setMethod("extract_quadratic_coeffs", "CoeffExtractor", function(object, affine_
     } else {
       var_offset <- affine_id_map[[var_id]][[1]]
       var_size <- as.integer(prod(affine_var_dims[[var_id]]))
+      ind_seq <- seq.int(from = var_offset + 1, length.out = var_size)
       if(var_id %in% names(coeffs)) {
         coeffs[[var_id]]$P <- coeffs[[var_id]]$P + sparseMatrix(i = c(), j = c(), dims = c(var_size, var_size))
-        coeffs[[var_id]]$q <- coeffs[[var_id]]$q + c[1, (var_offset + 1):(var_offset + var_size)]
+        ##coeffs[[var_id]]$q <- coeffs[[var_id]]$q + c[1, (var_offset + 1):(var_offset + var_size)]
+        coeffs[[var_id]]$q <- coeffs[[var_id]]$q + c[1L, ind_seq]
       } else {
         coeffs[[var_id]] <- list()
         coeffs[[var_id]]$P <- sparseMatrix(i = c(), j = c(), dims = c(var_size, var_size))
-        coeffs[[var_id]]$q <- c[1, (var_offset + 1):(var_offset + var_size)]
+        ## coeffs[[var_id]]$q <- c[1, (var_offset + 1):(var_offset + var_size)]
+        coeffs[[var_id]]$q <- c[1L, ind_seq]
       }
     }
   }
@@ -470,4 +489,4 @@ restore_quad_forms <- function(expr, quad_forms) {
 #   }
 #   Ps <- lapply(Ps, function(P) { Matrix(P, sparse = TRUE) })
 #   list(Ps, Matrix(Q, sparse = TRUE), R)
-# }
+# }

R/conic_solvers.R

@@ -156,6 +156,14 @@ ConicSolver.get_spacing_matrix <- function(dim, spacing, offset) {
 ##   return(mat)
 ## }
 
+reshape_sparse_mat_by_row <- function(M, new_dim) {
+  ## Note, no check on new_dim for efficiency
+  m <- t(M)
+  dim(m) <- rev(new_dim)
+  t(m)
+}
+
+
 #' @param constr A \linkS4class{Constraint} to format.
 #' @param exp_cone_order A list indicating how the exponential cone arguments are ordered.
 #' @describeIn ConicSolver Return a list representing a cone program whose problem data tensors
@@ -191,11 +199,22 @@ setMethod("reduction_format_constr", "ConicSolver", function(object, problem, co
     #   coeffs[[2]][gap:(2*(gap-1)),]
     # TODO: Keep X_coeff sparse while reshaping!
     X_coeff <- coeffs[[2]]
-    reshaped <- matrix(t(X_coeff), nrow = nrow(coeffs[[1]]), byrow = TRUE)
-    stacked <- -cbind(coeffs[[1]], reshaped)
-    stacked <- matrix(t(stacked), nrow = nrow(coeffs[[1]]) + nrow(X_coeff), ncol = ncol(coeffs[[1]]), byrow = TRUE)
-
-    offset <- cbind(offsets[[1]], matrix(t(offsets[[2]]), nrow = nrow(offsets[[1]]), byrow = TRUE))
+    
+    ## reshaped <- matrix(t(X_coeff), nrow = nrow(coeffs[[1]]), byrow = TRUE)
+    cfs1 <- coeffs[[1]]
+    new_r <- nrow(cfs1)
+    new_c <- prod(dim(X_coeff)) / new_r
+    reshaped <- reshape_sparse_mat_by_row(X_coeff, new_dim = c(new_r, new_c))
+    
+    stacked <- -cbind(cfs1, reshaped)
+    ## stacked <- matrix(t(stacked), nrow = nrow(coeffs[[1]]) + nrow(X_coeff), ncol = ncol(coeffs[[1]]), byrow = TRUE)
+    stacked <- reshape_sparse_mat_by_row(
+      stacked,
+      new_dim = c(new_r + nrow(X_coeff), ncol(cfs1))
+    )
+
+    offset <- cbind(offsets[[1]],
+                    matrix(t(offsets[[2]]), nrow = nrow(offsets[[1]]), byrow = TRUE))
     offset <- matrix(t(offset), ncol = 1)
     return(list(Matrix(stacked, sparse = TRUE), as.vector(offset)))
   } else if(inherits(constr, "ExpCone")) {
@@ -2088,7 +2107,8 @@ setMethod("perform", signature(object = "MOSEK", problem = "Problem"), function(
   coeff_offs <- ConicSolver.get_coeff_offset(problem@objective@args[[1]])
   c <- coeff_offs[[1]]
   constant <- coeff_offs[[2]]
-  data[[C_KEY]] <- as.vector(c)
+  #data[[C_KEY]] <- as.vector(c)
+  data[[C_KEY]] <- c
   inv_data$n0 <- length(data[[C_KEY]])
   data[[OBJ_OFFSET]] <- constant[1]
   data[[DIMS]] <- list()

R/dcp2cone.R

@@ -58,14 +58,15 @@ setMethod("stuffed_objective", signature(object = "ConeMatrixStuffing", problem
   C <- CR[[1]]
   R <- CR[[2]]
 
-  c <- matrix(C, ncol = 1)   # TODO: Check if converted to dense matrix and flattened like in CVXPY
+  ##c <- matrix(C, ncol = 1)   # TODO: Check if converted to dense matrix and flattened like in CVXPY
+
   boolint <- extract_mip_idx(variables(problem))
   boolean <- boolint[[1]]
   integer <- boolint[[2]]
   # x <- Variable(extractor@N, boolean = boolean, integer = integer)
   x <- Variable(extractor@N, 1, boolean = boolean, integer = integer)
 
-  new_obj <- t(c) %*% x + 0
+  new_obj <- C %*% x + 0
 
   return(list(new_obj, x, R[1]))
 })

R/dgp2dcp.R

@@ -489,6 +489,8 @@ setMethod("names", signature(x = "DgpCanonMethods"), function(x) { names(Dgp2Dcp
 
 # TODO: How to implement this with S4 setMethod? Signature is x = "DgpCanonMethods", i = "character", j = "missing".
 '[[.DgpCanonMethods' <- function(x, i, j, ..., exact = TRUE) { do.call("$", list(x, i)) }
+## Here's an implementation using S4!
+## setMethod(f = "[[", signature = c(object = "DgpCanonMethods", i = "character", j = "missing") function(x, i, j) [email protected][[i]] )
 
 #' @param name The name of the atom or expression to canonicalize.
 #' @describeIn DgpCanonMethods Returns either a canonicalized variable or

R/generics.R

@@ -220,6 +220,12 @@ setGeneric("is_dcp", function(object) { standardGeneric("is_dcp") })
 #' @export
 setGeneric("is_dgp", function(object) { standardGeneric("is_dgp") })
 
+## #' DGP canonical methods indexing
+## #' @param A \linkS4class{DgpCanonMethods} object
+## #' @param a variable index of type character
+## #' @export
+## setGeneric("[[", function(object, i, j) standardGeneric("[["))
+
 #'
 #' Size of Expression
 #'

man/CVXR-package.Rd

@@ -1,4 +1,5 @@
 #include "CVXR.h"
+using namespace Rcpp;
 
 // Source:  Rcpp Introduction
 
@@ -14,3 +15,61 @@ Rcpp::NumericVector cpp_convolve(Rcpp::NumericVector xa, Rcpp::NumericVector xb)
       iab[i + j] += ia[i] * ib[j];
   return xab;
 }
+
+// [[Rcpp::export(.sweep_dgCmat_in_place)]]
+void multiply_dgCMatrix_vector(S4 A, NumericVector c_part) {
+  // Get the slot values from the dgCMatrix object
+  IntegerVector i = A.slot("i");
+  IntegerVector p = A.slot("p");
+  NumericVector x_values = A.slot("x");
+  
+  // Get the number of columns in the matrix
+  int n = p.length() - 1;
+  int k = c_part.length();
+  
+  // Check if the length of the vector matches the number of columns or for a scalar
+  if (k != n && k != 1) {
+    stop("mutiply_dgCMatrix_vector: Incompatible dimensions");
+  }
+
+  if (k == 1) { // c_part is a scalar
+    for (int i = 0; i < x_values.length(); ++i) {
+      x_values[i] *= c_part[0];
+    }
+  } else {
+    // Perform in-place multiplication
+    for (int col = 0; col < n; ++col) {
+      int start = p[col];
+      int end = p[col + 1];
+      
+      for (int idx = start; idx < end; ++idx) {
+	x_values[idx] *= c_part[col];
+      }
+    }
+  }
+  
+}
+
+// [[Rcpp::export(.sweep_in_place)]]
+void sweep_in_place(Rcpp::NumericMatrix P, Rcpp::NumericVector c_part) {
+  int nrow = P.nrow();
+  int ncol = P.ncol();
+  int k = c_part.length();
+  if (ncol != k && k != 1) {
+    Rcpp::stop("sweep_in_place: Incompatible dimensions");
+  }
+
+  if (k == 1) { // x is a scalar
+    for (int j = 0; j < ncol; j++) {
+      for (int i = 0; i < nrow; i++) {
+	P(i, j) = P(i, j) * c_part[0];
+      }
+    }
+  } else {
+    for (int j = 0; j < ncol; j++) {
+      for (int i = 0; i < nrow; i++) {
+	P(i, j) = P(i, j) * c_part[j];
+      }
+    }
+  }
+}

src/RcppConv.cpp

@@ -49,6 +49,28 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// multiply_dgCMatrix_vector
+void multiply_dgCMatrix_vector(S4 A, NumericVector c_part);
+RcppExport SEXP _CVXR_multiply_dgCMatrix_vector(SEXP ASEXP, SEXP c_partSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< S4 >::type A(ASEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type c_part(c_partSEXP);
+    multiply_dgCMatrix_vector(A, c_part);
+    return R_NilValue;
+END_RCPP
+}
+// sweep_in_place
+void sweep_in_place(Rcpp::NumericMatrix P, Rcpp::NumericVector c_part);
+RcppExport SEXP _CVXR_sweep_in_place(SEXP PSEXP, SEXP c_partSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type P(PSEXP);
+    Rcpp::traits::input_parameter< Rcpp::NumericVector >::type c_part(c_partSEXP);
+    sweep_in_place(P, c_part);
+    return R_NilValue;
+END_RCPP
+}
 // LinOp__new
 SEXP LinOp__new();
 RcppExport SEXP _CVXR_LinOp__new() {
@@ -415,6 +437,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_CVXR_build_matrix_0", (DL_FUNC) &_CVXR_build_matrix_0, 2},
     {"_CVXR_build_matrix_1", (DL_FUNC) &_CVXR_build_matrix_1, 3},
     {"_CVXR_cpp_convolve", (DL_FUNC) &_CVXR_cpp_convolve, 2},
+    {"_CVXR_multiply_dgCMatrix_vector", (DL_FUNC) &_CVXR_multiply_dgCMatrix_vector, 2},
+    {"_CVXR_sweep_in_place", (DL_FUNC) &_CVXR_sweep_in_place, 2},
     {"_CVXR_LinOp__new", (DL_FUNC) &_CVXR_LinOp__new, 0},
     {"_CVXR_LinOp__get_sparse", (DL_FUNC) &_CVXR_LinOp__get_sparse, 1},
     {"_CVXR_LinOp__set_sparse", (DL_FUNC) &_CVXR_LinOp__set_sparse, 2},

src/RcppExports.cpp

@@ -52,7 +52,7 @@ test_that("Test non-negative least squares", {
     coeff <- cbind(b, beta_ols, beta_nnls)
     colnames(coeff) <- c("Actual", "OLS", "NNLS")
     rownames(coeff) <- paste("beta", 1:length(b)-1, sep = "")
-    barplot(t(coeff), ylab = "Coefficients", beside = TRUE, legend = TRUE)
+    barplot(t(coeff), ylab = "Coefficients", beside = TRUE, legend.text = TRUE)
 })
 
 test_that("Test censored regression", {