Logistic regression with LBFGS solver

Question

Is there any open source library or code which implements Logistic Regression using L-BFGS solver?

I would prefer Python, but other languages are welcome, too.

Answer 1

Here is an example of logistic regression estimation using the limited memory BFGS [L-BFGS] optimization algorithm. I will be using the optimx function from the optimx library in R, and SciPy's scipy.optimize.fmin_l_bfgs_b in Python.

Python

The example that I am using is from Sheather (2009, pg. 264). The following Python code shows estimation of the logistic regression using the BFGS algorithm:

# load required libraries
import numpy as np
import scipy as sp
import scipy.optimize
import pandas as pd
import os

# hyperlink to data location
urlSheatherData = "http://www.stat.tamu.edu/~sheather/book/docs/datasets/MichelinNY.csv"

# read in the data to a NumPy array
arrSheatherData = np.asarray(pd.read_csv(urlSheatherData))

# slice the data to get the dependent variable
vY = arrSheatherData[:, 0].astype('float64')

# slice the data to get the matrix of predictor variables
mX = np.asarray(arrSheatherData[:, 2:]).astype('float64')

# add an intercept to the predictor variables
intercept = np.ones(mX.shape[0]).reshape(mX.shape[0], 1)
mX = np.concatenate((intercept, mX), axis = 1)

# the number of variables and obserations
iK = mX.shape[1]
iN = mX.shape[0]

# logistic transformation
def logit(mX, vBeta):
    return((np.exp(np.dot(mX, vBeta))/(1.0 + np.exp(np.dot(mX, vBeta)))))

# stable parametrisation of the cost function
def logLikelihoodLogitStable(vBeta, mX, vY):
    return(-(np.sum(vY*(np.dot(mX, vBeta) -
    np.log((1.0 + np.exp(np.dot(mX, vBeta))))) +
                    (1-vY)*(-np.log((1.0 + np.exp(np.dot(mX, vBeta))))))))

# score function
def likelihoodScore(vBeta, mX, vY):
    return(np.dot(mX.T,
                  (logit(mX, vBeta) - vY)))

#====================================================================
# optimize to get the MLE using the BFGS optimizer (numerical derivatives)
#====================================================================
optimLogitBFGS = sp.optimize.minimize(logLikelihoodLogitStable,
                                  x0 = np.array([10, 0.5, 0.1, -0.3, 0.1]),
                                    args = (mX, vY), method = 'BFGS',
                                    options={'gtol': 1e-3, 'disp': True})

print(optimLogitBFGS) # print the results of the optimisation

And this can easily be adapted to the scipy.optimize.fmin_l_bfgs_b function:

#====================================================================
# optimize to get the MLE using the L-BFGS optimizer (analytical derivatives)
#====================================================================
optimLogitLBFGS = sp.optimize.fmin_l_bfgs_b(logLikelihoodLogitStable,
                                  x0 = np.array([10, 0.5, 0.1, -0.3, 0.1]),
                                    args = (mX, vY), fprime = likelihoodScore,
                                    pgtol =  1e-3, disp = True)

print(optimLogitLBFGS) # print the results of the optimisation

R

Using the L-BFGS-B optimizer in R is just as simple. First the version with the BFGS algorithm:

library(optimx)

# read in the data
urlSheatherData = "http://www.stat.tamu.edu/~sheather/book/docs/datasets/MichelinNY.csv"
dfSheatherData = as.data.frame(read.csv(urlSheatherData, header = T))

# create the design matrices
vY = as.matrix(dfSheatherData['InMichelin'])
mX = as.matrix(dfSheatherData[c('Service','Decor', 'Food', 'Price')])

# add an intercept to the predictor variables
mX = cbind(rep(1, nrow(mX)), mX)

# the number of variables and observations
iK = ncol(mX)
iN = nrow(mX)

# define the logistic transformation
logit = function(mX, vBeta) {
  return(exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) )
}

