Introduction

Good day! 👋

In this post, the logistic regression model is covered in detail. This model is used to predict binary outcomes and is well known for its ability to drive interpretability. The objective is to start with the basics of the model formulation, consider how to interpret the model parameters, introduce the concept of how to regularise the model, give the analytical forms of the gradient vector and Hessian matrix, and define how to extend the class to the multi-class classification form.

This post is very involved. I have tried to lay out each component in a simple, consistent, and digestible format. If you spot any issues or mistakes, please reach out to me at ryanbalshaw81@gmail.com or create a merge request .

Kind regards,

Ryan Balshaw 🦮

1. The logistic regression model

Logistic regression is a foundational approach to binary classification problem used to model dichotomous outcome variables. This model defines a linear discriminative classification model, i.e., a model of $p(y\vert \mathbf{x})$ with $y \in [0, 1], \mathbf{x} \in \mathbb{R}^{d}$. To define this model, a linear mapping of data $\mathbf{x}$ to a single variable $z$, which is directly a one layer neural network mapping to a single node with a learnable bias parameter, is required. This linear mapping is given by

where $\boldsymbol{\zeta}\in \mathbb{R}^{d+1}=[\mathbf{w}^T, b]^T$ is a vector containing $\mathbf{w} \in \mathbb{R}^{d}$, a learnable weight vector, and $b$, a learnable scalar offset parameter, respectively.

Note: This is just one type of linear model, other choices exist, e.g., kernel-based methods such as radial basis functions!

To convert this linear model to a classification setting, a discriminative conditional distribution for a specific label $y$, i.e., $p(y=1\vert \mathbf{x}, \boldsymbol{\zeta})$, is constructed. This representation is given by

where $\sigma(u) = \frac{1}{1 + e^{-u}}$ represents the sigmoid function with a property that $1 - \sigma(u)=\sigma(-u)$. This conditional model is assumed to follow the following format

such that a Bernoulli distribution can be used to define an estimator for the free parameters $\widehat{\boldsymbol{\zeta}}$ via conditional maximum likelihood estimation. $\left(y_i, \mathbf{x}_i \right)$, where $i \in \left[1,N \right]$ and $y_i \in [0, 1]$, the conditional label distribution can be described as

Using this distribution, the probability mass function, a probability of a discrete random variable, for a single observation becomes

Given a dataset $\mathcal{D}$ consisting of independent and identically distributed (i.i.d.) samples from some true data distribution $p_{data} \left(\mathbf{x}, y \right)$, where $\mathbf{X} \in \mathbb{R}^{N \times d} = \left[\mathbf{x}_1, \cdots, \mathbf{x}_N \right]^T$, the likelihood function of the conditional distribution is given by

where $p_i = p \left(y=1 \vert \mathbf{x}_i, {\boldsymbol{\zeta}} \right)$. The log-likelihood (LL) function $L(\boldsymbol{\zeta}, \mathbf{X}, \mathbf{y})$ then becomes

The LL function can be used to define the conditional maximum likelihood estimator for the estimate $\widehat{\boldsymbol{\zeta}}$ in numerical, i.e., non closed-form, format

The negative LL (NLL) is commonly referred to as the cross-entropy loss due to its direct relationship to cross entropy in information theory. Moreover, for clarity, $\mathbf{x}, y$ are observed and not necessarily in our control, hence the best estimate $\widehat{\boldsymbol{\zeta}}$ is identified based on the observations captured in $\mathcal{D}$.

2. Interpreting the model formulation

The logistic regression model can be interpreted through two avenues, $(i)$ its mathematical formulation and $(ii)$ fitting the model to a simple but informative classification problem. Both avenues are explored in this section.

2.1 Mathematical interpretation

To interpret the model, some work needs to be done to figure out what linear function is being modelled. Currently, only $z(\mathbf{x}, \boldsymbol{\zeta})$ contains linearity, while $p(y = 1\vert \mathbf{x}, {\boldsymbol{\zeta}})$ is non-linear due to $\sigma(u)$. One way to interpret the linearity of the model is to try and invert $\sigma(u)$, which conveniently is given by the logit function

Thus, the logit function can be applied as follows to

This is primarily useful as it provides an indication of what is being learnt by the model: the mapping $z$ represents the log-odds between the two classes. What is the implication of representing the log-odds? Let’s consider these three scenarios, but I challenge you to think about it yourself before viewing the answers:

1. $p(y = 1\vert \mathbf{x}, {\boldsymbol{\zeta}}) = p(y = 0\vert \mathbf{x}, {\boldsymbol{\zeta}})$
2. $p(y = 1\vert \mathbf{x}, {\boldsymbol{\zeta}}) > p(y = 0\vert \mathbf{x}, {\boldsymbol{\zeta}})$
3. $p(y = 1\vert \mathbf{x}, {\boldsymbol{\zeta}}) < p(y = 0\vert \mathbf{x}, {\boldsymbol{\zeta}})$

2.2 Visual interpretation

Consider a 2D problem where we have data sampled from two Gaussians:

Training a model on this loss with $L_1$ regularisation produces the following loss function:

To interpret the logistic regression model, the log-odds and discriminative distribution for class one, i.e., $p(y=1\vert \mathbf{x}, \boldsymbol{\zeta})$ can be inspected. These two functions are given by

We can see the log-odds is a plane that splits the data into two classes.

3. Combating over-fitting via regularisation

Naively maximising the LL using gradient ascent is problematic due to over-fitting if you observe a sub-set of the true conditional distribution being approximated. There are two general solutions to this problem:

1. Penalised maximum likelihood estimation with a regularisation term $R(\boldsymbol{\zeta})$.
2. Using informative priors, $p(\boldsymbol{\zeta})$, over the parameters.

For the sake of this discussion, I will only focus on solution one as this is the process used by the paper being discussed. However, solution one and two are linked, so I will make sure that this link is clear.

The penalised maximum likelihood estimator manifests as

where $\alpha$ is a regularisation term used to control the enforcement of the parameter penalty term $R(\boldsymbol{\zeta})$. Note that the sign given to the penalisation term is negative as we are considering maximisation. Two common choices exist for $R(\boldsymbol{\zeta})$, $L_1$ regularisation and $L_2$ regularisation, otherwise known as lasso regularisation and ridge regression, respectively. An alternative is elastic net regularisation, which linearly combines $L_1$ and $L_2$ regularisation with a interpolation parameter $p$. The $L_1$ and $L_2$ regularisation terms manifest as the Euclidean and Manhattan norms

and

respectively. The choice between $R_1(\boldsymbol{\zeta})$ and $R_2(\boldsymbol{\zeta})$ comes down the (1) ease of gradient computation, and (2) parameter interpretability. $R_2(\boldsymbol{\zeta})$ has a simpler derivative, while $R_1(\boldsymbol{\zeta})$ drives weight interpretability by favouring sparsity. The Bayesian perspective of the prior $p(\boldsymbol{\zeta})$ over the parameters can identified, whereby $R_2(\boldsymbol{\zeta})$ is akin to a zero-mean, unit variance Gaussian prior distribution and $R_1(\boldsymbol{\zeta})$ is akin to a Laplace prior distribution centred around zero, which is considered sparse. Recall that Bayes theorem states

Bayes theorem represents the likelihood multiplied by the prior (the joint distribution $p(y, \boldsymbol{\zeta} \vert \mathbf{x})$) normalised by the marginal likelihood/model evidence. The form of the penalised maximum likelihood estimator can be recovered by considering the log of the numerator in isolation.

In summary, at this point in time we have a model formulation that maps from $\mathbf{x}$ to $y$ and an objective function with terms to combat overfitting. Now, all we need is data, an optimisation procedure, and code! This is the fun part, so let’s gather what we need for optimisation.

4. Estimating the model parameters $\widehat{\boldsymbol{\zeta}}$

