1
$\begingroup$

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 <armadillo>
#include <cmath>

using namespace boost::math;
arma::mat getW(
    arma::mat& beta,
    arma::mat& X,
    std::string family,
    std::string link
)
{
    arma::mat w;
    if(family == "poisson") {
        if(link == "identity") {
            w = arma::diagmat(1/(X * beta));
        }
    }
    else if(family == "binomial") {
        if(link == "logit") {
            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,
    std::string link
)
{
    arma::mat z;
    if(family == "poisson") {
        if(link == "identity") {
            z = y;
        }
    }
    else if(family == "binomial") {
        if(link=="logit") {
            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,
    std::string link
)
{

    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);

    if(family=="binomial" and link=="logit")
    {
        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?

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

You are compiling without optimizations enabled:

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

Instead, try using:

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,
    std::string link
)

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,
    const std::string& link
)

and similarly for the glmMat() function.

It may not matter much in this case (the optimization flag is likely more important), but keep this in mind in the future, e.g., if you're passing around large matrices by value and they're constant, they can be passed around by const reference instead.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.