Fun with Logistic Regression¶

This notebook was derived from the Caltech "Learning from Data" course, specifically Lecture 9 on the logistic regression model. It is a simple Python implementation of the logistic regression model using NumPy.

Table of Contents¶

• Derivation
  • Error Measure
  • Learning Algorithm
• Example
  • Data Set
  • Learning a Hypothesis
  • Learning with Less Data
  • Learning Curve
  • Visualizing the Error Surface
  • Nonlinear Transformation

Notation¶

The problem structure is the classic classification problem. Our data set $\mathcal{D}$ is composed of $N$ samples. Each sample is a tuple containing a feature vector and a label. For any sample $n$ the feature vector is a $d+1$ dimensional column vector denoted by ${\bf x}_n$ with $d$ real-valued components known as features. Samples are represented in homogeneous form with the first component equal to $1$: $x_0=1$. Vectors are bold-faced. The associated label is denoted $y_n$ and can take on only two values: $+1$ or $-1$.

$$\mathcal{D} = \lbrace ({\bf x}_1, y_1), ({\bf x}_2, y_2), ..., ({\bf x}_N, y_N) \rbrace \\ {\bf x}_n = \begin{bmatrix} 1 & x_1 & ... & x_d \end{bmatrix}^T$$

Derivation¶

Despite the name logistic regression this is actually a probabilistic classification model. It is also a linear model which can be subjected to nonlinear transforms.

All linear models make use of a "signal" $s$ which is a linear combination of the input vector ${\bf x}$ components weighed by the corresponding components in a weight vector ${\bf w}$.

$${\bf w} = \begin{bmatrix} w_0 & w_1 & ... & w_d \end{bmatrix}^T \\ s = w_0 + w_1 x_1 + \;...\; + w_d x_d = \sum_{i=0}^d w_i x_i = {\bf w} \cdot {\bf x} = {\bf w}^T {\bf x}$$

Note that the homogeneous representation (with the $1$ at the first component) allows us to include a constant offset using a more compact vector-only notation (instead of ${\bf w}^T {\bf x}+b$).

Linear classification passes the signal through a harsh threshold:

$$h({\bf x}) = \operatorname{sign}(s)$$

Linear regression uses the signal directly without modification:

$$h({\bf x}) = s$$

Logistic regression passes the signal through the logistic/sigmoid but then treats the result as a probability:

$$h({\bf x}) = \theta(s)$$

The logistic function is

$$\theta(s) = \frac{e^s}{1+e^s} = \frac{1}{1+e^{-s}}$$

There are many other formulas that can achieve a soft threshold such as the hyperbolic tangent, but this function results in some nice simplification.

We say that the data is generated by a noisy target.

$$P(y\mid{\bf x})=\begin{cases} f({\bf x}) & \text{for }y=+1 \\ 1-f({\bf x}) & \text{for }y=-1 \end{cases}$$

With this noisy target we want to learn a hypothesis $h({\bf x})$ that best fits the above noisy target according to some error function.

$$h({\bf x})=\theta({\bf w}^T {\bf x})\approx f({\bf x})$$

It's important to note that the data does not tell you the probability of a label but rather what label the sample has after being generated by the target distribution.

Error Measure¶

To learn a good hypothesis we want to find a hypothesis parameterization ${\bf w}$ (the weight vector) that minimizes some in-sample error measure $E_\text{in}$.

$${\bf w}_h = \underset{{\bf w}}{\operatorname{argmin}} \; E_\text{in}({\bf w})$$

The error measure we will use is both plausible and nice. It is based on likelihood which is the probability of generating the data given a model.

If our hypothesis is close to our target distribution ($h\approx f$) then we expect that probability of generating the data to be high.

There is some controversy with using likelihood. We are really looking for the most probable hypothesis given the data: $\underset{h}{\operatorname{argmax}} P(h\mid{\bf x})$. The likelihood approach is looking for the hypothesis that makes the data most probable: $\underset{h}{\operatorname{argmax}} P({\bf x}\mid h)$.

The Bayesian approach tackles this issue using Bayes' Theorem but introduces other issues such as choosing priors.

