Logistic Regression

Notebook version: 1.0 (Oct 16, 2015)

Author: Jesús Cid Sueiro (jcid@tsc.uc3m.es)

Changes: v.1.0 - First version

1. Introduction

1.1. Binary classification and decision theory

The goal of a classification problem is to assign a class or category to every instance or observation of a data collection. Here, we will assume that every instance ${\bf x}$ is an $N$-dimensional vector in $\mathbb{R}^N$, and that the class $y$ of sample ${\bf x}$ is an element of a binary set ${\mathcal Y} = \{0, 1\}$. The goal of a classifier is to predict the true value of $y$ after observing ${\bf x}$. We will denote as $\hat{y}$ the classifier output or decision. If $y=\hat{y}$, the decision is an hit, otherwise $y\neq \hat{y}$ and the decision is an error.

Decision theory provides a solution to the classification problem in situations where the relation between instance ${\bf x}$ and its class $y$ is given by a known probabilistic model: assume that every tuple $({\bf x}, y)$ is an outcome of a random vector $({\bf X}, Y)$ with joint distribution $p_{{\bf X},Y}({\bf x}, y)$. A natural criteria for classification is to select predictor $\hat{Y}=f({\bf x})$ in such a way that the probability or error, $P\{\hat{Y} \neq Y\}$ is minimum. Noting that $$ P\{\hat{Y} \neq Y\} = \int P\{\hat{Y} \neq Y | {\bf x}\} p_{\bf X}({\bf x}) d{\bf x} $$ the optimal decision is got if, for every sample ${\bf x}$, we make decision minimizing the condicional error probability: \begin{align} \hat{y}^* &= \arg\min_{\hat{y}} P\{\hat{y} \neq Y |{\bf x}\} \\ &= \arg\max_{\hat{y}} P\{\hat{y} = Y |{\bf x}\} \\ \end{align} Thus, the optimal decision rule can be expressed as $$ P_{Y|X}(1|{\bf x}) \quad\mathop{\gtrless}^{\hat{y}=1}_{\hat{y}=0}\quad P_{Y|X}(0|{\bf x}) $$ or, equivalently $$ P_{Y|X}(1|{\bf x}) \quad\mathop{\gtrless}^{\hat{y}=1}_{\hat{y}=0}\quad \frac{1}{2} $$ The classifier implementing this decision rule is usually named MAP (Maximum A Posterior).

1.2. Parametric classification.

Classical decision theory is grounded on the assumption that the probabilistic model relating the observed sample ${\bf X}$ and the true hypothesis $Y$ is known. Unfortunately, this is unrealistic in many applications, where the only available information to construct the classifier is a dataset $\mathcal S = \{({\bf x}^{(k)}, y^{(k)}), \,k=1,\ldots,K\}$ of intances and their respective class labels.

A more realistic formulation of the classification problem is the following: given a dataset $\mathcal S = \{({\bf x}^{(k)}, y^{(k)}) \in {\mathbb{R}}^N \times {\mathcal Y}, \, k=1,\ldots,K\}$ of independent and identically distributed (i.i.d.) samples from an unknown distribution $p_{{\bf X},Y}({\bf x}, y)$, predict the class $y$ of a new sample ${\bf x}$ with the minimum probability of error.

Since the probabilistic model generating the data is unknown, the MAP decision rule cannot be applied. However, many classification algorithms use the dataset to obtain an estimate of the posterior class probabilities, and apply it to implement an approximation to the MAP decision maker.

Parametric classifiers based on this idea assume, additionally, that the posterior class probabilty satisfies some parametric formula: $$ P_{Y|X}(1|{\bf x},{\bf w}) = f_{\bf w}({\bf x}) $$ where ${\bf w}$ is a vector of parameters. Therefore, the dataset ${\mathcal S}$ is used to compute an estimate $\hat{\bf w}$ of ${\bf w}$ and then apply the decision rule $$ f_{\hat{\bf w}}({\bf x}) \quad\mathop{\gtrless}^{\hat{y}=1}_{\hat{y}=0}\quad \frac{1}{2} $$

In this lesson we explore one of the most popular model-based parametric classification methods: logistic regression.

2. Logistic regression.

2.1. The logistic function

