<div dir="ltr">Hi there,<div><br></div><div>I am new to RcppArmadillo and very excited about the sparse matrix functionality it offers. I have a large optim that I wish to break down into stochastic descent in Rcpp to speed things up.</div>
<div><br></div><div>To start small, I used a simple logisitic regression function that I have working in R -</div><div><br></div><div><div>Norm <- function(w1, w2){</div><div>  sqrt(sum((w2 - w1)^2))</div><div>}</div><div>
<br></div><div>xentropy <- function(x, y, w) {</div><div>  - ((y * x) / (1 + exp(y * t(w) %*% x)) )</div><div>}</div><div><br></div><div>stochDesc <- function(param, eta = 0.01, maxit = 10, reltol = .01, x, y) {</div>
<div>  # Initialize</div><div>  w = param</div><div>  iter = 0</div><div>  rel_err = reltol + 1</div><div>  x_dash = cbind(rep(1, nrow(x)), x) # Add bias</div><div>  </div><div>  while(iter < maxit & rel_err > reltol ) { </div>
<div>    iter = iter + 1</div><div>    </div><div>    prev_w <- w</div><div>    for(curidx in sample(nrow(x))) { # Random shuffle</div><div>      grad <- xentropy(x_dash[curidx, ], y[curidx], w)</div><div>      w <- w - eta *  grad # Update weights</div>
<div>    }</div><div>    rel_err = Norm(w, prev_w)<br></div><div>    # could also use: </div><div>    # rel_err = base::norm(matrix(prev_w - w, nrow = 1), type = "F")</div><div>    cat(paste("Epoch ", iter, " Relative cost:", rel_err, "\n"))</div>
<div>  }</div><div>  list(w = w, iter = iter)</div><div>}</div></div><div><br></div><div>This does every operation that my more complex gradient descent will do, but it is pretty straightforward so I started converting it to Rcpp -</div>
<div><br></div><div><div>#include <RcppArmadillo.h></div><div>#include <iostream></div><div>using namespace std;</div><div>using namespace Rcpp ;</div><div>// [[Rcpp::depends(RcppArmadillo)]]</div><div><br></div>
<div>double Norm(NumericVector w1, NumericVector w2) { </div><div>  // Compute the Forbius Norm</div><div>  double op = pow(sum(pow(w2 - w1, 2)), .5); </div><div>  return(op);</div><div>}</div><div><br></div><div><br></div>
<div>NumericVector xentropy(NumericVector x, NumericVector y, NumericVector w) { </div><div>  // Return the cross entropy</div><div>  // NumericVector op = - ((y * x) / (1 + exp(y * trans(w) %*% x)) ); // How to transpose?</div>
<div>  NumericVector op = w; // Not correct</div><div>  return (op);</div><div>}</div><div><br></div><div><br></div><div>// [[Rcpp::export]]</div><div>NumericVector stochDescCpp(NumericVector param, NumericMatrix x, NumericVector y,</div>
<div>                          double eta = 0.01, int maxit = 10, double reltol = .01) {</div><div>  // Run stochastic Descent to compute logistic regression</div><div>  // Initialize</div><div>  NumericVector w (param);</div>
<div>  NumericVector prev_w (w);</div><div>  int epoch = 0;</div><div>  NumericVector grad (param.length());</div><div>  </div><div>  double rel_err = reltol + 1;</div><div>  // Don't know how to add bias, will do that in R before passing it here.</div>
<div>  //NumericMatrix x_dash(x.nrow(), x.ncol() + 1, rep(1, x.nrow()), x.begin());</div><div>  </div><div>  // Loop through for each epoch</div><div>  while(epoch < maxit && rel_err > reltol) { </div><div>    epoch++;</div>
<div>    </div><div>    prev_w = w;</div><div>    for(int curidx = 0; curidx < x.nrow(); ++curidx) { // How to random shuffle?</div><div>      grad = xentropy(x(curidx, _), y[curidx], w);</div><div>      w = w - eta * grad; // Update weights</div>
<div>    }</div><div>    rel_err = Norm(w, prev_w);</div><div>    cout << "Epoch " << epoch << " Relative cost:" << rel_err << endl;</div><div>  }</div><div>  return(w);</div>
<div>}</div></div><div><br></div><div>I am running into many issues. Here they are ranked in order of importance) -</div><div>1) When I run this (incorrect) program -</div><div>stochDescCpp(param = rep(0,3), x = X, y = Y, maxit = 100)<br>
</div><div><br></div><div>I get -</div><div><div>Error in .Primitive(".Call")(<pointer: 0x0000000069d42560>, param, x,  : </div><div>  negative length vectors are not allowed</div></div><div><br></div><div>
All the inputs are valid and work fine with the R version. So I don't know where this is coming from? I put a bunch of couts and I know that it is happening inside the for loop but don't know why. In general, what is the best way to debug code written this way? I am using RStudio.</div>
<div><br></div><div>2) I don't know how to transpose the matrix in rCpp. I read that I could that using rCppArmadillo, but then I don't know how to convert from NumericMatrix to arma matrix. I don't want to copy as this operation will happen 1000x of times, so this has to be fast.</div>
<div><br></div><div>3) I don't know how to randomly shuffle the rows. I tried <span class="" style="background-color:transparent;color:inherit;font-family:Menlo,Monaco,'Courier New',monospace;font-size:12.025px;line-height:18px;white-space:pre-wrap">RcppArmadillo</span><span class="" style="background-color:transparent;color:inherit;font-family:Menlo,Monaco,'Courier New',monospace;font-size:12.025px;line-height:18px;white-space:pre-wrap;font-weight:bold">::</span><span class="" style="background-color:transparent;color:inherit;font-family:Menlo,Monaco,'Courier New',monospace;font-size:12.025px;line-height:18px;white-space:pre-wrap">sample but haven't gotten the output to work. E.g. do I access elements as x(rowOrd[curIdx],_)?</span></div>
<div><span class="" style="background-color:transparent;color:inherit;font-family:Menlo,Monaco,'Courier New',monospace;font-size:12.025px;line-height:18px;white-space:pre-wrap"><br></span></div><div><span class="" style="background-color:transparent;color:inherit;font-family:Menlo,Monaco,'Courier New',monospace;font-size:12.025px;line-height:18px;white-space:pre-wrap">4) (minor issue) How do I add the bias unit? In other words, add a new column to the input matrix filled with all 1s. Do I have to create a new matrix, copy it column by column? this is not a dealbreaker since I can do this in R before I send in the inputs.</span></div>
<div><br></div><div>Sorry if my questions are too basic and don't seem well researched. I am trying to go for a practical example. </div><div><br></div><div>Any help will be greatly appreciated.</div><div><br></div><div>
Sincerely,</div><div>Saurabh</div></div>