Logistic regression is a workhorse in data mining. Like decision trees and SVMs, it is a very standard classifier. If you have a labeled data, logistic regression definitely is one of the classifiers that should tried.
In this post, I’d like to show how to implement a logistic regression using Microsoft Solver Foundation in F#. Because we use the optimization in Microsoft Solver Foundation, the core part of a logistic regression contains only 20 – 30 lines of code.
If you are only interested in F# and L-BFGS, but not logistic regression, please only read the second section: L-BFGS solver in Microsoft Solver Foundation.
Background for logistic regression
Logistic regression model predicts the probability that a data instance x being labeled as y
![clip_image002[5]](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgo2PDkBzgDI1uMO4oWcP4hwBw-xnbTSrB9x9bKXqQSFeHNicoW3QE7gF8BIsp9UjHEOfRWBLL6dO-15grKCM1uU9ePCjVHDC1EGUe89DHHg2rNLcD7tVDYkRc_lhUzoF_FBdsdSeSBTP-P/s1600-h/clip_image002%5B5%5D%5B4%5D.png "clip_image002[5]"){kind=small}
A data instance is a vector with each dimension representing a value for some feature. For example, the Iris flower data set we are going to use in this post has 4 features: Sepal Length, Sepal Width, Petal Length and Petal Width. Each instance has a label y indicating the category. The above model only models binary labels. Logistic regression also has multiple categories version, see [Mitchell]. In this post, we focus on the binary version; however multiple versions should be straightforward once we know how to implement the binary one.
The model contains two parameters weight vector w and offset b, which should be estimated/learned on the training data set![clip_image004[3]](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhvMRU3VmhzYgEXU0a5HXRt_-GyfaU8D0CFa2A80uoVO3-hHJo1x3KEQfIV52XinCQTgUROKzHm9uGpCpu21_j9CaAP5qA0iDGIQ0qnacsc3j9KZo2XgyJRkxdX-dX4pAy6ozVGlsxcVcF0/s1600-h/clip_image004%5B3%5D%5B3%5D.png "clip_image004[3]"){kind=small}
If we assume that all the instances are independent and identically distributed (i.i.d.). Then the regularized joint probability is
![clip_image006[3]](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjzto2mO3bnBcpLtO5q8JQrd26hH86IY4_4eOJSzWlrtidbRZlt5bLJTbNkMaKxJ4B5lzBDdpov6zAKwGEYWiGfbcJosiw4ENaOp2wl4U7QgH33QcweiH99f9zIDn4M4Qgmvv70espLhAOY/s1600-h/clip_image006%5B3%5D%5B3%5D.png "clip_image006[3]")
where ![clip_image008[1]](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhBR4s76Tj_DrQar4J1F_P0aIgZWXz0gvjDZCWQagRcX-UVDcopqWKycBMGhNYZ-eN0Oj8EzcBQzDJ1cLOjHjXKgZexHygPgxicyWQMnB-U7tPRJyuV-qhrtF4Q4MBrLzbW9qxNtj-Ice3M/s1600-h/clip_image008%5B1%5D%5B3%5D.png "clip_image008[1]"){kind=small}
The log version of this likelihood is
![clip_image010[1]](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjGs6SEKhNaoLUbc3ARMeijjnTzS9_qnPdRfM9lqq-2uyTHd1zU19gwebaFFgfF2wpWjowL6cmEzMsUJozOM0GT1SaZ7gOGceLudA7r2VJACAewxgTeMyaVUZC-HlrqVeJzzhtFOHaNTT2R/s1600-h/clip_image010%5B1%5D%5B4%5D.png "clip_image010[1]")
We want to find the model parameter w that maximizes this likelihood. P(X) is a concave function with respect to w, so it has a global maximal.
What left is how to optimize this function. BFGS is a Quasi-Newton method that only needs first order gradient while a Newton family method needs second order explicitly. Microsoft Solver Foundation has a limited memory version of BFGS (L-BFGS). (L-BFGS is more popular than BFGS, sometimes people mention BFGS to mean L-BFGS.) There is a memory parameter m in L-BFGS, the space complexity of this algorithm is O(md), where d is the dimension of the variable space. Say, if you want to optimize a function that has 1,000,000 free variables, and set m to 10, then you need only a little more than 1M*10*8/1024^2=76M bytes.
The input for an L-BFGS optimizer is: 1) the function evaluator L(w) and 2) the gradient evaluator 
The output is the best 
L-BFGS solver in Microsoft Solver Foundation
Microsoft Solver Foundation already implements an L-BFGS solver. The reader is suggested to read pages 24-27 & 72-75 of MSF-SolverProgrammingPrimer.pdf (on my machine: C:\Program Files\Microsoft Solver Foundation\3.0.1.10599\Documents) for detailed reference. Because the code for optimizing Rosenbrock function is in C#, here I translated it to F#. Note how I translate between delegate type System.Func<> and F# functions.
Btw, I find that using the Microsoft.Solver.Foundatoin.dll in Sho would be more convenient as it packages everything into one DLL.
#r @"C:\Program Files (x86)\Sho 2.0 for .NET 4\packages\Optimizer\Microsoft.Solver.Foundation.dll"
openMicrosoft.SolverFoundation.Common
openMicrosoft.SolverFoundation.Solvers
openMicrosoft.SolverFoundation.Services
letsolverParams = newCompactQuasiNewtonSolverParams()
letsolver = newCompactQuasiNewtonSolver()
//add variables
let_, vidVaribaleX = solver.AddVariable(null)
let_, vidVaribaleY = solver.AddVariable(null)
//add a row and set it as the goal
let_, vidRow = solver.AddRow(null)
solver.AddGoal(vidRow, 0, true)
letOriginalRosenbrockFunction =
letf (model:INonlinearModel) (rowVid:int) (values:ValuesByIndex) (newValues:bool) =
Math.Pow(1. - values.[1], 2.) + 100. * (Math.Pow(values.[2] - (values.[1] * values.[1]), 2.))
newSystem.Func<INonlinearModel, int, ValuesByIndex, bool, float> (f)
letOriginalRosenbrockGradient =
letf (model:INonlinearModel) (rowVid:int) (values:ValuesByIndex) (newValues:bool) (gradient:ValuesByIndex) =
gradient.[1] <- -2. * (1. - values.[1]) - 400. * values.[1] * (values.[2] - (values.[1] * values.[1]))
gradient.[2] <- 200. * (values.[2] - (values.[1] * values.[1]))
newSystem.Action<INonlinearModel, int, ValuesByIndex, bool, ValuesByIndex> (f)
solver.FunctionEvaluator <- OriginalRosenbrockFunction
solver.GradientEvaluator <- OriginalRosenbrockGradient
solver.Solve(solverParams);
Console.WriteLine(solver.ToString());
and the output:
>
Minimize problem
Dimensions = 2
Variable indexes:
1, 2
Starting point = 0, 0
Solution quality is: LocalOptima
Number of iterations performed: 21
Number of evaluation calls: 27
Finishing point =0.999999997671584, 0.999999995285991
Finishing value =5.74844710320033E-18
val it : unit = ()
>
Implementation for logistic regression
We first get a sample data set online:
module Net =
open System.Net
let fetchUrlSimple (url:string) =
let req = WebRequest.Create(url)
let response = req.GetResponse()
use stream = response.GetResponseStream()
use streamreader = new System.IO.StreamReader(stream)
streamreader.ReadToEnd()
type dataset =
{ features: float array array; // (instance = float array) array
mutable labels: int array; //
}
let iris =
let page = Net.fetchUrlSimple @"http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data"
let lines = page.Split([|'\n'|], StringSplitOptions.RemoveEmptyEntries)
let getIrisId (s:string) =
if s = "Iris-setosa" then
2
elif s = "Iris-versicolor" then
-1
else
1
let instances =
lines
|> Array.map (fun line ->
let s = line.Split ','
let features = [| "1.0"; s.[0]; s.[1]; s.[2]; s.[3] |] |> Array.map float // add a dumy 1.0 into the fatures
let label = getIrisId s.[4]
features, label
)
|> Array.filter (fun (_, label) -> label <> 2)
let F, L = Array.unzip instances
{
features = F;
labels = L;
}
The logistic regression implementation:
module LogReg =
let mutable lambda = 0.1
let dotProduct (x:float array) (g:ValuesByIndex) =
let mutable dot = 0.0
for i=0 to x.Length-1 do
dot <- dot + x.[i] * g.[i+1]
dot
let sigmoid (x:float array) (y:int) (g:ValuesByIndex) =
let mutable dot = 0.0
for i=0 to x.Length-1 do
dot <- dot + x.[i] * g.[i+1]
let z = (float y) * dot
(* code to deal with extream numbers, but with some issues with l(w) calculation
if z > 30.0 then
1.0
elif z < -30.0 then
0.0
else
*)
1.0 / (1.0 + exp (- z))
let sigmoid2 (x:float array) (y:int) (g:float array) =
let mutable dot = 0.0
for i=0 to x.Length-1 do
dot <- dot + x.[i] * g.[i]
1.0 / (1.0 + exp (- (float y) * dot))
let logregValue(ds:dataset) =
let dim = ds.features.[0].Length
let f (model:INonlinearModel) (rowVid:int) (values:ValuesByIndex) (newValues:bool) =
let mutable L = 0.0
for i=1 to dim do
L <- L + values.[i]*values.[i]
L <- - (L * lambda / 2.0)
for i=0 to ds.features.Length-1 do
//L <- L - log (1.0 + exp (- (float ds.labels.[i]) * dotProduct ds.features.[i] values ))
L <- L + log (sigmoid ds.features.[i] ds.labels.[i] values)
// printfn "L = %.10f" L
L
new System.Func<INonlinearModel, int, ValuesByIndex, bool, float> (f)
let logregGradient(ds:dataset) =
// printfn "gradient calculated"
let dim = ds.features.[0].Length
let f (model:INonlinearModel) (rowVid:int) (values:ValuesByIndex) (newValues:bool) (gradient:ValuesByIndex) =
for j=1 to dim do
gradient.[j] <- -lambda * values.[j]
for i=0 to ds.features.Length-1 do
let coef = (1.0 - sigmoid ds.features.[i] ds.labels.[i] values) * (float ds.labels.[i])
for j=1 to dim do
gradient.[j] <- gradient.[j] + coef * ds.features.[i].[j-1]
new System.Action<INonlinearModel, int, ValuesByIndex, bool, ValuesByIndex> (f)
let makeSolver(ds:dataset) =
// set the solver parameters
let solverParams = new CompactQuasiNewtonSolverParams()
//solverParams.IterationLimit <- 11
//solverParams.Tolerance <- 1e-4
// the memory parameter m, default 17
solverParams.IterationsToRemember <- 10
// solver
let solver = new CompactQuasiNewtonSolver()
let dim = ds.features.[0].Length
// add variables
for i=1 to dim do
solver.AddVariable(null) |> ignore
// add a row and set it as the goal
let _, vidRow = solver.AddRow(null)
solver.AddGoal(vidRow, 0, false) |> ignore
// set funcation evaluator and gradient evaluator
solver.FunctionEvaluator <- logregValue(ds)
solver.GradientEvaluator <- logregGradient(ds)
// solve!
solver.Solve(solverParams) |> ignore
// get the optimized point (w) in the solver
let w = Array.create dim 0.0
for i=1 to dim do
w.[i-1] <- solver.GetValue(i).ToDouble()
w
Try our logistic regression model on the Iris dataset:
let sol = LogReg.makeSolver(iris)
let nWrongs = ref 0
for i=0 to iris.features.Length-1 do
let prob = LogReg.sigmoid2 iris.features.[i] iris.labels.[i] sol
printfn "L = %A sig value = %A label = %A" iris.labels.[i] prob (if prob > 0.5 then "correct" else incr nWrongs; "wrong")
Reference
Tom Mitchell, Generative and Discriminative Classifiers: Naive Bayes and Logistic Regression, new book chapter. This chapter gives a good introduction to logistic regression. For a more theoretic comparison, see: Ng & Jordan: on discriminative vs generative classifiers a comparison of logistic regression and naïve bayes.
Thomas P. Minka, A comparison of numerical optimizers for logistic regression.
It is well written, the first page ready gives a concise introduction to logistic regression. The rest of the paper has a lot of discussion on optimization methods. Tom also has a Matlab package.
Galen Andrew and Jianfeng Gao. Scalable training of L1-regularized log-linear models. In Proceedings of the 24th International Conference on Machine Learning (ICML 2007), pp. 33-40, 2007.
L1-regularized logistic regression with an implementation in C++.
libLBFGS(LBFGS in C): http://www.chokkan.org/software/liblbfgs/
Note
I quite like this post because it combines so many good things together: F#, logistic regression, optimization and playing with sample datasets.