The logistic regression model assumes that the binary class label $Y \in \{0,1\}$ of observation $X\in \mathbb{R}^N$ satisfies the expression. $$p_{Y|{\bf X},W}(1|{\bf x}, {\bf w}) = f({\bf w}^\intercal{\bf x})$$ $$p_{Y|{\bf,X},W}(0|{\bf x}, {\bf w}) = 1-f({\bf w}^\intercal{\bf x})$$ where ${\bf w}$ is a parameter vector and $f(·)$ is the logistic function, which is defined by $$f(t) = \frac{1}{1+\exp(-t)}$$

It is straightforward to see that the logistic function has the following properties:

• P1: Probabilistic output: $\quad 0 \le f(t) \le 1$
• P2: Symmetry: $\quad f(-t) = 1-f(t)$
• P3: Monotonicity: $\quad f'(t) = f(t)·(1-f(t)) \ge 0$

In the following we define a logistic function in python, and use it to plot a graphical representation.

Exercise 1: Verify properties P2 and P3.

2.2. Classifiers based on the logistic model.

The MAP classifier under a logistic model will have the form $$p_{Y|{\bf X},W}(1|{\bf x}, {\bf w}) \quad\mathop{\gtrless}^{\hat{y}=1}_{\hat{y}=0} \quad \frac{1}{2} $$ that is $$f({\bf w}^\intercal{\bf x}) \quad\mathop{\gtrless}^{\hat{y}=1}_{\hat{y}=0} \quad \frac{1}{2} $$ therefore $$ 2 \quad\mathop{\gtrless}^{\hat{y}=1}_{\hat{y}=0} \quad 1 + \exp(-{\bf w}^\intercal{\bf x}) $$ which is equivalent to $${\bf w}^\intercal{\bf x} \quad\mathop{\gtrless}^{\hat{y}=1}_{\hat{y}=0}\quad 0 $$ Therefore, the classifiers based on the logistic regression model are given by linear decision boundaries passing through the origin, ${\bf x} = {\bf 0}$.

3.3. Nonlinear classifiers.

The logistic model can be extended to construct non-linear classifiers by using non-linear data transformations. A general form for a nonlinear logistic regression model is

$$p_{Y|{\bf X},W}(1|{\bf x}, {\bf w}) = f({\bf w}^\intercal{\bf z}({\bf x})) $$

where ${\bf z}({\bf x})$ is an arbitrary nonlinear transformation. The boundary decision in that case is given by equation

$$ {\bf w}^\intercal{\bf z}({\bf x}) = 0 $$

Exercise 2: Modify the code above to generate a 3D surface plot of the polynomial logistic regression model given by

$$ P(1|{\bf x}) = f(1 + 10 x_0 + 10 x_1 - 20 x_0^2 + 5 x_0 x_1 + x_1^2) $$

3. Inference

Remember that the starting idea of parametric classification is to use data in a training set $\mathcal S = \{({\bf x}_k, y_k) \in {\mathbb{R}}^N \times \{0,1\}, k=1,\ldots,K\}$ to estimate the unknown parameter ${\bf w}$, and use the estimate $\hat{\bf w}$ to compute the label prediction for any new observation as

$$\hat{y} = \arg\max_y P_{Y|{\bf X},{\bf W}}(y|{\bf x},\hat{\bf w}).$$

To do so, we will make the following assumptions:

• A1. The samples in ${\mathcal S}$ are i.i.d.
• A2. Target $Y^{(k)}$ only depends on ${\bf x}^{(k)}$, but not on ${\bf x}^{(l)}$ for any $l\neq k$.
• A3. (Logistic Regression): target $Y^{(k)}$ satisfies a logistic model: $$P_{Y|{\bf X},{\bf W}}(1|{\bf x}, {\bf w}) = f({\bf w}^\intercal{\bf z}({\bf x})).$$

In the following we discuss different approaches to the estimation of ${\bf w}$.

• Maximum Likelihood (ML) estimat: ${\bf w}_{\text{ML}} = \arg\max_{\bf w} p_{{\mathcal S}|{\bf W}}({\mathcal S}|{\bf w})$
• Maximum A Posteriori estimate (MAP): ${\bf w}_{\text{MAP}} = \arg\max_{\bf w} p_{{\bf W}|{\mathcal S}}({\bf w}|{\mathcal S})$
• Bayesian estimate: ${\bf w}_{\text{B}} = \mathbb{E}\{{\bf W}|{\mathcal S}\}$