# stable parametrisation of the log-likelihood function
# Note: The negative of the log-likelihood is being returned, since we will be
#       /minimising/ the function.
logLikelihoodLogitStable = function(vBeta, mX, vY) {
  return(-sum(
    vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta)))
    + (1-vY)*(-log(1 + exp(mX %*% vBeta)))
  )  # sum
  )  # return
}

# score function
likelihoodScore = function(vBeta, mX, vY) {
  return(t(mX) %*% (logit(mX, vBeta) - vY) )
}

# initial set of parameters
vBeta0 = c(10, -0.1, -0.3, 0.001, 0.01)  # arbitrary starting parameters

#====================================================================
# optimize to get the MLE using the BFGS optimizer (numerical derivatives)
#====================================================================
optimLogitBFGS = optim(vBeta0, logLikelihoodLogitStable,
                    mX = mX, vY = vY, method = 'BFGS', hessian=TRUE)
optimLogitBFGS # get the results of the optimisation

and then the version with the L-BFGS-B from the optimx package:

#====================================================================
# optimize to get the MLE using the L-BFGS optimizer (analytical derivatives)
#====================================================================
optimLogitLBFGS = optimx(vBeta0, logLikelihoodLogitStable, method = 'L-BFGS-B',
                            gr = likelihoodScore, mX = mX, vY = vY, hessian=TRUE)

summary(optimLogitLBFGS)

Answer 2

If you're worrying about memory I guess you're either working with embedded hardware or expecting to have a big model. I'm going to guess that it's the latter and that you have a high dimensional text or bioinformatics classification problem of some sort. If so you should ponder Mallet's Java implementation, since that plugs into their relevant logistic regression (a.k.a. maxent) models most easily.

L-BFGS as a standalone algorithm is available in Java, Python, C and fortran implementations, handily linked from the L-BFGS wikipedia page. The Python (SciPy) version will presumably be of most interest to you. Applying this to a logistic regression models is relatively straightforward, except perhaps for the part where you choose a regulariser. Full disclosure: I do not use SciPy.

In logistic regression applications fancy regularisation and a limited-memory optimisation process, while conceptually separate, are often needed together due to the nature of the problem. Hence there's some reason to choose a library that bundles the two together in a sensible manner.

Answer 3

The Apache Spark compute engine is open source and has great performance on very large datasets. As of version 1.2 (I think) from 2014, Spark MLlib supports LogisticRegressionWithLBFGS. The API has bindings for Python, Scala or Java.

It uses feature scaling and L2-Regularization by default, unlike the gsm method in R.

There is an explanation with example code at Linear Methods - MLlib - Spark Documentation. The documentation license is CC BY-SA 3.0 US, so here is a snippet.

from pyspark.mllib.regression import LabeledPoint, LinearRegressionWithSGD
from numpy import array

# Load and parse the data
def parsePoint(line):
values = [float(x) for x in line.replace(',', ' ').split(' ')]
return LabeledPoint(values[0], values[1:])

data = sc.textFile("data/mllib/ridge-data/lpsa.data")
parsedData = data.map(parsePoint)

# Build the model
model = LinearRegressionWithSGD.train(parsedData)

# Evaluate the model on training data
valuesAndPreds = parsedData.map(lambda p: (p.label, model.predict(p.features)))
MSE = valuesAndPreds.map(lambda (v, p): (v - p)**2).reduce(lambda x, y: x + y) / valuesAndPreds.count()
print("Mean Squared Error = " + str(MSE))

Answer 4

Sk-learn has an excellent Logistic Regression implementation. It's just a wrapper around [LIBLINEAR], but LIBLINEAR is state-of-the-art and although it doesn't use LBFGS, it uses something else called dual-coordinate descent, which according to this paper is even better in many situations.

An alternate supposedly Python friendly implementation that does include LBFGS is Le Zhang's Maximum Entropy Toolkit although I haven't used it yet.

Answer 5

From here, http://www.kazanovaforanalytics.com/download.html, You can download a .jar with an implementation of logistic regression via newton raphson method that minimizes the -2 log likelihood. A comprehensive example can be found here :

http://www.kazanovaforanalytics.com/example_classes.txt

Provided via Apache licence 2.0, so you can include it in commercial applications.