Logistic Regression

Notebook version: 1.0 (Oct 12, 2016)

Author: Jesús Cid Sueiro (jcid@tsc.uc3m.es)
        Jerónimo Arenas García (jarenas@tsc.uc3m.es)

Changes: v.1.0 - First version
         v.1.1 - Typo correction. Prepared for slide presentation

Logistic Regression

1. Introduction

1.1. Binary classification and decision theory. The MAP criterion

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 conditional 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|{\bf X}}(1|{\bf x}) \quad\mathop{\gtrless}^{\hat{y}=1}_{\hat{y}=0}\quad P_{Y|{\bf X}}(0|{\bf x}) $$

or, equivalently

$$ P_{Y|{\bf 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 Posteriori).

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 instances 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. Given the expression of the MAP decision maker, classification consists in comparing the value of $f_{\bf w}({\bf x})$ with the threshold $\frac{1}{2}$, and each parameter vector would be associated to a different decision maker.

In practice, the dataset ${\mathcal S}$ is used to select a particular parameter vector $\hat{\bf w}$ according to certain criterion. Accordingly, the decision rule becomes

$$ 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.

<img src="figs/parametric_decision.png", width=300>

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}}(1|{\bf x}, {\bf w}) = g({\bf w}^\intercal{\bf x})$$$$P_{Y|{\bf,X}}(0|{\bf x}, {\bf w}) = 1-g({\bf w}^\intercal{\bf x})$$

where ${\bf w}$ is a parameter vector and $g(·)$ is the logistic function, which is defined by

$$g(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 g(t) \le 1$
• P2: Symmetry: $\quad g(-t) = 1-g(t)$
• P3: Monotonicity: $\quad g'(t) = g(t)·[1-g(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.

Exercise 2: Implement a function to compute the logistic function, and use it to plot such function in the inverval $[-6,6]$.

2.2. Classifiers based on the logistic model.

The MAP classifier under a logistic model will have the form

$$P_{Y|{\bf X}}(1|{\bf x}, {\bf w}) = g({\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 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}}(1|{\bf x}, {\bf w}) = g[{\bf w}^\intercal{\bf z}({\bf x})] $$

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

$$ {\bf w}^\intercal{\bf z} = 0 $$

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

$$ P_{Y|{\bf X}}(1|{\bf x}, {\bf w}) = g(1 + 10 x_0 + 10 x_1 - 20 x_0^2 + 5 x_0 x_1 + x_1^2) $$

3. Inference

Remember that the idea of parametric classification is to use the training data set $\mathcal S = \{({\bf x}^{(k)}, y^{(k)}) \in {\mathbb{R}}^N \times \{0,1\}, k=1,\ldots,K\}$ to set the parameter vector ${\bf w}$ according to certain criterion. Then, the estimate $\hat{\bf w}$ can be used to compute the label prediction for any new observation as

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

<img src="figs/parametric_decision.png", width=300>

In the following, 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): We assume a logistic model for the a posteriori probability of ${Y=1}$ given ${\bf X}$, i.e.,

$$P_{Y|{\bf X}}(1|{\bf x}, {\bf w}) = g[{\bf w}^\intercal{\bf z}({\bf x})].$$

We need still to choose a criterion to optimize with the selection of the parameter vector. In the notebook, we will discuss two different approaches to the estimation of ${\bf w}$:

• Maximum Likelihood (ML): $\hat{\bf w}_{\text{ML}} = \arg\max_{\bf w} P_{{\mathcal S}|{\bf W}}({\mathcal S}|{\bf w})$
• Maximum A Posteriori (MAP): $\hat{\bf w}_{\text{MAP}} = \arg\max_{\bf w} p_{{\bf W}|{\mathcal S}}({\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}}(0|{\bf x}, {\bf w}) = 1-g[{\bf w}^\intercal{\bf z}({\bf x})] = g[-{\bf w}^\intercal{\bf z}({\bf x})]$$

we can write

$$P_{Y|{\bf X}}(y|{\bf x}, {\bf w}) = g[\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

$$\hat{\bf w}_{\text{ML}} = \arg\max_{\bf w} P_{{\mathcal S}|{\bf W}}({\mathcal S}|{\bf w}) = \arg\min_{\bf w} L({\bf w}) $$

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

$$ L({\bf w}) = - \log P_{{\mathcal S}|{\bf W}}({\mathcal S}|{\bf w}) = - \log\left[P\left(y^{(1)},\ldots,y^{(K)}| {\bf x}^{(1)},\ldots, {\bf x}^{(K)},{\bf w}\right)\right] $$

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}}\left(y^{(k)}|{\bf x}^{(k)},{\bf w}\right)\right] \\ &= - \sum_{k=1}^K\log\left[P_{Y|{\bf X}}\left(y^{(k)}|{\bf x}^{(k)},{\bf w}\right)\right] \end{align}

Using A3 (the logistic model)

\begin{align} L({\bf w}) &= - \sum_{k=1}^K\log\left[g\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)}-g(\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 search for 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)}-g({\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 &= [g({\bf w}_n^\intercal {\bf z}^{(1)}), \ldots, g({\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} \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 certain conditions, the gradient descent method can be shown to converge asymptotically (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 gets very slow 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 large enough $\rho$, the gradient descent method does not converge. Can you estimate (through manual observation) an approximate value of $\rho$ stating a boundary between convergence and divergence?

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$ (but not too 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 logistic regression model can be potentially reduced by using polynomial transformations.

To compute the polynomial 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 transformation 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_{\bf w} 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)} $$

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

The numerator of the above expression is the product of two terms:

• The likelihood $P_{{\mathcal S}|{\bf W}}({\mathcal S}|{\bf w})$, which takes large values for parameter vectors $\bf w$ that fit well the training data
• The prior distribution of weights $p_{\bf W}({\bf w})$, which expresses our a priori preference for some solutions. Usually, we recur to prior distributions that take large values when $\|{\bf w}\|$ is small (associated to soft classification borders).

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}$.

Therefore, the MAP criterion prefers solutions that simultaneously fit well the data and our a priori belief about which solutions should be preferred.

$$\hat{\bf w}_{\text{MAP}} = \arg\max_{\bf w} P_{{\mathcal S}|{\bf W}}({\mathcal S}|{\bf w}) \cdot p_{\bf W}({\bf w})$$

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.

We can check that the MAP criterion adds a penalty term to the ML objective, that penalizes parameter vectors for which the prior distribution of weights takes small values.

4.1 MAP estimation with Gaussian prior

If we assume that ${\bf W}$ is a zero-mean Gaussian random variable with variance matrix $v{\bf I}$,

$$ p_{\bf W}({\bf w}) = \frac{1}{(2\pi v)^{N/2}} \exp\left(-\frac{1}{2v}\|{\bf w}\|^2\right) $$

the MAP estimate becomes

\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} \left({\bf y}-\hat{\bf p}_n\right) + \frac{2}{C}{\bf w}, $$

we obtain the following gradient descent rule for MAP estimation

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

4.2 MAP estimation with Laplacian prior

If we assume that ${\bf W}$ follows 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}$), 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.

5.2. 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$. Taking the gradient of $C({\bf w})$, and setting the result to ${\bf 0}$, 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.

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.