For the mathematical derivation of the logistic regression algorithm, the following representation of the logistic model will be useful: noting that

$$p_{Y|{\bf,X},W}(0|{\bf x}, {\bf w}) = 1-f({\bf w}^\intercal{\bf z}({\bf x})) = f(-{\bf w}^\intercal{\bf z}({\bf x}))$$

we can write

$$p_{Y|{\bf X},W}(y|{\bf x}, {\bf w}) = f(\overline{y}{\bf w}^\intercal{\bf z}({\bf x}))$$

where $\overline{y} = 2y-1$ is a symmetrized label ($\overline{y}\in\{-1, 1\}$).

3.1. ML estimation.

The ML estimate is defined as

$${\bf w}_{\text{ML}} = \arg\min_{\bf w} L({\bf w}) $$

where $L({\bf w})$ is the negative log-likelyhood function, given by

$$ L({\bf w}) = - \log\left(p\left(y^{(1)},\ldots,y^{(K)}| {\bf x}^{(1)},\ldots, {\bf x}^{(K)},{\bf w}\right)\right) $$

This function can be computed as follows. Using assumption A1,

$$ L({\bf w}) = - \log\left(\prod_{k=1}^K p\left(y^{(k)}|{\bf x}^{(1)},\ldots,{\bf x}^{(K)},{\bf w}\right)\right), $$

using A2,

\begin{align} L({\bf w}) &= - \log\left(\prod_{k=1}^K p_{Y|{\bf X},{\bf W}}\left(y^{(k)}|{\bf x}^{(k)},{\bf w}\right)\right) \\ &= - \sum_{k=1}^K\log\left(p_{Y|{\bf X},{\bf W}}\left(y^{(k)}|{\bf x}^{(k)},{\bf w}\right)\right) \end{align}

and, using A3 (the logistic model)

\begin{align} L({\bf w}) &= - \sum_{k=1}^K\log\left(f\left(\overline{y}^{(k)}{\bf w}^\intercal {\bf z}^{(k)}\right)\right) \\ &= \sum_{k=1}^K\log\left(1+\exp\left(-\overline{y}^{(k)}{\bf w}^\intercal {\bf z}^{(k)}\right)\right) \end{align}

where ${\bf z}^{(k)}={\bf z}({\bf x}^{(k)})$.

It can be shown that $L({\bf w})$ is a convex and differentiable function of ${\bf w}$. Therefore, its minimum is a point with zero gradient.

\begin{align} \nabla_{\bf w} L(\hat{\bf w}_{\text{ML}}) &= - \sum_{k=1}^K \frac{\exp\left(-\overline{y}^{(k)}\hat{\bf w}_{\text{ML}}^\intercal {\bf z}^{(k)}\right) \overline{y}^{(k)} {\bf z}^{(k)}} {1+\exp\left(-\overline{y}^{(k)}\hat{\bf w}_{\text{ML}}^\intercal {\bf z}^{(k)} \right)} = \\ &= - \sum_{k=1}^K \left(y^{(k)}-f(\hat{\bf w}_{\text{ML}}^T {\bf z}^{(k)})\right) {\bf z}^{(k)} = 0 \end{align}

Unfortunately, $\hat{\bf w}_{\text{ML}}$ cannot be taken out from the above equation, and some iterative optimization algorithm must be used to estimate the minimum.

3.2. Gradient descent.

A simple iterative optimization algorithm is gradient descent.

\begin{align} {\bf w}_{n+1} = {\bf w}_n - \rho_n \nabla_{\bf w} L({\bf w}_n) \end{align}

where $\rho_n >0$ is the learning step.

Applying the gradient descent rule to logistic regression, we get the following algorithm:

\begin{align} {\bf w}_{n+1} &= {\bf w}_n + \rho_n \sum_{k=1}^K \left(y^{(k)}-f({\bf w}_n^\intercal {\bf z}^{(k)})\right) {\bf z}^{(k)} \end{align}

Defining vectors

\begin{align} {\bf y} &= (y^{(1)},\ldots,y^{(K)})^\intercal \\ \hat{\bf p}_n &= (f({\bf w}_n^\intercal {\bf x}^{(1)}), \ldots, f({\bf w}_n^\intercal {\bf z}^{(K)}))^\intercal \end{align}

and matrix \begin{align} {\bf Z} = \left({\bf z}^{(1)},\ldots,{\bf z}^{(K)}\right)^\intercal \end{align}

we can write

\begin{align} {\bf w}_{n+1} &= {\bf w}_n + \rho_n {\bf Z}^\intercal \left({\bf y}-\hat{\bf p}_n\right) \end{align}

In the following, we will explore the behavior of the gradient descend method using the Iris Dataset.

3.2.1 Example: Iris Dataset.

As an illustration, consider the Iris dataset , taken from the UCI Machine Learning repository. This data set contains 3 classes of 50 instances each, where each class refers to a type of iris plant (setosa, versicolor or virginica). Each instance contains 4 measurements of given flowers: sepal length, sepal width, petal length and petal width, all in centimeters.

We will try to fit the logistic regression model to discriminate between two classes using only two attributes.

First, we load the dataset and split them in training and test subsets.

Now, we select two classes and two attributes.

3.2.2. Data normalization

Normalization of data is a common pre-processing step in many machine learning algorithms. Its goal is to get a dataset where all input coordinates have a similar scale. Learning algorithms usually show less instabilities and convergence problems when data are normalized.

We will define a normalization function that returns a training data matrix with zero sample mean and unit sample variance.

Now, we can normalize training and test data. Observe in the code that the same transformation should be applied to training and test data. This is the reason why normalization with the test data is done using the means and the variances computed with the training set.

The following figure generates a plot of the normalized training data.

In order to apply the gradient descent rule, we need to define two methods:

• A fit method, that receives the training data and returns the model weights and the value of the negative log-likelihood during all iterations.
• A predict method, that receives the model weight and a set of inputs, and returns the posterior class probabilities for that input, as well as their corresponding class predictions.

We can test the behavior of the gradient descent method by fitting a logistic regression model with ${\bf z}({\bf x}) = (1, {\bf x}^\intercal)^\intercal$.

3.2.3. Free parameters

Under certaing conditions, the gradient descent method can be shown to converge assymtotically (i.e. as the number of iterations goes to infinity) to the ML estimate of the logistic model. However, in practice, the final estimate of the weights ${\bf w}$ depend on several factors:

• Number of iterations
• Initialization
• Learning step

Exercise: Visualize the variability of gradient descent caused by initializations. To do so, fix the number of iterations to 200 and the learning step, and execute the gradient descent 100 times, storing the training error rate of each execution. Plot the histogram of the error rate values.

Note that you can do this exercise with a loop over the 100 executions, including the code in the previous code slide inside the loop, with some proper modifications. To plot a histogram of the values in array p with nbins, you can use plt.hist(p, n)

3.2.3.1. Learning step

The learning step, $\rho$, is a free parameter of the algorithm. Its choice is critical for the convergence of the algorithm. too large values of $\rho$ make the algorithm diverge. For too small values, the convergence is too slown and more iterations are required for a good convergence.

Exercise 3: Observe the evolution of the negative log-likelihood with the number of iterations, for different values of $\rho$. It is easy to check that, for $\rho$ large enough, the gradient descent method does not converge. Can you estimate (through manual observation) and approximate value of $\rho$ stating a boundary between convergence and non-convergence?

Exercise 4: In this exercise we explore the influence of the learning step more sistematically. Use the code in the previouse exercises to compute, for every value of $\rho$, the average error rate over 100 executions. Plot the average error rate vs. $\rho$.

Note that you should explore the values of $\rho$ in a logarithmic scale. For instance, you can take $\rho = 1, 1/10, 1/100, 1/1000, \ldots$

In practice, the selection of $\rho$ may be a matter of trial an error. Also there is some theoretical evidence that the learning step should decrease along time up to cero, and the sequence $\rho_n$ should satisfy two conditions:

• C1: $\sum_{n=0}^{\infty} \rho_n^2 < \infty$ (decrease slowly)
• C2: $\sum_{n=0}^{\infty} \rho_n = \infty$ (do not decrease too much slowly)

For instance, we can take $\rho_n= 1/n$. Another common choice is $\rho_n = \alpha/(1+\beta n)$ where $\alpha$ and $\beta$ are also free parameters that can be selected by trial and error with some heuristic method.

3.2.4. Visualizing the posterior map.

We can also visualize the posterior probability map estimated by the logistic regression model for the estimated weights.

3.2.5. Polynomial Logistic Regression

The error rates of the logitic regression model can be potentially reduce by using polynomial transformations.

To compute the polinomial transformation up to a given degree, we can use the PolynomialFeatures method in sklearn.preprocessing.

Visualizing the posterior map we can se that the polynomial transformations produces nonlinear decision boundaries.

4. Regularization and MAP estimation.

An alternative to the ML estimation of the weights in logistic regression is Maximum A Posteriori estimation. Modelling the logistic regression weights as a random variable with prior distribution $p_{\bf W}({\bf w})$, the MAP estimate is defined as

$$ \hat{\bf w}_{\text{MAP}} = \arg\max p({\bf w}|{\mathcal S}) $$

The posterior density $p({\bf w}|{\mathcal S})$ is related to the likelihood function and the prior density of the weights, $p_{\bf W}({\bf w})$ through the Bayes rule

$$ p({\bf w}|{\mathcal S}) = \frac{p\left(y^{(1)},\ldots,y^{(K)}|{\bf x}^{(1)},\ldots, {\bf x}^{(K)},{\bf w}\right) p_{\bf W}({\bf w})} {p\left(y^{(1)},\ldots,y^{(K)}|{\bf x}^{(1)},\ldots, {\bf x}^{(K)}\right)} $$

In general, the denominator in this expression cannot be computed analytically. However, it is not required for MAP estimation because it does not depend on ${\bf w}$. Thus, we can compute the MAP estimate as \begin{align} \hat{\bf w}{\text{MAP}} &= \arg\max{\bf w} p\left(y^{(1)},\ldots,y^{(K)}|{\bf x}^{(1)},\ldots, {\bf x}^{(K)},{\bf w}\right) p{\bf W}({\bf w}) \ &= \arg\max{\bf w} \left{ \log\left(p\left(y^{(1)},\ldots,y^{(K)}|{\bf x}^{(1)},\ldots, {\bf x}^{(K)},{\bf w}\right) \right)

     + \log\left( p_{\bf W}({\bf w})\right)
     \right\}  \\

&= \arg\min{\bf w} \left{L({\bf w}) - \log\left( p{\bf W}({\bf w})\right) \right} \end{align} where $L(·)$ is the negative log-likelihood function. Two common choices for the prior are

• Gaussian: ${\bf W}$ is a zero-mean Gaussian random variable with variance matrix $v{\bf I}$. Thus: $$ p_{\bf W}({\bf w}) = \frac{1}{(2\pi v)^{N/2}} \exp\left(-\frac{1}{2v}\|{\bf w}\|^2\right) $$ and the MAP estimate is \begin{align} \hat{\bf w}{\text{MAP}} &= \arg\min{\bf w} \left{L({\bf w}) + \frac{1}{C}|{\bf w}|^2

     \right\}
  

  \end{align} where $C = 2v$. Noting that $$\nabla_{\bf w}\left\{L({\bf w}) + \frac{1}{C}\|{\bf w}\|^2\right\} = - {\bf Z}^\intercal \left({\bf y}-\hat{\bf p}_n\right) + \frac{2}{C}{\bf w}, $$ the gradient descent rule for MAP estimation becomes \begin{align} {\bf w}_{n+1} &= \left(1-\frac{2\rho_n}{C}\right){\bf w}_n

  • \rho_n {\bf Z}^\intercal \left({\bf y}-\hat{\bf p}_n\right) \end{align}

• Laplacian: ${\bf W}$ is a multivariate zero-mean Laplacian distribution given by $$ p_{\bf W}({\bf w}) = \frac{1}{(2 C)^{N}} \exp\left(-\frac{1}{C}\|{\bf w}\|_1\right) $$ (where $\|{\bf w}\|=|w_1|+\ldots+|w_N|$ is the $L_1$ norm of ${\bf w}$) and the MAP estimate is

\begin{align} \hat{\bf w}_{\text{MAP}} &= \arg\min_{\bf w} \left\{L({\bf w}) + \frac{1}{C}\|{\bf w}\|_1 \right\} \end{align}

The additional term introduced by the prior in the optimization algorithm is usually named the regularization term. It is usually very effective to avoid overfitting when the dimension of the weight vectors is high. Parameter $C$ is named the inverse regularization strength.

Exercise 5: Derive the gradient descent rules for MAP estimation of the logistic regression weights with Laplacian prior.

5. Other optimization algorithms

5.1. Stochastic Gradient descent.

Stochastic gradient descent (SGD) is based on the idea of using a single sample at each iteration of the learning algorithm. The SGD rule for ML logistic regression is

\begin{align} {\bf w}_{n+1} &= {\bf w}_n + \rho_n {\bf z}^{(n)} \left(y^{(n)}-\hat{p}^{(n)}_n\right) \end{align}

Once all samples in the training set have been applied, the algorith can continue by applying the training set several times.

The computational cost of each iteration of SGD is much smaller than that of gradient descent, though it usually needs more iterations to converge.

Exercise 5: Modify logregFit to implement an algorithm that applies the SGD rule one to every sample in the training set.

5.2. The Newton's method

Assume that the function to be minimized, $C({\bf w})$, can be approximated by its second order Taylor series expansion around ${\bf w}_0$

$$ C({\bf w}) \approx C({\bf w}_0) + \nabla_{\bf w}^\intercal C({\bf w}_0)({\bf w}-{\bf w}_0) + \frac{1}{2}({\bf w}-{\bf w}_0)^\intercal{\bf H}({\bf w}_0)({\bf w}-{\bf w}_0) $$

where ${\bf H}({\bf w}_k)$ is the *Hessian* matrix of $C$ at ${\bf w}_k$. Then, the minimum of C around ${\bf w}_0$ can be approximated as $$ {\bf w}^* = {\bf w}_0 - {\bf H}({\bf w}_0)^{-1} \nabla_{\bf w}^\intercal C({\bf w}_0) $$

Since the second order polynomial is only an approximation to $C$, ${\bf w}^*$ is only an approximation to the optimal weight vector, but we can expect ${\bf w}^*$ to be closer to the minimizer of $C$ than ${\bf w}_0$. Thus, we can repeat the process, computing a second order approximation around ${\bf w}^*$ and a new approximation to the minimizer.

The Newton's method is based on this idea. At each optization step, the function to be minimized is approximated by a second order approximation using a Taylor series expansion around the current estimate. As a result, the learning rules becomes

$$\hat{\bf w}_{n+1} = \hat{\bf w}_{n} - \rho_n {\bf H}({\bf w}_k)^{-1} \nabla_{{\bf w}}C({\bf w}_k) $$

For instance, for the MAP estimate with Gaussian prior, the Hessian matrix becomes

$$ {\bf H}({\bf w}) = \frac{2}{C}{\bf I} + \sum_{k=1}^K f({\bf w}^T {\bf z}^{(k)}) \left(1-f({\bf w}^T {\bf z}^{(k)})\right){\bf z}^{(k)} ({\bf z}^{(k)})^\intercal $$

Defining diagonal matrix

$$ {\mathbf S}({\bf w}) = \text{diag}\left(f({\bf w}^T {\bf z}^{(k)}) \left(1-f({\bf w}^T {\bf z}^{(k)})\right)\right) $$

the Hessian matrix can be written in more compact form as

$$ {\bf H}({\bf w}) = \frac{2}{C}{\bf I} + {\bf Z}^\intercal {\bf S}({\bf w}) {\bf Z} $$

Therefore, the Newton's algorithm for logistic regression becomes

\begin{align} \hat{\bf w}_{n+1} = \hat{\bf w}_{n} + \rho_n \left(\frac{2}{C}{\bf I} + {\bf Z}^\intercal {\bf S}(\hat{\bf w}_{n}) {\bf Z} \right)^{-1} {\bf Z}^\intercal \left({\bf y}-\hat{\bf p}_n\right) \end{align}

Some variants of the Newton method are implemented in the Scikit-learn package.

6. Logistic regression in Scikit Learn.

The scikit-learn package includes an efficient implementation of logistic regression. To use it, we must first create a classifier object, specifying the parameters of the logistic regression algorithm.