# 3.6. Implementation of Softmax Regression from Scratch¶

Just as we implemented linear regression from scratch, we believe that
multiclass logistic (softmax) regression is similarly fundamental and
you ought to know the gory details of how to implement it from scratch.
As with linear regression, after doing things by hand we will breeze
through an implementation in Gluon for comparison. To begin, let’s
import our packages (only `autograd`

, `nd`

are needed here because
we will be doing the heavy lifting ourselves.)

```
import sys
sys.path.insert(0, '..')
%matplotlib inline
import d2l
from mxnet import autograd, nd
```

We will work with the Fashion-MNIST dataset just introduced, cuing up an iterator with batch size 256.

```
batch_size = 256
train_iter, test_iter = d2l.load_data_fashion_mnist(batch_size)
```

## 3.6.1. Initialize Model Parameters¶

Just as in linear regression, we represent each example as a vector. Since each example is a \(28 \times 28\) image, we can flatten each example, treating them as \(784\) dimensional vectors. In the future, we’ll talk about more sophisticated strategies for exploiting the spatial structure in images, but for now we treat each pixel location as just another feature.

Recall that in softmax regression, we have as many outputs as there are categories. Because our dataset has \(10\) categories, our network will have an output dimension of \(10\). Consequently, our weights will constitute a \(784 \times 10\) matrix and the biases will constitute a \(1 \times 10\) vector. As with linear regression, we will initialize our weights \(W\) with Gaussian noise and our biases to take the initial value \(0\).

```
num_inputs = 784
num_outputs = 10
W = nd.random.normal(scale=0.01, shape=(num_inputs, num_outputs))
b = nd.zeros(num_outputs)
```

Recall that we need to *attach gradients* to the model parameters. More
literally, we are allocating memory for future gradients to be stored
and notifiying MXNet that we want gradients to be calculated with
respect to these parameters in the first place.

```
W.attach_grad()
b.attach_grad()
```

## 3.6.2. The Softmax¶

Before implementing the softmax regression model, let’s briefly review
how operators such as `sum`

work along specific dimensions in an
NDArray. Given a matrix `X`

we can sum over all elements (default) or
only over elements in the same column (`axis=0`

) or the same row
(`axis=1`

). Note that if `X`

is an array with shape `(2, 3)`

and
we sum over the columns (`X.sum(axis=0`

), the result will be a (1D)
vector with shape `(3,)`

. If we want to keep the number of axes in the
original array (resulting in a 2D array with shape `(1,3)`

), rather
than collapsing out the dimension that we summed over we can specify
`keepdims=True`

when invoking `sum`

.

```
X = nd.array([[1, 2, 3], [4, 5, 6]])
X.sum(axis=0, keepdims=True), X.sum(axis=1, keepdims=True)
```

```
(
[[5. 7. 9.]]
<NDArray 1x3 @cpu(0)>,
[[ 6.]
[15.]]
<NDArray 2x1 @cpu(0)>)
```

We are now ready to implement the softmax function. Recall that softmax
consists of two steps: First, we exponentiate each term (using `exp`

).
Then, we sum over each row (we have one row per example in the batch) to
get the normalization constants for each example. Finally, we divide
each row by its normalization constant, ensuring that the result sums to
\(1\). Before looking at the code, let’s recall what this looks
expressed as an equation:

The denominator, or normalization constant, is also sometimes called the partition function (and its logarithm the log-partition function). The origins of that name are in statistical physics where a related equation models the distribution over an ensemble of particles).

```
def softmax(X):
X_exp = X.exp()
partition = X_exp.sum(axis=1, keepdims=True)
return X_exp / partition # The broadcast mechanism is applied here
```

As you can see, for any random input, we turn each element into a non-negative number. Moreover, each row sums up to 1, as is required for a probability. Note that while this looks correct mathematically, we were a bit sloppy in our implementation because failed to take precautions against numerical overflow or underflow due to large (or very small) elements of the matrix, as we did in Section 2.5.