The first choice that needs to be made, at least in principle, is how we want to solve for $\boldsymbol{\zeta}$. For those coming from an optimisation background, the decision of what to do is trivial; just maximise $L(\boldsymbol{\zeta})$ with or without regularisation. However, I want to take an alternative route, just because I can. This route is called empirical risk minimisation (ERM), which sounds quite fancy but it is simple in practice. ERM introduces the risk

which manifests as the expectation over the negative of the log of the conditional label distribution. ERM approximates the risk by the empirical risk given a sample set $\mathcal{D}$, the observed dataset, as

which is just the normalised negative log-likelihood (NLL). This empirical risk function is then minimised, hence the name: ERM. To minimise this function, often gradient-based optimisation is used. This is a more “statistical” approach to obtain an optimal estimate $\widehat{\boldsymbol{\zeta}}$, but in actual fact this is just a round-about way to get to saying: just minimise the sample-normalised NLL function. However, I find it provides a suitable approach to reproduce the terms typically seen in other derivations of the logistic regression model. A standard ERM procedure is mini-batch gradient descent, which uses a batch of data $\mathbf{D}_b$, a subset of $\mathcal{D}$. Assuming $\mathcal{L}(\boldsymbol{\zeta})$ is differentiable, we can identify different terms for different optimisation procedures. The first of these is the gradient vector.

Note that while I have not included the regularisation in this formulation, you can add it in by adding it to $\mathcal{L}(\boldsymbol{\zeta})$: $\mathcal{L}_R(\boldsymbol{\zeta}) = \mathcal{L}(\boldsymbol{\zeta}) + \alpha R(\boldsymbol{\zeta})$. If the decision is made to use gradient descent, a first-order iterative method, then the constraint is that $R(\boldsymbol{\zeta})$ must be differentiable. Otherwise,

4.1 The gradient vector

The gradient vector $\nabla_{\boldsymbol{\zeta}} \mathcal{L}(\boldsymbol{\zeta})$ is given by

where $\mathbf{y} \in \mathbb{R}^N$ is a vector of the expected class labels, $\mathbf{s}=\sigma(\mathbf{X}\boldsymbol{\zeta})\in \mathbb{R}^N$, and $\overline{\mathbf{X}} \in \mathbb{R}^{N \times d + 1}$ is represented as

Note: This is the gradient vector of $\mathcal{L}(\boldsymbol{\zeta})$ and it does not include the regularisation term. To include it, you will just need to add $\alpha \cdot \nabla_{\boldsymbol{\zeta}}R(\boldsymbol{\zeta})$ to $\nabla_{\boldsymbol{\zeta}} \mathcal{L}(\boldsymbol{\zeta})$.

4.2 The Hessian

Certain solvers, e.g., second-order methods such as Newton’s method in Optimisation, require Hessian information to perform parameter updates. The Hessian is expressed as the transpose of the Jacobian of the gradient vector, $\mathbf{H}(f(\mathbf{x})) = \mathbf{J}\left( \nabla f(\mathbf{x}) \right)^T.$

The Hessian matrix $\mathbf{H}(\boldsymbol{\zeta})$ is given by

Note: This is the Hessian matrix of $\mathcal{L}(\boldsymbol{\zeta})$ and it does not include the regularisation term. To include it, you will just need to add $\alpha \cdot \mathbf{J}\left( \nabla R(\boldsymbol{\zeta}) \right)^T$ to $\mathbf{H}(\boldsymbol{\zeta})$.

5. Implementing this in code

I realise that this post is growing in length, so I recommend that you look at my Python-based implementation inside the sjmelck_pages repository itself or inside the random_paper_implementations repository. In both repos, I constructed the model to be compatible with scikit-learn, just for the challenge. I am fully aware that it already exists within scikit-learn’s ecosystem.

6. Multinomial logistic regression

To wrap this discussion up, I will also discuss how multinomial logistic regression can be scaled up to multiple classes. I will not progress further than the model formulation and how to compute one of the more complex terms in the gradient vector. The first step is to introduce a mapping to $\mathbf{z} \in \mathcal{R}^K$, where $K$ refers to the number of classes considered. This is achieved using

where $\mathbf{W}=[\mathbf{w}_k, \cdots, \mathbf{w}_K]^T$ represents $K$ independent mappings and $\mathbf{b} = [b_k, \cdots, b_K]^T$. To derive the form of $p(y_i=k\vert \mathbf{x})$, there are two prominent forms in the literature, namely via $(i)$ minimal parametrisation and $(ii)$ over-parametrisation of the weight vectors. I will discuss each in turn.

6.1 Minimal parametrisation

To derive this formation, the constraint $\sum_{l=1}^{K}p(y=l \vert \mathbf{x})=1$ is key. Specifically, we can write

and then parametrise the log-odds between $p(y=k \vert \mathbf{x})$ and $p(y=K \vert \mathbf{x})$ with the linear mapping to $z_k=\mathbf{w}_k^T \mathbf{x} + b_k$. This is given as

for the log-odds and as

for the odds ratio. Hence, under this formulation, $p(y=k \vert \mathbf{x}) = e^{\mathbf{w}_k^T \mathbf{x} + b_k}p(y=K \vert \mathbf{x})$. Substituting this into the probability sum constraint returns

which can then be used to obtain

In this solution, the fact that one classes probability is fully parametrised by the probabilities of the other classes is key.

6.2 Over-parametrisation

In the over-parametrised form, the Gibbs measure is used to define

where $Z$ is a partition function that can be used to enforce that the Gibbs measure outputs a valid probability. Specifically, the sum over discrete random variables constraint can be used to define $Z$

Hence, the conditional distribution $p(y=k \vert \mathbf{x})$ is given by

where $\boldsymbol{\sigma}_s(\mathbf{z}): \mathbb{R}^K \rightarrow \mathbb{R}^K$ is commonly referred to as the softmax function. One problem with this term is that it is invariant to translation, i.e., $\sigma_s^{(k)}(\mathbf{z} + \mathbf{c}) = \sigma_s^{(k)}(\mathbf{z})$, and if we set $\mathbf{c} = \boldsymbol{\zeta}_k$, we can recover the minimal parametrisation form, effectively removing the weights for the $K^{th}$ class.

6.3 The multinomial objective function

For both forms of $p(y=k\vert \mathbf{x})$, the observed class labels $y_i \in {1, \cdots, K}$ for $i=1, \cdots, N$ are considered as samples of categorically distributed random variables $Y_1, \cdots, Y_K$ leading to the probability mass function

where $p_{ij}=p(y=j \vert \mathbf{x}_i)$, $y_i$ is the class label for sample $x_i$, and $[y=k]$ is the Iverson bracket

The likelihood function is defined by

One interesting take-away is that this loss only considers the probability given to the expected class label, and on the surface it appears that the parameter updates will only care about the correct class prediction. To determine if this is truly the case, it is of interest to look at how the parameters are updated. Specifically, iterative parameter updates, e.g., gradient descent, require the gradient vector, which now relies on the term $\partial \boldsymbol{\sigma} / \partial \mathbf{z}$, i.e., the Jacobian of the softmax function with respect to $\mathbf{z}$. The form of this term is discussed to drive intuition behind the loss.

6.4 The Jacobian of the softmax function

Each row in the Jacobian $\mathbf{J}_\boldsymbol{\sigma} , \boldsymbol{\sigma}(\mathbf{z}) \in \mathbb{R}^{K \times K}$ is given by

For simplicity, lets consider one term in this vector $\frac{\partial \sigma_s^{(i)}(\mathbf{z})}{\partial z_j}$, which represents $J_{\boldsymbol{\sigma}, ij} , \sigma(\mathbf{z})$. To help this process, the log-derivative trick can be used, which is given by

Hence, this trick yields $\frac{\partial f}{\partial u}= f \frac{\partial}{\partial u} \log f$. Applying this to the considered term yields

Hence, the full Jacobian can be written neatly as

This relation indicates that the gradient vector will contain sensitivity information on both the correct class predictions and the incorrect class predictions.

7. Conclusion

If you get this far, thank you for reading this post! I hope you found parts of this write-up informative and that you ended off with a better understanding than when you started reading the post! If you did not, please let me know!