# Part 3: Softmax Regression


```
# 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 second part of the lab we saw how to make a linear binary classifier using logisitic regression. In this  part of the lab we'll turn our attention to multi-class classification.

Softmax regression (or multinomial logistic regression) is a generalisation of logistic regression to the case where we want to handle multiple classes. In logistic regression we assumed that the labels were binary: $y_i\in \{0,1\}$. We used such a classifier to distinguish between two kinds of hand-written digits. Softmax regression allows us to handle $y_i \in \{1,\dots,K\}$ where $K$ is the number of classes.

Recall that in logistic regression, we had a training set $\{(\mathbf{x}_1,y_1),\dots,(\mathbf{x}_m,y_m)\}$ of $m$ labeled examples, where the input features are $\mathbf{x}_i \in \mathbb{R}^n$. In logistic regression, our hypothesis took the form:

\begin{align}
h_\theta(\mathbf{x}) &= \frac{1}{1 + \exp(-\mathbf{x}^\top\theta)} \equiv \sigma(\mathbf{x}^\top\theta)
\end{align}

and the model parameters $\theta$ were trained to minimise the cost function

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

In the softmax regression setting, we are interested in multi-class classification, and so the label $y$
 can take on $K$ different values, rather than only two. Thus, in our training set $\{(\mathbf{x}_1,y_1),\dots,(\mathbf{x}_m,y_m)\}$, we now have that $y_i \in \{1,\dots,K\}$.

Given a test input $\mathbf{x}$, we want our hypothesis to estimate the probability that $P(y=k|\mathbf{x})$ for each value of $k=1,\dots,K$. That is to say, we want to estimate the probability of the class label taking on each of the $K$ different possible values. Thus, our hypothesis will output a $K$-dimensional vector (whose elements sum to 1) giving us our $K$ estimated probabilities. Concretely, our hypothesis $h_\theta(\mathbf{x})$ takes the form:

\begin{align}
h_\theta(\mathbf{x}) =
\begin{bmatrix}
P(y = 1 | \mathbf{x}; \theta) \\
P(y = 2 | \mathbf{x}; \theta) \\
\vdots \\
P(y = K | \mathbf{x}; \theta)
\end{bmatrix}
=
\frac{1}{ \sum_{j=1}^{K}{\exp(\theta^{(j)\top} \mathbf{x}) }}
\begin{bmatrix}
\exp(\theta^{(1)\top} \mathbf{x} ) \\
\exp(\theta^{(2)\top} \mathbf{x} ) \\
\vdots \\
\exp(\theta^{(K)\top} \mathbf{x} ) \\
\end{bmatrix}
\end{align}

Here $\theta^{(1)},\theta^{(2)},\dots,\theta^{(K)} \in \mathbb{R}^n$ are the parameters of our model. Notice that the term $\frac{1}{\sum_{j=1}^K exp(\theta^{(j)\top} \mathbf{x})}$ normalizes the distribution, so that it sums to one.

For convenience, we will also write $\theta$ to denote all the parameters of our model. When you implement softmax regression, it is usually convenient to represent $\theta$ as a $n$-by-$K$ matrix obtained by concatenating $\theta_{(1)},\theta^{(2)},\dots,\theta^{(K)}$ into columns, so that

\begin{align}
\theta = \left[\begin{array}{cccc}| & | & | & | \\
\theta^{(1)} & \theta^{(2)} & \cdots & \theta^{(K)} \\
| & | & | & |
\end{array}\right].
\end{align}


## Cost Function

We now describe the cost function that we’ll use for softmax regression. In the equation below, $1\{\cdot\}$
 is an "indicator function", such that $1\{\mathrm{a true statement}\}=1$, and $1\{\mathrm{a false statement}\}=0$. For example, $1\{2+2=4\}$ evaluates to $1$; whereas $1\{1+1=5\}$ evaluates to $0$. Our cost function will be:

\begin{align}
J(\theta) = - \left[ \sum_{i=1}^{m} \sum_{k=1}^{K}  1\left\{y_{i} = k\right\} \log \frac{\exp(\theta^{(k)\top} \mathbf{x}_i)}{\sum_{j=1}^K \exp(\theta^{(j)\top} \mathbf{x}_i)}\right]
\end{align}
 
Notice that this generalises the logistic regression cost function, which could also have been written:

\begin{align}
J(\theta) &= - \left[ \sum_{i=1}^m   (1-y^{(i)}) \log (1-h_\theta(\mathbf{x}_i)) + y^{(i)} \log h_\theta(\mathbf{x}_i) \right] \\
&= - \left[ \sum_{i=1}^{m} \sum_{k=0}^{1} 1\left\{y^{(i)} = k\right\} \log P(y^{(i)} = k | \mathbf{x}_i ; \theta) \right]
\end{align}


The softmax cost function is similar, except that we now sum over the $K$ different possible values of the class label. Note also that in softmax regression, we have that

\begin{equation}
P(y_i = k | \mathbf{x}_i ; \theta) = \frac{\exp(\theta^{(k)\top} \mathbf{x}_i)}{\sum_{j=1}^K \exp(\theta^{(j)\top} \mathbf{x}_i) }
\end{equation}

We cannot solve for the minimum of $J(\theta)$ analytically, and thus we'll resort to using gradient descent as before. Taking derivatives, one can show that the gradient is:

\begin{align}
\nabla_{\theta^{(k)}} J(\theta) = - \sum_{i=1}^{m}{ \left[ \mathbf{x}_i \left( 1\{ y_i = k\}  - P(y_i = k | \mathbf{x}_i; \theta) \right) \right]  }
\end{align}

Armed with this formula for the derivative, one can then use it directly with a gradient descent solver (or any other 1st-order gradient based optimiser).

__Use the code box below to complete the implementation of the functions that return the gradients of the softmax loss function, $\nabla_{\theta^{(k)}} J(\theta) \,\, \forall k$ and the loss function itself, $J(\theta)$:__


```
import torch
import torch.nn.functional as f

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

