# Part 2: Logistic Regression and Gradient Checking


```
# Execute this code block to install dependencies when running on colab
try:
    import torch
except:
    from os.path import exists
    from wheel.pep425tags import get_abbr_impl, get_impl_ver, get_abi_tag
    platform = '{}{}-{}'.format(get_abbr_impl(), get_impl_ver(), get_abi_tag())
    cuda_output = !ldconfig -p|grep cudart.so|sed -e 's/.*\.\([0-9]*\)\.\([0-9]*\)$/cu\1\2/'
    accelerator = cuda_output[0] if exists('/dev/nvidia0') else 'cpu'

    !pip install -q http://download.pytorch.org/whl/{accelerator}/torch-1.0.0-{platform}-linux_x86_64.whl torchvision
```

In the first part of the lab we saw how to make predictions of continously varying values with a linear regression model. Lets now turn our focus to binary classification using a simple classification algorithm known as Logistic regression.

In linear regression we tried to predict the value of $y$ for an example $\mathbf{x}$ using a linear function $y=\mathbf{x}^\top\theta$ (where $\mathbf{x}$ and $\theta$ are column-vectors). This will clearly not be a great solution for predicting binary-valued labels ($y\in\{0,1\}$). In logistic regression we use a different hypothesis class to try to predict the probability that a given example belongs to the "1" class versus the probability that it belongs to the "0" class. Specifically, we will try to learn a function of the form:

\begin{align}
P(y=1|\mathbf{x}) &= \frac{1}{1 + \exp(-\mathbf{x}^\top\theta)} \equiv \sigma(\mathbf{x}^\top\theta),\\
P(y=0|\mathbf{x}) &= 1 - P(y=1|\mathbf{x}).
\end{align}
 
The function $\sigma(z) \equiv \frac{1}{1 + \exp(-z)}$ is called the "sigmoid" or "logistic" function. The sigmoid function squashes any real valued input into the range $[0,1]$ enabling us to interprete the output as a probability. Our goal is to search for a value of $\theta$ so that the probability $P(y=1|\mathbf{x})=\sigma(\mathbf{x}^\top\theta)$ is large when $\mathbf{x}$ belongs to the "1" class and small when $\mathbf{x}$ belongs to the "0" class (so that $P(y=0|\mathbf{x})$ is large). 

With Linear Regression, the natural cost function was one that measured the sum of squared residuals (the difference between the predicted value and true value). With logisitic regression we have a probabilisitic model, so it makes sense that we use a function that measures the likelihood of the data given the model (note that we want to maximise this function rather than minimise it). As an aside, note that in the case of linear regression if we assume that the data has errors that are IID (independently and identically distributed) according to a Normal distribution, then it can be shown that the maximising the likelihood is exactly the same as minimising the sum of squared residuals. For logistic regression, the likelihood function for a single data point is:

\begin{align}
p(y|\mathbf{x}; \theta) &= \sigma(\mathbf{x}^\top\theta)^y(1-\sigma(\mathbf{x}^\top\theta)^{(1-y)}.
\end{align}

For a complete dataset of points $(y_i, \mathbf{x}_i)$, then the complete likelihood is:

\begin{align}
L(\theta) &= \prod_i \sigma(\mathbf{x}_i^\top\theta)^{y_i}(1-\sigma(\mathbf{x}_i^\top\theta)^{(1-y_i)}
\end{align}

However, it is considerably easier to maximise the log-likelihood function:

\begin{align}
\mathcal{l}(\theta) &= \log L(\theta) \\
                    & = \log \prod_i \sigma(\mathbf{x}_i^\top\theta)^{y_i}(1-\sigma(\mathbf{x}_i^\top\theta)^{(1-y_i)} \\
                    & = \sum_i y_i \log(\sigma(\mathbf{x}_i^\top\theta)) + (1-y_i) \log(1-\sigma(\mathbf{x}_i^\top\theta))
\end{align}

Clearly, maximising the log-likelihood is equivalent to minimising the negative log-likelihood. The negative of the log-likelihood function having the form $-\sum_i y_i \log(p) + (1-y_i) \log(p)$, where p is a function returning the predicted probability of class "1", is often called the __"Binary Cross Entropy"__ function, __"Binary Cross Entropy Loss"__ or sometimes the __"log loss"__.

For conciseness and computational efficiency, we can write the negative logistic regression log-likelihood function in matrix form. Assuming the $y_i$ are stored in a column vector $\mathbf{y}$ and the data vectors $x_i$ in the __rows__ of a matrix $\mathbf{X}$, then: 

\begin{align}
\mathrm{NLL}(\theta) & = -(\mathbf{y} \log(\sigma(\mathbf{X}\theta)) + (1-\mathbf{y}) \log(1-\sigma(\mathbf{X}\theta)))
\end{align}

The gradients of this function are given by:

\begin{align}
\nabla_\theta \mathrm{NLL}(\theta) & = \mathbf{X}^T(\sigma(\mathbf{X}\theta) - \mathbf{y})
\end{align}

__Use the box below to compute the gradients of the negative log-likelihood function $\nabla_\theta \mathrm{NLL}(\theta)$. You can use `torch.sigmoid()` to apply the sigmoid function.__


```
import torch

# we wouldn't normally do this, but for this lab we want to work in double precision
# as we'll need the numerical accuracy later on for doing checks on our gradients:
torch.set_default_dtype(torch.float64) 

def logistic_regression_loss_grad(theta, X, y):
    # YOUR CODE HERE
    grad=X.t()@(torch.sigmoid(X@theta)-y)
#     raise NotImplementedError()
    return grad
```


```
theta = torch.zeros(1)
X = torch.Tensor([[1]])
y = torch.Tensor([[0]])
assert(logistic_regression_loss_grad(theta, X, y) == 0.5)

```

## Training a Logistic Regressor with real data

We'll now try gradient descent using our gradient function on a real dataset from `scikit-learn` called `digits`. 

The `digits` dataset contains handwritten characters (much like the `MNIST` dataset that you may have heard of - we'll explore `MNIST` in a future lab). As logistic regression is a binary classifier, we'll just use the first 2 characters (0 and 1) from the `digits` dataset, and make our own training and test splits:


```
from sklearn.datasets import load_digits

X, y = tuple(torch.Tensor(z) for z in load_digits(2, True)) #convert to pytorch Tensors
X = torch.cat((X, torch.ones((X.shape[0], 1))), 1) # append a column of 1's to the X's
y = y.reshape(-1, 1) # reshape y into a column vector

# We're also going to break the data into a training set for computing the regression parameters
# and a test set to evaluate the predictive ability of those parameters
perm = torch.randperm(y.shape[0])
X_train = X[perm[0:260], :]
y_train = y[perm[0:260]]
X_test = X[perm[260:], :]
y_test = y[perm[260:]]
```

Now we have the data, we can use our loss function to try and estimate the optimal parameters for the two-digit classification problem. We'll use `PyTorch`s `torch.nn.functional.binary_cross_entropy_with_logits` function to print out the Binary Cross Entropy of the training data at each iteration, and of the test data once the optimisation is complete. 

Note: `logits` refers to unscaled probabilities before the sigmoid is applied, so in the `binary_cross_entropy_with_logits` function we just pass in $\mathbf{X}\theta$. `PyTorch` does also have a `torch.nn.binary_cross_entropy` method that takes in probabilities, however, as we'll see when implementing neural networks in a later lab, we'll often choose to work with logits as they provide better numerical stability thanks to the _log-sum-exp_ trick. 


```
alpha = 0.001
theta_gd = torch.rand((X_train.shape[1], 1))
for e in range(0, 10):
    gr = logistic_regression_loss_grad(theta_gd, X_train, y_train)
    theta_gd -= alpha * gr
    print("Epoch:", e, " BCE of training data:", torch.nn.functional.binary_cross_entropy_with_logits(X_train @ theta_gd, y_train))

print("Gradient Descent Theta:", theta_gd.t())
print("BCE of test data:", torch.nn.functional.binary_cross_entropy_with_logits(X_test @ theta_gd, y_test))
```

    Epoch: 0  BCE of training data: tensor(82.0769)
    Epoch: 1  BCE of training data: tensor(19.0741)
    Epoch: 2  BCE of training data: tensor(32.0213)
    Epoch: 3  BCE of training data: tensor(0.2375)
    Epoch: 4  BCE of training data: tensor(0.1180)
    Epoch: 5  BCE of training data: tensor(0.0898)
    Epoch: 6  BCE of training data: tensor(0.0755)
    Epoch: 7  BCE of training data: tensor(0.0612)
    Epoch: 8  BCE of training data: tensor(0.0469)
    Epoch: 9  BCE of training data: tensor(0.0326)
    Gradient Descent Theta: tensor([[ 0.6219,  0.4870, -0.4789, -0.6284,  0.0851,  0.9381,  0.2797,  0.7038,
              0.6723,  0.2665, -2.1668, -0.3051,  0.7746,  0.1612,  0.4266,  0.4605,
              0.8650, -0.6009, -2.2096,  2.2712,  2.9163, -1.2616, -0.2625,  0.0411,
              0.3620, -0.8273, -0.9611,  2.9302,  3.5129, -0.0618, -0.9063,  0.4057,
              0.5006, -1.1616, -0.9266,  3.0836,  3.2987, -0.5377, -1.5661,  0.1528,
              0.1053,  0.0253, -1.5991,  1.8560,  3.0142, -1.3664, -0.4140,  0.6402,
              0.9595,  0.6807, -1.4236, -0.1138,  0.4964, -0.9779,  0.0708,  0.5491,
              0.4488,  0.2829,  0.1307, -0.9637,  0.3821,  1.0497,  0.9260,  0.8819,
              0.8054]])
    BCE of test data: tensor(2.4425e-17)


What do you observe from running the above? Write your answer in the following block:

- **Answer**:

**With the increasing of the  number of iterations, the Binary Cross Entropy of the data will go down even is quite near to zero, meaning the predicted probability of class "1" is near to 100% which is the maximum. **

## Gradient Checking

How can we be sure that our gradient function was correct? We might have made an error in the analytical derivation or in the implementation of that derivation into code. Even though we observed the optimisation process on real data converge (hopefully!), you might have made a subtle error in the implementation...

So far we have worked with relatively simple algorithms where it is straightforward to compute the objective function and its gradient with pen-and-paper, and then implement the necessary computations in PyTorch. For more complex models that we will see later, the gradient computation can be notoriously difficult to debug and get right. Sometimes a subtly buggy implementation will manage to learn something that can look surprisingly reasonable (while performing less well than a correct implementation). Thus, even with a buggy implementation, it may not at all be apparent that anything is amiss. In this section, we describe a method for numerically checking the derivatives computed by your code to make sure that your implementation is correct. Carrying out the derivative checking procedure described here will significantly increase your confidence in the correctness of your code.

Suppose we want to minimize $J(\theta)$ as a function of $\theta$. For this example, suppose $J:\mathbb{R}\mapsto\mathbb{R}$, so that $\theta \in \mathbb{R}$. If we are using gradient descent (or other gradient-based optimisation function), then we usually have implemented some function $g(\theta)$ that purportedly computes $\frac{d}{d\theta}J(\theta)$.

How can we check if our implementation of $g$ is correct?

Recall the mathematical definition of the derivative is:

\begin{align}
\frac{d}{d\theta}J(\theta) = \lim_{\epsilon \rightarrow 0}
\frac{J(\theta+ \epsilon) - J(\theta-\epsilon)}{2 \epsilon}.
\end{align}

Thus, at any specific value of $\theta$, we can numerically approximate the derivative as follows:

\begin{align}
\frac{J(\theta+{\rm EPSILON}) - J(\theta-{\rm EPSILON})}{2 \times {\rm EPSILON}}
\end{align}
 
In practice, we set ${\rm EPSILON}$ to a small constant, say around $10^{−4}$. (There is a large range of values of EPSILON values that should work well, but we don’t set ${\rm EPSILON}$ to be "extremely" small, say $10^{−20}$, as that would lead to numerical roundoff errors.)

Thus, given a function $g(\theta)$ that is supposedly computing $\frac{d}{d\theta}J(\theta)$, we can now numerically verify its correctness by checking that

\begin{align}
g(\theta) \approx
\frac{J(\theta+{\rm EPSILON}) - J(\theta-{\rm EPSILON})}{2 \times {\rm EPSILON}}.
\end{align}
 
The degree to which these two values should approximate each other will depend on the details of $J$. But assuming ${\rm EPSILON}=10^{−4}$, you’ll usually find that the left- and right-hand sides of the above will agree to at least 4 significant digits (and often many more).

Now, consider the case where $\theta \in \mathbb{R}^n$ is a vector rather than a single real number (so that we have $n$ parameters that we want to learn), and $J: \mathbb{R}^n \mapsto \mathbb{R}$. We now generalize our derivative checking procedure to the case where $\theta$ may be a vector (as in our linear regression and logistic regression examples).

Suppose we have a function $g_i(\theta)$ that purportedly computes $\frac{\partial}{\partial\theta_i}J(\theta)$; we’d like to check if $g_i$ is outputting correct derivative values. Let $\textstyle \theta^{(i+)} = \theta + {\rm EPSILON} \times \vec{e}_i$, where

\begin{align}
\vec{e}_i = \begin{bmatrix}0 \\ 0 \\ \vdots \\ 1 \\ \vdots \\ 0\end{bmatrix}
\end{align}

is the $i$-th basis vector (a vector of the same dimension as $\theta$, with a "1" in the $i$-th position and "0"s everywhere else). So, $\theta^{(i+)}$ is the same as $\theta$, except its $i$-th element has been incremented by ${\rm EPSILON}$. Similarly, let $\theta^{(i−)}=\theta−{\rm EPSILON} \times \vec{e}_i$ be the corresponding vector with the $i$-th element decreased by ${\rm EPSILON}$.

We can now numerically verify $g_i(\theta)$'s correctness by checking, for each $i$, that:

\begin{align}
g_i(\theta) \approx
\frac{J(\theta^{(i+)}) - J(\theta^{(i-)})}{2 \times {\rm EPSILON}}.
\end{align}

### Gradient checker code

The following code block contains an implementation of the gradient checking proceedure described above. 


```
from random import randrange

def grad_check(f, x, analytic_grad, num_checks=10, h=1e-5):
    sum_error = 0
    for i in range(num_checks):
        ix = tuple([randrange(m) for m in x.shape]) #randomly sample value to change
        print("ix:",ix)
        oldval = x[ix].item()
        print("oldval:",oldval)
        x[ix] = oldval + h # increment by h
        fxph = f(x) # evaluate f(x + h)
        x[ix] = oldval - h # increment by h
        fxmh = f(x) # evaluate f(x - h)
        x[ix] = oldval # reset

        grad_numerical = (fxph - fxmh) / (2 * h)
        grad_analytic = analytic_grad[ix]
        rel_error = abs(grad_numerical - grad_analytic) / (abs(grad_numerical) + abs(grad_analytic) + 1e-8)
        sum_error += rel_error
        print('numerical: %f\tanalytic: %f\trelative error: %e' % (grad_numerical, grad_analytic, rel_error))
    return sum_error / num_checks
```

To use the gradient checker, we provide our analytical gradients, together with a function that computes the actual loss (rather than the gradients of the loss) and the parameters at which the gradient was computed:


```
#we'll use random parameters:
theta = torch.rand_like(theta_gd)*0.001  #rand_like Returns a tensor with the same size as input that is filled with random numbers from a uniform distribution on the interval [0, 1)
# and compute the analytic gradient (w.r.t the test data we loaded in this case)
grad = logistic_regression_loss_grad(theta, X_test, y_test)

# we need a function that computes the loss for a given theta (and implicitly the data)
def func(th):
    sigm = torch.sigmoid(X_test @ th)
    f = -(y_test.t() @ torch.log(sigm) + (1 - y_test.t()) @ torch.log(1 - sigm)); #Binary Cross Entropy
    return f

# and run the gradient checker
relerr = grad_check(func, theta, grad)
print("average error:", relerr)

assert relerr < 1e-5

```

    ix: (10, 0)
    oldval: 0.00013600983552179024
    numerical: 180.345802	analytic: 180.345802	relative error: 1.297988e-10
    ix: (44, 0)
    oldval: 0.0006849358283198793
    numerical: -297.998654	analytic: -297.998654	relative error: 8.906050e-11
    ix: (46, 0)
    oldval: 7.922916556820725e-05
    numerical: 155.510853	analytic: 155.510853	relative error: 2.962880e-11
    ix: (12, 0)
    oldval: 7.236484496028096e-05
    numerical: -41.644761	analytic: -41.644761	relative error: 1.030502e-09
    ix: (12, 0)
    oldval: 7.236484496028096e-05
    numerical: -41.644761	analytic: -41.644761	relative error: 1.030502e-09
    ix: (59, 0)
    oldval: 0.0002575770439390047
    numerical: 106.996047	analytic: 106.996047	relative error: 3.059648e-10
    ix: (21, 0)
    oldval: 0.000387605472563184
    numerical: 131.219843	analytic: 131.219843	relative error: 1.909091e-10
    ix: (27, 0)
    oldval: 0.000492927963875524
    numerical: -297.855481	analytic: -297.855480	relative error: 1.041454e-10
    ix: (30, 0)
    oldval: 0.00044110459900564136
    numerical: 162.684347	analytic: 162.684347	relative error: 1.729017e-11
    ix: (61, 0)
    oldval: 0.0001568127053011459
    numerical: -57.242986	analytic: -57.242986	relative error: 3.349651e-10
    average error: tensor([[3.2628e-10]])


Running the above, you should have a very small average error, and the relative error for each trial should also be a very small value.

f08362d1684e73075801bb4464791233cc2b6f34
ipynb
Jupyter Notebook
1_2_logistic_regression_submit_29299675.ipynb
mjjackey/DL_Lab_Soton
5df0dc3124e6fae6c27bfb99d70c457dd77935c5
Apache-2.0
2021-11-09T09:49:16.000Z
2021-11-09T09:49:16.000Z
1_2_logistic_regression_submit_29299675.ipynb
mjjackey/DL_Lab_Soton
5df0dc3124e6fae6c27bfb99d70c457dd77935c5
Apache-2.0
1_2_logistic_regression_submit_29299675.ipynb
mjjackey/DL_Lab_Soton
5df0dc3124e6fae6c27bfb99d70c457dd77935c5
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES