这是使用有限内存 BFGS [L-BFGS] 优化算法进行逻辑回归估计的示例。我将在 R 中使用库中的optimx函数，在 Python 中使用optimxSciPy 中的函数。scipy.optimize.fmin_l_bfgs_b

Python

我使用的示例来自Sheather (2009, pg. 264)。以下 Python 代码显示了使用 BFGS 算法对逻辑回归的估计：

# 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

这可以很容易地适应scipy.optimize.fmin_l_bfgs_b功能：

#====================================================================
# 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

在 R 中使用 L-BFGS-B 优化器同样简单。首先是 BFGS 算法的版本：

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

然后是optimx包中带有 L-BFGS-B 的版本：

#====================================================================
# 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)

如果您担心内存，我猜您要么使用嵌入式硬件，要么期望拥有一个大模型。我猜是后者，你有某种高维文本或生物信息学分类问题。如果是这样，您应该考虑Mallet 的Java 实现，因为它最容易插入其相关的逻辑回归（又名 maxent）模型。

L-BFGS 作为独立算法可用于 Java、Python、C 和 fortran 实现，可从L-BFGS 维基百科页面轻松链接。Python (SciPy) 版本可能是您最感兴趣的。将其应用于逻辑回归模型相对简单，除了您选择正则化器的部分。完全披露：我不使用 SciPy。

在逻辑回归应用中，花哨的正则化和有限内存优化过程虽然在概念上是分开的，但由于问题的性质，通常需要一起使用。因此，有理由选择一个以合理的方式将两者捆绑在一起的库。

Apache Spark计算引擎是开源的，在非常大的数据集上具有出色的性能。从 2014 年的 1.2 版（我认为）开始，Spark MLlib 支持 LogisticRegressionWithLBFGS。该 API 具有 Python、Scala 或 Java 的绑定。

与 R 中的 gsm 方法不同，它默认使用特征缩放和 L2-Regularization 。

线性方法 - MLlib - Spark 文档中有示例代码的解释。文档许可证是 CC BY-SA 3.0 US，所以这里是一个片段。

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

Sk-learn具有出色的逻辑回归实现。它只是 [LIBLINEAR] 的一个包装器，但 LIBLINEAR 是最先进的，虽然它不使用 LBFGS，但它使用了另一种称为双坐标下降的东西，根据本文，这在许多情况下甚至更好。

另一个包含 LBFGS 的 Python 友好实现是 Le Zhang 的最大熵工具包，尽管我还没有使用它。