```
X = nd.random.normal(shape=(2, 5))
X_prob = softmax(X)
X_prob, X_prob.sum(axis=1)
```

```
(
[[0.21324193 0.33961776 0.1239742 0.27106097 0.05210521]
[0.11462264 0.3461234 0.19401033 0.29583326 0.04941036]]
<NDArray 2x5 @cpu(0)>,
[1.0000001 1. ]
<NDArray 2 @cpu(0)>)
```

## 3.6.3. The Model¶

Now that we have defined the softmax operation, we can implement the
softmax regression model. The below code defines the forward pass
through the network. Note that we flatten each original image in the
batch into a vector with length `num_inputs`

with the `reshape`

function before passing the data through our model.

```
def net(X):
return softmax(nd.dot(X.reshape((-1, num_inputs)), W) + b)
```

## 3.6.4. The Loss Function¶

Next, we need to implement the cross entropy loss function, introduced in Section 3.4. This may be the most common loss function in all of deep learning because, at the moment, classification problems far outnumber regression problems.

Recall that cross entropy takes the negative log likelihood of the
predicted probability assigned to the true label \(-\log p(y|x)\).
Rather than iterating over the predictions with a Python `for`

loop
(which tends to be inefficient), we can use the `pick`

function which
allows us to select the appropriate terms from the matrix of softmax
entries easily. Below, we illustrate the `pick`

function on a toy
example, with 3 categories and 2 examples.

```
y_hat = nd.array([[0.1, 0.3, 0.6], [0.3, 0.2, 0.5]])
y = nd.array([0, 2], dtype='int32')
nd.pick(y_hat, y)
```

```
[0.1 0.5]
<NDArray 2 @cpu(0)>
```

Now we can implement the cross-entropy loss function efficiently with just one line of code.

```
def cross_entropy(y_hat, y):
return -nd.pick(y_hat, y).log()
```

## 3.6.5. Classification Accuracy¶

Given the predicted probability distribution `y_hat`

, we typically
choose the class with highest predicted probability whenever we must
output a *hard* prediction. Indeed, many applications require that we
make a choice. Gmail must catetegorize an email into Primary, Social,
Updates, or Forums. It might estimate probabilities internally, but at
the end of the day it has to choose one among the categories.

When predictions are consistent with the actual category `y`

, they are
coorect. The classification accuracy is the fraction of all predictions
that are correct. Although we cannot optimize accuracy directly (it is
not differentiable), it’s often the performance metric that we care most
about, and we will nearly always report it when training classifiers.

To compute accuracy we do the following: First, we execute
`y_hat.argmax(axis=1)`

to gather the predicted classes (given by the
indices for the largest entires each row). The result has the same shape
as the variable `y`

. Now we just need to check how frequently the two
match. Since the equality operator `==`

is datatype-sensitive (e.g. an
`int`

and a `float32`

are never equal), we also need to convert both
to the same type (we pick `float32`

). The result is an NDArray
containing entries of 0 (false) and 1 (true). Taking the mean yields the
desired result.

```
def accuracy(y_hat, y):
return (y_hat.argmax(axis=1) == y.astype('float32')).mean().asscalar()
```

We will continue to use the variables `y_hat`

and `y`

defined in the
`pick`

function, as the predicted probability distribution and label,
respectively. We can see that the first example’s prediction category is
2 (the largest element of the row is 0.6 with an index of 2), which is
inconsistent with the actual label, 0. The second example’s prediction
category is 2 (the largest element of the row is 0.5 with an index of
2), which is consistent with the actual label, 2. Therefore, the
classification accuracy rate for these two examples is 0.5.

```
accuracy(y_hat, y)
```

```
0.5
```

Similarly, we can evaluate the accuracy for model `net`

on the data
set (accessed via `data_iter`

).