'''
Theta: f*K(20*10)
X: num_items*f(100*20)
y: num_items*1(100*1) value:0~9
out: f*K
'''
def softmax_regression_loss_grad(Theta, X, y):
    '''Implementation of the gradient of the softmax loss function.
    
    Theta is the matrix of parameters, with the parameters of the k-th class in the k-th column
    X contains the data vectors (one vector per row)
    y is a column vector of the targets
    '''
    # YOUR CODE HERE
    num_classes = Theta.size()[1]  #K
    features_dim = Theta.size()[0]  #f
    num_items = y.size()[0]
    y_label = torch.zeros(num_items,num_classes,dtype=torch.float64).scatter_(1,y,1)  #num_items * K
    prob = f.softmax(X @ Theta, dim=1) #num_items*K, dim=1 for column
    grad = - X.t() @ (y_label - prob)
    return grad

def exp_prob(Theta_k,Theta,x_i):
    sum=0
    for j in range(Theta.shape[1]):  #K
        sum+=torch.exp(x_i @ Theta[:,j])
    prob=torch.div(torch.exp(x_i @ Theta_k),sum) 
    return prob  
  
def softmax_regression_loss(Theta, X, y):
    '''Implementation of the softmax loss function.
        
    Theta is the matrix of parameters, with the parameters of the k-th class in the k-th column
    X contains the data vectors (one vector per row)
    y is a column vector of the targets
    '''
    # YOUR CODE HERE
    loss=0.0
    for i in range(X.shape[0]): #num_items
        x_i=X[i,:]  #the i-th row,1*f
        y_i=y[i]
        class_indicator=0
        for k in range(Theta.shape[1]): #K 
            Theta_k = Theta[:,k] #the k-th column,f*1
            class_inidicator = int(y_i.numpy()==k)  #numpy
            prob=exp_prob(Theta_k,Theta,x_i)
            loss += float(class_inidicator) * torch.log(prob)  #float,* seperate
    return -loss  #-

def softmax_regression_loss2(Theta, X, y):
    num_classes = Theta.size()[1]  #K
    features_dim = Theta.size()[0]  #f
    num_items = y.size()[0]
    y_label=torch.zeros(num_items,num_classes,dtype=torch.float64).scatter_(1,y,1)  #num_items * K
    prob=-f.log_softmax(X @ Theta, dim=1)  #num_items * K
    loss=torch.sum(y_label*prob)
    return loss
   
```

__Use the following code block to confirm that your implementation is correct using gradient checking. If there are problems with your gradient or loss, go back and fix them!:__


```
 # from torch.autograd import gradcheck
from random import randrange

def grad_check(f, x, analytic_grad, num_checks=10, h=1e-3):
    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)
        print("x[ix]:",x[ix])
        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

# Create some test data:
num_classes = 10   #K
features_dim = 20  #f
num_items = 100
Theta = torch.randn((features_dim, num_classes)) #features_dim * num_classes(f*K), filled with random numbers from a normal distribution
X = torch.randn((num_items,features_dim))  #num_items * features_dim(num_items*f)
y = torch.torch.randint(0, num_classes, (num_items, 1)) #0~num_classes-1, size:num_items*1

# compute the analytic gradient
grad = softmax_regression_loss_grad(Theta, X, y)
    
# run the gradient checker    
grad_check(lambda th: softmax_regression_loss(th, X, y), Theta, grad)
```

    ix: (7, 3)
    x[ix]: tensor(0.0487)
    numerical: -3.040584	analytic: -3.040584	relative error: 8.887084e-09
    ix: (13, 4)
    x[ix]: tensor(-1.0416)
    numerical: -3.542292	analytic: -3.542292	relative error: 1.051560e-08
    ix: (1, 7)
    x[ix]: tensor(-2.0798)
    numerical: -6.047417	analytic: -6.047417	relative error: 6.418247e-09
    ix: (3, 4)
    x[ix]: tensor(1.8285)
    numerical: 1.309945	analytic: 1.309945	relative error: 7.116945e-08
    ix: (2, 5)
    x[ix]: tensor(1.9467)
    numerical: 14.097420	analytic: 14.097420	relative error: 3.242193e-09
    ix: (13, 0)
    x[ix]: tensor(0.9674)
    numerical: 9.549632	analytic: 9.549632	relative error: 4.199948e-09
    ix: (1, 7)
    x[ix]: tensor(-2.0798)
    numerical: -6.047417	analytic: -6.047417	relative error: 6.418247e-09
    ix: (12, 9)
    x[ix]: tensor(0.0319)
    numerical: 1.614008	analytic: 1.614008	relative error: 1.278455e-08
    ix: (4, 9)
    x[ix]: tensor(0.4010)
    numerical: 0.039318	analytic: 0.039318	relative error: 3.090074e-06
    ix: (6, 5)
    x[ix]: tensor(0.4945)
    numerical: 4.357503	analytic: 4.357503	relative error: 6.102861e-08





    tensor(3.2747e-07)




```
Theta = torch.Tensor([[1, 0], [0, 1]])
X = torch.Tensor([[1, 0], [0, 1]])
y = torch.LongTensor([[0], [1]])
loss=softmax_regression_loss(Theta, X, y) 
print("loss:",loss)
assert torch.abs(softmax_regression_loss(Theta, X, y) - 0.6265) < 0.0001
grad = softmax_regression_loss_grad(Theta, X, y)
assert torch.torch.allclose(torch.abs(grad/0.2689), torch.ones_like(grad), atol=0.001)  
#This function checks if all self and other satisfy the condition:∣self−other∣≤atol+rtol×∣other∣, rtol=1e-08

```

    loss: tensor(0.6265)


## Training Softmax regression with gradient descent on real data

We'll now try gradient descent with our softmax regression using the digits dataset. As before, when we looked at logistic regression, we load the data and create test and training sets. Note that this time we'll use all the classes:


```
from sklearn.datasets import load_digits

X, y = (torch.Tensor(z) for z in load_digits(10, 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
X /= 255
y = y.reshape(-1, 1) # reshape y into a column vector
y = y.type(torch.LongTensor)

# 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], :]  #261*f
y_train = y[perm[0:260]]
X_test = X[perm[260:], :]
y_test = y[perm[260:]]
```

We now define a simple gradient descent loop to train the model:


```
alpha = 0.1
theta_gd = torch.rand((X_train.shape[1], 10)) #Returns a tensor filled with random numbers from a uniform distribution on the interval [0, 1)
#num_items*10

for e in range(0, 1000):
    gr = softmax_regression_loss_grad(theta_gd, X_train, y_train)
    theta_gd -= alpha * gr
    if e%100 == 0:
        print("Training Loss: ", softmax_regression_loss(theta_gd, X_train, y_train))

# Compute the accuracy of the test set
proba = torch.softmax(X_test @ theta_gd, 1)
print(float((proba.argmax(1)-y_test[:,0]==0).sum()) / float(proba.shape[0]))  #argmax(1) returns the indices of the maximum values of prob
print()

```

    Training Loss:  tensor(592.3777)
    Training Loss:  tensor(270.9790)
    Training Loss:  tensor(174.5995)
    Training Loss:  tensor(132.5848)
    Training Loss:  tensor(109.0211)
    Training Loss:  tensor(93.7592)
    Training Loss:  tensor(82.9511)
    Training Loss:  tensor(74.8239)
    Training Loss:  tensor(68.4447)
    Training Loss:  tensor(63.2739)
    0.9232270657124269
    


Running the above, you should observe that the training loss decreases over time. The final accuracy on the test set is also printed and should be around 90% (it will depend on the particular training/test splits you generated as well as the initial parameters for the softmax).

# Overparameterisation in softmax regression

Softmax regression has an unusual property that it has a "redundant" set of parameters. To explain what this means, suppose we take each of our parameter vectors $\theta^{(j)}$, and subtract some fixed vector $\psi$ from it, so that every $\theta^{(j)}$ is now replaced with $\theta^{(j)}−\psi$ (for every $j=1,\dots,k$). Our hypothesis now estimates the class label probabilities as

\begin{align}
P(y^{(i)} = k | x^{(i)} ; \theta)
&= \frac{\exp((\theta^{(k)}-\psi)^\top x^{(i)})}{\sum_{j=1}^K \exp( (\theta^{(j)}-\psi)^\top x^{(i)})}  \\
&= \frac{\exp(\theta^{(k)\top} x^{(i)}) \exp(-\psi^\top x^{(i)})}{\sum_{j=1}^K \exp(\theta^{(j)\top} x^{(i)}) \exp(-\psi^\top x^{(i)})} \\
&= \frac{\exp(\theta^{(k)\top} x^{(i)})}{\sum_{j=1}^K \exp(\theta^{(j)\top} x^{(i)})}.
\end{align}

__In other words, subtracting $\psi$ from every $\theta^{(j)}$ does not affect our hypothesis’ predictions at all!__ This shows that softmax regression’s parameters are "redundant". More formally, we say that our softmax model is "overparameterised" meaning that for any hypothesis we might fit to the data, there are multiple parameter settings that give rise to exactly the same hypothesis function $h_\theta$ mapping from inputs $\mathbf{x}$ to the predictions.

Further, if the cost function $J(\theta)$ is minimized by some setting of the parameters $(\theta^{(1)},\theta^{(2)},\dots,\theta^{(k)})$, then it is also minimised by $\theta^{(1)}-\psi,\theta^{(2)}-\psi,\dots,\theta^{(k)}-\psi)$ for any value of $\psi$. Thus, the minimiser of $J(\theta)$ is not unique. 

(Interestingly, $J(\theta)$ is still convex, and thus gradient descent will not run into local optima problems. The Hessian is however singular/non-invertible, which causes a straightforward implementation of Newton's method (a second-order optimiser) to run into numerical problems.)

Notice also that by setting $\psi=\theta^{(K)}$, one can always replace $\theta^{(K)}$ with $\theta^{(K)}-\psi=\mathbf{0}$ (the vector of all $0$’s), without affecting the hypothesis. Thus, one could "eliminate" the vector of parameters $\theta^{(K)}$ (or any other $\theta^{(k)}$, for any single value of $k$), without harming the representational power of our hypothesis. Indeed, rather than optimising over the $K \cdot n$ parameters $(\theta^{(1)},\theta^{(2)},\dots,\theta^{(k)})$  (where $\theta^{(k)} \in \mathbb{R}^n$, one can instead set $\theta^{(K)}=\mathbf{0}$  and optimize only with respect to the $(K-1) \cdot n$ remaining parameters.

__Use the following block to implement the softmax gradients for the case where the final column of the parameters theta is fixed to be zero:__


```
import torch
import torch.nn.functional as f

'''
Theta: f*K(20*10)
X: num_items*f(100*20)
y: num_items*1(100*1) value:0~9
out: f*K
'''
def softmax_regression_loss_grad_0(Theta, X, y):
    '''Implementation of the gradient of the softmax loss function, with the parameters of the
    last class fixed to be zero.
    
    Theta is the matrix of parameters, with the parameters of the k-th class in the k-th column; 
            K-1 classes are included, and the parameters of the last class are implicitly zero.
    X contains the data vectors (one vector per row)
    y is a column vector of the targets
    '''
    
    # add the missing column of zeros:
    Theta = torch.cat((Theta, torch.zeros(Theta.shape[0],1)), 1)

    # YOUR CODE HERE
    num_classes = Theta.size()[1]  #(K+1)
    features_dim = Theta.size()[0]  #f
    num_items = y.size()[0]
    y_label = torch.zeros(num_items,num_classes,dtype=torch.float64).scatter_(1,y,1)  #num_items * (K+1)
    prob = f.softmax(X @ Theta, dim=1) #num_items*K, dim=1 for column
    grad = - X.t() @ (y_label - prob)  
  
    # remove the last column from the gradients
    grad = grad[0:grad.shape[0], 0:grad.shape[1]-1]
    return grad
```


```
Theta = torch.Tensor([[1, 0], [0, 0]])
X = torch.Tensor([[1, 0], [0, 1]])
y = torch.LongTensor([[0], [1]])

grad = softmax_regression_loss_grad(Theta, X, y)
grad0 = softmax_regression_loss_grad_0(Theta[:,0:grad.shape[1]-1], X, y)
assert torch.torch.allclose(grad[:,0:grad.shape[1]-1], grad0)

```

Finally, we can run gradient descent with our reduced paramter gradient function, and confirm that the results are similar to before:


```
alpha = 0.1
theta_gd = torch.rand((X_train.shape[1], 9))

for e in range(0, 1000):
    gr = softmax_regression_loss_grad_0(theta_gd, X_train, y_train)
    theta_gd -= alpha * gr

theta_gd = torch.cat((theta_gd, torch.zeros(theta_gd.shape[0], 1)), 1)
proba = torch.softmax(X_test @ theta_gd, 1)
print(float((proba.argmax(1)-y_test[:,0]==0).sum()) / float(proba.shape[0]))
print()
```

    0.9108653220559532
    



```

```

8c92d073da9901855b4194df2df12f93ef2a920a
ipynb
Jupyter Notebook
1_3_softmax_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_3_softmax_regression_submit_29299675.ipynb
mjjackey/DL_Lab_Soton
5df0dc3124e6fae6c27bfb99d70c457dd77935c5
Apache-2.0
1_3_softmax_regression_submit_29299675.ipynb
mjjackey/DL_Lab_Soton
5df0dc3124e6fae6c27bfb99d70c457dd77935c5
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES