# Part 2: Logistic Regression and Gradient Checking


```python
# 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.__


```python
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
    raise NotImplementedError()
    return grad
```


```python
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:


```python
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. 


```python
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))
```

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


```python
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

        oldval = x[ix].item()
        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:


```python
#we'll use random parameters:
theta = torch.rand_like(theta_gd)*0.001
# 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));
    return f

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

assert relerr < 1e-5

```

81347abb0812222b273e81249e0a681bc16fdcb8
ipynb
Jupyter Notebook
Pytorch Practical Tasks/1_2_logistic_regression.ipynb
VladimirsHisamutdinovs/deep-learning-pytorch
252a2d9e496c444fad1cb77a9e200afca9590aef
Apache-2.0
Pytorch Practical Tasks/1_2_logistic_regression.ipynb
VladimirsHisamutdinovs/deep-learning-pytorch
252a2d9e496c444fad1cb77a9e200afca9590aef
Apache-2.0
Pytorch Practical Tasks/1_2_logistic_regression.ipynb
VladimirsHisamutdinovs/deep-learning-pytorch
252a2d9e496c444fad1cb77a9e200afca9590aef
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES