# Notebook 7: Logistic Regression and SoftMax for MNIST

## Learning Goal

The goal of this notebook is to familiarize the reader with SoftMax regression (a generalization of logistic regression to more than two categories), categorical predictions, and the MNIST handwritten dataset. The reader will understand how to use the Scikit Logistic regression package and visualize learned weights.

## Overview
### The MNIST dataset:
The MNIST classification problem is one of the classical ML problems for learning classification on high-dimensional data with a fairly sizable number of examples (60000). Yann LeCun and collaborators collected and processed $70000$ handwritten digits (60000 are used for training and 10000 for testing) to produce what became known as one of the most widely used datasets in ML: the [MNIST](http://yann.lecun.com/exdb/mnist/) dataset. Each handwritten digit comes in a grayscale square image in the shape of a $28\times 28$ pixel grid. Every pixel takes a value in the range $[0,255]$, representing $256$ nuances of the gray color. The problem of image classification finds applications in a wide range of fields and is important for numerous industry applications of ML. 
### SoftMax regression:
We will use SoftMax regression, which can be thought of as a statistical model which assigns a probability that a given input image corresponds to any of the 10 handwritten digits. The model is a generalization of the logistic regression and reads:
\begin{align}
p(y=i|\boldsymbol{x};W) = \frac{e^{\boldsymbol{w}_i^T \boldsymbol{x}}}{\sum_{j=0}^9 e^{\boldsymbol{w}_j^T}},
\end{align}
Where $p(y=i|\boldsymbol{x};W)$ is the probability that input $\boldsymbol{x}$ is the $i$-th digit, $i\in[0,9]$.
The model also has 10 weight vectors $\boldsymbol{w}_i$ which we will train below. Finally, one can use this information for prediction by taking the value of $y$ for which this probability is maximized:
\begin{align}
y_{pred}=\arg\max_i p(y=i|\boldsymbol{x})
\end{align}

## Numerical Experiments

The reader is invited to check out the code below to build up their intuition about SoftMax regression. The following notebook is a slight modification of [this Scikit tutorial](http://scikit-learn.org/dev/auto_examples/linear_model/plot_sparse_logistic_regression_mnist.html) by Arthur Mensch on studying the MNIST problem using Logistic Regression.


```python
import time
import numpy as np
import random

from sklearn.datasets import fetch_openml # MNIST data
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.utils import check_random_state

print(__doc__)
```

    Automatically created module for IPython interactive environment



```python
# Turn down for faster convergence
t0 = time.time()
train_size = 50000
test_size = 10000

### load MNIST data from https://www.openml.org/d/554
X, y = fetch_openml('mnist_784', version=1, return_X_y=True)

# shuffle data
random_state = check_random_state(0)
permutation = random_state.permutation(X.shape[0])
X = X[permutation]
y = y[permutation]
X = X.reshape((X.shape[0], -1))

# pick training and test data sets 
X_train, X_test, y_train, y_test = train_test_split(X,y,train_size=train_size,test_size=test_size)

X_test_unscaled=  X_test.copy()

# scale data to have zero mean and unit variance [required by regressor]
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

# apply logistic regressor with 'sag' solver, C is the inverse regularization strength
clf = LogisticRegression(C=1e5,
                         multi_class='multinomial',
                         penalty='l2', solver='sag', tol=0.1)
# fit data
clf.fit(X_train, y_train)
# percentage of nonzero weights
sparsity = np.mean(clf.coef_ == 0) * 100
# compute accuracy
score = clf.score(X_test, y_test)

#display run time
run_time = time.time() - t0
print('Example run in %.3f s' % run_time)

print("Sparsity with L2 penalty: %.2f%%" % sparsity)
print("Test score with L2 penalty: %.4f" % score)
```

    Example run in 38.831 s
    Sparsity with L2 penalty: 9.18%
    Test score with L2 penalty: 0.8958



```python
%matplotlib inline
import matplotlib.pyplot as plt

# plot weights vs the pixel position
coef = clf.coef_.copy()
plt.figure(figsize=(10, 5))
scale = np.abs(coef).max()
for i in range(10):
    l2_plot = plt.subplot(2, 5, i + 1)
    l2_plot.imshow(coef[i].reshape(28, 28), interpolation='nearest',
                   cmap=plt.cm.Greys, vmin=-scale, vmax=scale)
    l2_plot.set_xticks(())
    l2_plot.set_yticks(())
    l2_plot.set_xlabel('Class %i' % i)
plt.suptitle('classification weights vector $w_j$ for digit class $j$')

plt.show()
```

### Exercises: ###  
<ul>
<li> The snippet above invoked the "sag" solver which only deals with $L2$ regularization. Try another solver, e.g. 'liblinear', in scikit that supports $L1$ regularization. You can find more details on <a href="http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html"> http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html.</a>
<li> Visualize the misclassified samples. Analyze on how logistic regression performs on MNIST, and note this since we will study the MNIST problem later on using deep neural nets.

<li> Try to vary the size of the training set and see how that affects the accuracy.
<li>  [Theory] In this notebook, we considered a $K=10$ multi-class logistic regression (MLR) also known as the SoftMax regression. Recall that the our model reads:
$$
p(y=i|\boldsymbol{x};W) = \frac{e^{\boldsymbol{w}_i^T \boldsymbol{x}}}{\sum_{j=0}^9 e^{\boldsymbol{w}_j^T \boldsymbol{x}}}
$$
Use this to construct the likelihood function $P(y|\boldsymbol{x})$.

<li> [Theory] Now consider we're given a set of i.i.d. data $\mathcal{D}=\{(\boldsymbol{x}_1,y_1),\cdots, (\boldsymbol{x}_n,y_n)\}$. As we argued in the main text, the cross-entropy error used in logistic regression is motivated by maximizing the likelihood (or equivalently, minimizing the negative log-likelihood), find the corresponding error function $E(\boldsymbol{w}_0,\cdots,\boldsymbol{w}_9)$ in MLR using the hypothesis above. 

<li> [Theory] Using the error function derived above, calculate its gradient $\frac{\partial E}{\partial \boldsymbol{w}_k}$. Can you find an elegant way of writing the gradient?
</ul>



```python
yhat=clf.predict(X_test)
fails = X_test_unscaled[yhat!=y_test]
fails_ys = y_test[yhat!=y_test]
fails_yhats = yhat[yhat!=y_test]
```


```python
numimages = 10

Is = [random.randint(0,fails.shape[0]) for j in range(numimages)]
for i in Is:
    plt.figure()
    plt.imshow(fails[i].reshape(28, 28), interpolation='nearest',
                   cmap='gray_r',vmin=0, vmax=255)
    plt.show()
    print("Predicted: "+fails_yhats[i]+", Actual: "+fails_ys[i])
#     vmin=np.abs(fails[i]).min(), vmax=np.abs(fails[i]).max()
```


```python
fails.shape
```




    (1042, 784)




```python
# # Turn down for faster convergence
# t0 = time.time()
# train_size = 50000
# test_size = 10000

# ### load MNIST data from https://www.openml.org/d/554
# X, y = fetch_openml('mnist_784', version=1, return_X_y=True)

# # shuffle data
# random_state = check_random_state(0)
# permutation = random_state.permutation(X.shape[0])
# X = X[permutation]
# y = y[permutation]
# X = X.reshape((X.shape[0], -1))

# # pick training and test data sets 
# X_train, X_test, y_train, y_test = train_test_split(X,y,train_size=train_size,test_size=test_size)

# X_test_unscaled=  X_test.copy()

# # scale data to have zero mean and unit variance [required by regressor]
# scaler = StandardScaler()
# X_train = scaler.fit_transform(X_train)
# X_test = scaler.transform(X_test)

# # apply logistic regressor with 'sag' solver, C is the inverse regularization strength
# clf = LogisticRegression(C=1e5,
# #                          multi_class='multinomial',
#                          penalty='l1', solver='liblinear', tol=0.1)
# # fit data
# clf.fit(X_train, y_train)
# # percentage of nonzero weights
# sparsity = np.mean(clf.coef_ == 0) * 100
# # compute accuracy
# score = clf.score(X_test, y_test)

# #display run time
# run_time = time.time() - t0
# print('Example run in %.3f s' % run_time)

# print("Sparsity with L2 penalty: %.2f%%" % sparsity)
# print("Test score with L2 penalty: %.4f" % score)
```

    Automatically created module for IPython interactive environment


    /Users/will/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/logistic.py:469: FutureWarning: Default multi_class will be changed to 'auto' in 0.22. Specify the multi_class option to silence this warning.
      "this warning.", FutureWarning)



```python

```

f290995fd97b1c3318b9b0bef9097d658ab6af95
ipynb
Jupyter Notebook
jupyter_notebooks/notebooks/NB7_CVII-logreg_mnist.ipynb
wberdanier/mlreview_notebooks
5434239133767aec39970350ed00b41f2402d49d
MIT
jupyter_notebooks/notebooks/NB7_CVII-logreg_mnist.ipynb
wberdanier/mlreview_notebooks
5434239133767aec39970350ed00b41f2402d49d
MIT
jupyter_notebooks/notebooks/NB7_CVII-logreg_mnist.ipynb
wberdanier/mlreview_notebooks
5434239133767aec39970350ed00b41f2402d49d
MIT
Qwen/Qwen-72B
1. YES
2. YES