To determine the likelihood we assume the data was generated with our hypothesis $h$:

$$P(y\mid{\bf x})=\begin{cases} h({\bf x}) & \text{for }y=+1 \\ 1-h({\bf x}) & \text{for }y=-1 \end{cases} \\ $$

where $h({\bf x})=\theta({\bf w}^T {\bf x})$.

We don't want to deal with cases so we take advantage of a nice property of the logistic function: $\theta(-s)=1-\theta(s)$.

$$\text{if } y = +1 \text{ then } h({\bf x}) = \theta({\bf w}^T {\bf x}) = \theta(y \; {\bf w}^T {\bf x}) \\ \text{if } y = -1 \text{ then } 1 - h({\bf x}) = 1 - \theta({\bf w}^T {\bf x}) = \theta(- {\bf w}^T {\bf x}) = \theta(y \; {\bf w}^T {\bf x}) \\ $$

Using this simplification,

$$P(y\mid{\bf x})=\theta(y\; {\bf w}^T {\bf x})$$

The likelihood is defined for a data set $\mathcal{D}$ with $N$ samples given a hypothesis (denoted arbitrarily $g$ here):

$$L(\mathcal{D} \mid g) = \prod_{n=1}^{N} P(y_n \mid {\bf x}_n) = \prod_{n=1}^{N} \theta(y_n \; {\bf w}_g^T {\bf x}_n)$$

Now finding a good hypothesis is a matter of finding a hypothesis parameterization ${\bf w}$ that maximizes the likelihood.

$${\bf w}_h = \underset{{\bf w}}{\operatorname{argmax}} \; L(\mathcal{D} \mid h) = \underset{{\bf w}}{\operatorname{argmax}} \; \theta(y_n \; {\bf w}^T {\bf x}_n)$$

Maximizing the likelihood is equivalent to maximizing the log of the function since the natural logarithm is a monotonically increasing function:

$$\underset{{\bf w}}{\operatorname{argmax}} \; \ln \left( \prod_{n=1}^{N} \theta(y_n \; {\bf w}^T {\bf x}_n) \right)$$

We can maximize the above proportional to a constant as well so we'll tack on a $\frac{1}{N}$:

$$\underset{{\bf w}}{\operatorname{argmax}} \; \frac{1}{N} \ln \left( \prod_{n=1}^{N} \theta(y_n \; {\bf w}^T {\bf x}_n) \right)$$

Now maximizing that is the same as minimizing its negative:

$$\underset{{\bf w}}{\operatorname{argmin}} \left[ -\frac{1}{N} \ln \left( \prod_{n=1}^{N} \theta(y_n \; {\bf w}^T {\bf x}_n) \right) \right]$$

If we move the negative into the log and the log into the product we turn the product into a sum of logs:

$$\underset{{\bf w}}{\operatorname{argmin}} \;\frac{1}{N} \sum_{n=1}^{N} \ln \left( \frac{1}{\theta(y_n \; {\bf w}^T {\bf x}_n)} \right)$$

Expanding the logistic function,

$$\underset{{\bf w}}{\operatorname{argmin}} \;\frac{1}{N} \sum_{n=1}^{N} \ln \left( 1 + e^{y_n \; {\bf w}^T {\bf x}_n} \right)$$

Now we have a much nicer form for the error measure known as the "cross-entropy" error.

$$E_\text{in}({\bf w}) = \frac{1}{N} \sum_{n=1}^{N} \ln \left( 1+e^{-y_n \; {\bf w}^T {\bf x}_n} \right)$$

This is nice because it can be interpreted as the average point error where the point error function is

$$e(h({\bf x}_n), y_n) = \ln \left( 1+e^{-y_n \; {\bf w}^T {\bf x}_n} \right) \\ E_\text{in}({\bf w}) = \frac{1}{N} \sum_{n=1}^{N} e(h({\bf x}_n), y_n)$$

So to learn a hypothesis we'll want to perform the following optimization:

$${\bf w}_h = \underset{{\bf w}}{\operatorname{argmin}} \; E_\text{in}({\bf w}) = \underset{{\bf w}}{\operatorname{argmin}} \;\frac{1}{N} \sum_{n=1}^{N} \ln \left( 1 + e^{y_n \; {\bf w}^T {\bf x}_n} \right)$$

Learning Algorithm¶

The learning algorithm is how we search the set of possible hypotheses (hypothesis space $\mathcal{H}$) for the best parameterization (in this case the weight vector ${\bf w}$). This search is an optimization problem looking for the hypothesis that optimizes an error measure.

There is no nice, closed-form solution like with least-squares linear regression so we will use gradient descent instead. Specifically we will use batch gradient descent which calculates the gradient from all data points in the data set.

Luckily, our "cross-entropy" error measure is convex so there is only one minimum. Thus the minimum we arrive at is the global minimum.

Gradient descent is a general method and requires twice differentiability for smoothness. It updates the parameters using a first-order approximation of the error surface.

$${\bf w}_{i+1} = {\bf w}_i + \nabla E_\text{in}({\bf w}_i)$$

To learn we're going to minimize the following error measure using batch gradient descent.

$$e(h({\bf x}_n), y_n) = \ln \left( 1+e^{-y_n \; {\bf w}^T {\bf x}_n} \right) \\ E_\text{in}({\bf w}) = \frac{1}{N} \sum_{n=1}^{N} e(h({\bf x}_n), y_n) = \frac{1}{N} \sum_{n=1}^{N} \ln \left( 1+e^{-y_n \; {\bf w}^T {\bf x}_n} \right)$$

We'll need to the derivative of the point loss function and possibly some abuse of notation.

$$\frac{d}{d{\bf w}} e(h({\bf x}_n), y_n) = \frac{-y_n \; {\bf x}_n \; e^{-y_n {\bf w}^T {\bf x}_n}}{1 + e^{-y_n {\bf w}^T {\bf x}_n}} = -\frac{y_n \; {\bf x}_n}{1 + e^{y_n {\bf w}^T {\bf x}_n}}$$

With the point loss derivative we can determine the gradient of the in-sample error:

$$\begin{align} \nabla E_\text{in}({\bf w}) &= \frac{d}{d{\bf w}} \left[ \frac{1}{N} \sum_{n=1}^N e(h({\bf x}_n), y_n) \right] \\ &= \frac{1}{N} \sum_{n=1}^N \frac{d}{d{\bf w}} e(h({\bf x}_n), y_n) \\ &= \frac{1}{N} \sum_{n=1}^N \left( - \frac{y_n \; {\bf x}_n}{1 + e^{y_n {\bf w}^T {\bf x}_n}} \right) \\ &= - \frac{1}{N} \sum_{n=1}^N \frac{y_n \; {\bf x}_n}{1 + e^{y_n {\bf w}^T {\bf x}_n}} \\ \end{align}$$

Our weight update rule per batch gradient descent becomes

$$\begin{align} {\bf w}_{i+1} &= {\bf w}_i - \eta \; \nabla E_\text{in}({\bf w}_i) \\ &= {\bf w}_i - \eta \; \left( - \frac{1}{N} \sum_{n=1}^N \frac{y_n \; {\bf x}_n}{1 + e^{y_n {\bf w}_i^T {\bf x}_n}} \right) \\ &= {\bf w}_i + \eta \; \left( \frac{1}{N} \sum_{n=1}^N \frac{y_n \; {\bf x}_n}{1 + e^{y_n {\bf w}_i^T {\bf x}_n}} \right) \\ \end{align}$$

where $\eta$ is our learning rate.

Example¶

First we'll get the dependencies out of the way.

For pretty visualization we'll define some custom color maps.

We'll define our logistic/sigmoid function using the latter definition since it uses only one exponentiation.

$$\theta(s)=\frac{e^s}{1+e^s}=\frac{1}{1+e^{-s}}$$

Data Set¶

For this example we'll create a toy data set.

To visualize this properly we'll set the input dimension size to $d=2$. For now we'll keep it linear ($\phi({\bf x})={\bf x}$).

For our data set we'll generate $N=100$ points.

We'll define our input distribution $P({\bf x})$ to be uniform within the intervals $\left[-1, 1\right]$ for visualization purposes. Note the homogeneous representation with $x_0=1$.

Now we need to create a fake, linear, noisy target $P(y\mid{\bf x})$. The target distribution is parameterized by a "target function" $f$ which will be drawn randomly from our hypothesis set (linear logistic probability models). The weights are randomly chosen for good visual effect (thus the hardness measure).

With our target distribution we can generate our toy data set. Note that $\mathcal{Z}$ is the same as $\mathcal{X}$ in this linear example.

Now let's visualize our data set. For our pretty fills and contours we need to define the entire grid to some level of detail. Here we'll use $s=300$ meaning 300 points per axis.

We'll plot our data set along with the background probability from the noisy target. Red points are classified as $+1$ and blue points are classified as $-1$. The background color denotes the probability of the data generating distribution. The black line is the decision boundary where the probability is equal to $0.5$.

Learning a Hypothesis¶

Let's initialize our hypothesis parameters to zero and perform gradient descent with a learning rate of $\eta=10.0$. We'll continue iteration for a maximum of $\text{10,000}$ iterations or until the change in the weight vector's norm is smaller than $0.1\%$, which ever comes first. As a reminder, our update rule is

$${\bf w}_{i+1} = {\bf w}_i + \eta \; \left( \frac{1}{N} \sum_{n=1}^N \frac{y_n \; {\bf x}_n}{1 + e^{y_n {\bf w}_i^T {\bf x}_n}} \right)$$

Now we'll visualize the final hypothesis $g$ in the input space $\mathcal{X}$.

Since we have the target distribution we can stochastically estimate the out-of-sample error. Classifying using the logistic regression is a matter of choosing the label that is most probable for a data point. So if the output of the hypothesis function $h$ is greater than $0.5$ we label the data point $y=+1$ and vice versa.

Note the relatively high error rates for a toy data set due to treating the decision boundary as a hard threshold despite the noisy nature of the target distribution. The higher our hardness factor on our toy data set the smaller this error rate will be.

Learning with Less Data¶

If we learn from a smaller set of data, say $N=10$, we should see a larger discrepancy in the out-of-sample error.

If the resulting hypothesis looks off we need to remember what the learning algorithm sees: just the data.

With the metrics below we can see that the quality of the final hypothesis depends on the number of samples we have as predicted by learning theory.

Learning Curve¶

Since we can stochastically estimate the out-of-sample error we can plot approximate learning curves. The learning curve will show us how the in-sample error and out-of-sample error behave as the data resources improve.

Warning: this can take a while to evaluate.

Visualizing the Error Surface¶

For curiosity sake, let's try and visualize the error surface. But before visualizing the surface let's look at how the error measure changes as the gradient descent proceeds.

The in-sample error is parameterized by our model parameters ${\bf w}$. Unfortunately our parameter is of dimension $d_{\bf w}=3$ which makes visualization a bit more difficult (a four dimensional visualization problem: three parameters and an error value). Instead of tackling that directly (volumetrics and other fun stuff) we'll just visualize the top, front, and left slices of the error surface centered on the final hypothesis parameters. Along with these slices the iterations we've taken on the error surface using gradient descent which will be colored by iteration number backwards through a rainbow color map. Using the history of the weight vector through the gradient descent iterations we can determine the extents of our visualization.

Nonlinear Transformation¶

Like with the other linear models (perceptron, linear regression) we can transform the input data from the original input space ($\mathcal{X}$) into a higher dimension space ($\mathcal{Z}$) and fit our model there.

In this example we'll use a second-order polynomial transformation:

$$\phi({\bf x}) = \begin{bmatrix} 1 & x_0 & x_1 & x_0 x_1 & x_0^2 & x_1^2 \end{bmatrix}$$

For fun we'll try visualizing the error surface again. However, since we're operating in the $\mathcal{Z}$ space with dimension $d_\mathcal{Z}=6$ our two-dimensional visualizations are even more woefully inadequate. Regardless, below is the gradient descent evolution and all 15 combinations of two axes and their two-dimensional error "surface" slice centered on the hypothesis weight vector.