```
# The function will be gradually improved: the complete implementation will be
# discussed in the "Image Augmentation" section
def evaluate_accuracy(data_iter, net):
acc_sum, n = 0.0, 0
for X, y in data_iter:
y = y.astype('float32')
acc_sum += (net(X).argmax(axis=1) == y).sum().asscalar()
n += y.size
return acc_sum / n
```

Because we initialized the `net`

model with random weights, the
accuracy of this model should be close to random guessing, i.e. 0.1 for
10 classes.

```
evaluate_accuracy(test_iter, net)
```

```
0.0925
```

## 3.6.6. Model Training¶

The training loop for softmax regression should look strikingly familiar
if you read through our implementation of linear regression earlier in
this chapter. Again, we use the mini-batch stochastic gradient descent
to optimize the loss function of the model. Note that the number of
epochs (`num_epochs`

), and learning rate (`lr`

) are both adjustable
hyper-parameters. By changing their values, we may be able to increase
the classification accuracy of the model. In practice we’ll want to
split our data three ways into training, validation, and test data,
using the validation data to choose the best values of our
hyperparameters.

```
num_epochs, lr = 5, 0.1
# This function has been saved in the d2l package for future use
def train_ch3(net, train_iter, test_iter, loss, num_epochs, batch_size,
params=None, lr=None, trainer=None):
for epoch in range(num_epochs):
train_l_sum, train_acc_sum, n = 0.0, 0.0, 0
for X, y in train_iter:
with autograd.record():
y_hat = net(X)
l = loss(y_hat, y).sum()
l.backward()
if trainer is None:
d2l.sgd(params, lr, batch_size)
else:
# This will be illustrated in the next section
trainer.step(batch_size)
y = y.astype('float32')
train_l_sum += l.asscalar()
train_acc_sum += (y_hat.argmax(axis=1) == y).sum().asscalar()
n += y.size
test_acc = evaluate_accuracy(test_iter, net)
print('epoch %d, loss %.4f, train acc %.3f, test acc %.3f'
% (epoch + 1, train_l_sum / n, train_acc_sum / n, test_acc))
train_ch3(net, train_iter, test_iter, cross_entropy, num_epochs,
batch_size, [W, b], lr)
```

```
epoch 1, loss 0.7900, train acc 0.745, test acc 0.803
epoch 2, loss 0.5736, train acc 0.812, test acc 0.825
epoch 3, loss 0.5291, train acc 0.823, test acc 0.833
epoch 4, loss 0.5058, train acc 0.830, test acc 0.838
epoch 5, loss 0.4901, train acc 0.834, test acc 0.837
```

## 3.6.7. Prediction¶

Now that training is complete, our model is ready to classify some images. Given a series of images, we will compare their actual labels (first line of text output) and the model predictions (second line of text output).

```
for X, y in test_iter:
break
true_labels = d2l.get_fashion_mnist_labels(y.asnumpy())
pred_labels = d2l.get_fashion_mnist_labels(net(X).argmax(axis=1).asnumpy())
titles = [truelabel + '\n' + predlabel
for truelabel, predlabel in zip(true_labels, pred_labels)]
d2l.show_fashion_mnist(X[0:9], titles[0:9])
```

## 3.6.8. Summary¶

With softmax regression, we can train models for multi-category classification. The training loop is very similar to that in linear regression: retrieve and read data, define models and loss functions, then train models using optimization algorithms. As you’ll soon find out, most common deep learning models have similar training procedures.

## 3.6.9. Exercises¶

In this section, we directly implemented the softmax function based on the mathematical definition of the softmax operation. What problems might this cause (hint - try to calculate the size of \(\exp(50)\))?

The function

`cross_entropy`

in this section is implemented according to the definition of the cross-entropy loss function. What could be the problem with this implementation (hint - consider the domain of the logarithm)?What solutions you can think of to fix the two problems above?

Is it always a good idea to return the most likely label. E.g. would you do this for medical diagnosis?

Assume that we want to use softmax regression to predict the next word based on some features. What are some problems that might arise from a large vocabulary?