# what went wrong with my logistic regression implementation in c++?

I have implemented a simple logistic regression function with IRLS algorithm using the armadillo linear algebra libray

#include <iostream>
#include <string>
#include <boost/math/distributions/normal.hpp>
#include <boost/math/distributions/students_t.hpp>
#include <cmath>

using namespace boost::math;
arma::mat getW(
arma::mat& beta,
arma::mat& X,
std::string family,
)
{
arma::mat w;
if(family == "poisson") {
w = arma::diagmat(1/(X * beta));
}
}
else if(family == "binomial") {
arma::colvec tmp = exp(X * beta);
w = arma::diagmat(tmp/pow(1+tmp, 2));
}
}
else {
throw 1;
}
return w;
}

arma::mat getz(
arma::mat& y,
arma::mat& beta,
arma::mat& X,
std::string family,
)
{
arma::mat z;
if(family == "poisson") {
z = y;
}
}
else if(family == "binomial") {
arma::mat tmp = exp(X * beta);
z = X*beta + y % (pow(1+tmp, 2)/tmp) - 1 - tmp;
}
}
else {
throw 1;
}
return z;
}

inline arma::mat glmMat(
arma::mat& y,
arma::mat& x,
std::string family,
)
{

int n = x.n_rows;
int k = x.n_cols;

// add a col of all ones
arma::mat allOne(n, 1, arma::fill::ones);
x.insert_cols(0, allOne);
++k;

arma::mat res(k, 4);

{
arma::mat coef(k, 1, arma::fill::zeros);
arma::mat W = getW(coef, x, family, link);
arma::mat z = getz(y, coef, x, family, link);

try {
arma::mat J = x.t() * W * x;
arma::colvec coef1 = arma::solve(J, x.t()*W*z);
double coefdiff = max(abs(coef - coef1));
while(coefdiff >= 0.00001) {
coef = coef1;
W = getW(coef, x, family, link);
z = getz(y, coef, x, family, link);
J = x.t() * W * x;
coef1 = arma::solve(J, x.t()*W*z);
coefdiff = max(abs(coef - coef1));
}

arma::mat coefVarMatrix = J.i();
arma::colvec coefVar = coefVarMatrix.diag();
arma::colvec coefSe = pow(coefVar, .5);
arma::colvec zscore = coef / coefSe;
res.col(0) = coef;
res.col(1) = coefSe;
res.col(2) = zscore;

// calculate p values
auto d = normal_distribution<>();
for(int i=0; i<k; i++) {
double p = 2 * (1 - cdf(d, fabs(res(i, 2))));
if(p < 0 or p > 1) {
std::cerr << "Pval is abnormal from glm, dumping data to /tmp/tmpx.csv and /tmp/tmpy.csv" << std::endl;
x.save("/tmp/tmpx.csv", arma::csv_ascii);
y.save("/tmp/tmpy.csv", arma::csv_ascii);
throw 1;
}
res(i, 3) = p;
}
}
catch(...) {
std::cout << "something wrong..." << std::endl;
}
}
else {
throw 1;
}
return res;
}

int main(int argc, char const* argv[])
{
{
int nr = 5000;
int ncx = 50;
arma::mat x(nr, ncx, arma::fill::randu);
arma::mat y = arma::randi<arma::mat>(nr, 1, arma::distr_param(0, 1));
arma::mat xcol;
arma::mat res(ncx, 4);
for(int i=0; i<ncx; i++) {
xcol = x(arma::span::all, i);
res.row(i) = (glmMat(y, xcol, "binomial", "logit")).row(1);
}
res.print("res..........");
}
return 0;
}


Compiled it like this:

g++ glm.cpp --std=c++11 -larmadillo -llapack -lblas -o bin


The main function simulates a 5000x50 data set and performs logistic regression on each of them, the whole process takes about 23 seconds on my laptop.

Doing largely the same thing in R, it takes about 2 seconds:

testglm = function() {
x = matrix(rnorm(5000*50), 5000)
y = matrix(sample(0:1, 5000, repl=T), 5000)
res = apply(x, 2, function(coli) summary(glm(y~coli, family=binomial))\$coef[2, ])
# print(res)
}

system.time(testglm())
user  system elapsed
2.049   0.000   2.049


I am wondering what went wrong with my implementation?

• Using conditional statements in functions that are called repeatedly is not efficient. Is it possible to move those if statements outside the while loop? – Biswajit Banerjee Apr 3 '14 at 23:48

You are compiling without optimizations enabled:

g++ glm.cpp --std=c++11 -larmadillo -llapack -lblas -o bin

g++ -O3 -DNDEBUG glm.cpp --std=c++11 -larmadillo -llapack -lblas -o bin

Also, consider using Clang as an alternative compiler to g++ and see which produces faster code. It accepts the same command-line flags as g++ so it should be a drop-in replacement.

Disclaimer: I have worked on LLVM (which Clang uses internally).

A secondary and minor point is that you are passing around std::string by value, which causes many copies to be made (which requires memory allocation, data copies, and memory deallocation):

arma::mat getz(
arma::mat& y,
arma::mat& beta,
arma::mat& X,
std::string family,
)


Instead, pass non-modified data by const reference to make it more efficient:

arma::mat getz(
arma::mat& y,
arma::mat& beta,
arma::mat& X,
const std::string& family,

and similarly for the glmMat() function.