**Project Due:** March 1, 2021 at 11:59pm

**Late Policy:** Up to two slip days can be used for the final submission.

Please submit all required documents to CMS.

This is a group project. You can either work alone, or work __ ONLY__ with the other members of your group, which may contain AT MOST three people in total. Failure to adhere to this rule (by e.g. copying code) may result in an Academic Integrity Violation.

Overview: In this project, you will be implementing and evaluating gradient descent for empirical risk minimization on a multinomial logistic regression task on the MNIST dataset. We'll start by reviewing multinomial logistic regression as a method. Consider a multi-class classification task with $c$ classes in which the predictions $\hat y$ are distributions over the $c$ classes. That is, \[ \hat y \in \R^c, \hspace{2em} \hat y_j \ge 0 \text{ for all } j \in \{1, \ldots, c\} \hspace{2em} \text{and} \hspace{2em} \sum_{j=1}^c \hat y = 1. \] The example labels $y$ are one-hot vectors identifying a single class. That is, $y \in \R^c$, each $y_j \in \{0, 1\}$, and there exists exactly one $j \in \{1, \ldots, c\}$ such that $y_j = 1$. The examples themselves are vectors $x \in \R^d$. Multinomial logistic regression uses a linear model hypothesis with a softmax function. Specifically, given a parameter matrix $W \in \R^{c \times d}$, the hypothesis is \[ h_W(x) = \operatorname{softmax}(Wx) \] where the softmax function maps from $\R^c$ to $\R^c$ and is defined by \[ \left( \operatorname{softmax}(u) \right)_j = \frac{\exp(u_j)}{\sum_{k=1}^c \exp(u_k)}. \] We want to perform empirical risk minimization using the cross entropy loss function \[ L(\hat y, y) = -\sum_{j=1}^c y_j \cdot \log(\hat y_j), \] where here $\log$ refers to the natural logarithm, and using $\ell_2$ regularization. Our regularized empirical risk for some dataset $\mathcal{D}$, using regularization parameter $\gamma > 0$ is defined as \[ R(h_W) = \frac{1}{|\mathcal{D}|} \sum_{(x,y) \in \mathcal{D}} L(h_W(x), y) + \frac{\gamma}{2} \| W \|_F^2, \] where $\| W \|_F$ denotes the Frobenius norm (also known as Euclidean norm) is just the square root of the sum of the squares of the entries of $W$. \[ \| W \|_F=\sqrt(\sum_{i=1}^m \sum_{j=1}^n |w_{i,j}|^2 ) \]

Part 1: The gradient of the loss function. A key part of gradient descent is computing the gradient of the loss. For this problem, the gradient of the (unregularized) loss for a single example with respect to the parameters $W$ is \[ \nabla_W L(h_W(x), y) = \left( \operatorname{softmax}(Wx) - y \right) x^T. \]

Derive an expression for the gradient of the total loss with regularization over a whole dataset with respect to the weights $W$. Also write out an expression for a single update step of gradient descent. Include these expressions (and the derivation) in your report.

Implement a function `multinomial_logreg_loss_i`

to compute the loss of a single example $(x,y)$ given parameters $W$.

Implement a function `multinomial_logreg_grad_i`

to compute the gradient of a single example $(x,y)$ given parameters $W$.

One quick and rough way to check if your gradient function is *actually* the derivative of your loss function is to use the loss function to compute a numerical approximation of the derivative, and compare it to the value the gradient function produces. Recall that, for a function $f: \mathbb{R}^d \rightarrow \mathbb{R}$ and vectors $x$ and $v$,
\[
v^T \nabla f(x) = \lim_{\eta \rightarrow 0} \frac{f(x + \eta v) - f(x)}{\eta},
\]
so for sufficiently small $\eta$,
\[
v^T \nabla f(x) \approx \frac{f(x + \eta v) - f(x)}{\eta}.
\]
How might you use this formula to give a rough test of whether your gradient and loss functions are correct?
Write code to run the test you developed, using the function `test_gradient`

. `test_gradient`

is only meant to be self-contained i.e. don't integrate it with the rest of the steps. The functions will run through the autograder to determine if you are checking accurately, this could consist of print statements or any other way to show that the two values are approximately the same.
Include the results of your test and an explanation in your report.

Part 2: Gradient descent.

Implement the function `gradient_descent`

, the gradient descent algorithm. Your implementation should use the function `multinomial_logreg_total_grad`

which computes the total gradient of the loss function for the whole dataset, using your implementation of `multinomial_logreg_grad_i`

.

Run the gradient descent algorithm with L2 regularization parameter $\gamma = 0.0001$ and step size $\alpha = 1.0$ for 10 iterations, and record the value of the parameters every 10 iterations (i.e. only at the beginning and end). Time how long it took to run the algorithm. How long do you expect it would take to run 1000 iterations with this method? Report these results in your report.

Part 3: Computing more efficiently

Re-implement the functions `multinomial_logreg_total_loss`

and `multinomial_logreg_total_grad`

to use NumPy matrix arithmetic rather than Python for-loops to express computations of the loss and gradients of the whole dataset more efficiently. You should be able to implement this using no Python for-loops at all.

Re-run the gradient descent algorithm with the same parameters as in Part 2, for 10 epochs as before. Time how long it took to run the algorithm. How does this compare to your measurements from Part 2? Which approach is faster? Can you explain why? Report these results in your report.

Part 4: Evaluating gradient descent

Implement a function `multinomial_logreg_error`

to compute the error of the classifier given parameters $W$.

Run your gradient descent algorithm with L2 regularization parameter $\gamma = 0.0001$ and step size $\alpha = 1.0$ for 1000 iterations, and record the value of the parameters every 10 iterations (i.e. a total of 101 records).

Evaluate the training error, test error, training loss, and test loss of each recorded model exactly, using the entire MNIST training or test dataset. Plot the resulting error and loss with your results in part in four figures, one for Training error, one for Test error, one for Training loss, and one for Test loss. Here, the $x$-axis of your figures should report iterations.

Now write a function `estimate_multinomial_logreg_error`

to estimate the error of the model using subsampling.
Estimate the training error and test error of each recorded model using a random subsample of the dataset (of sizes 100 and 1000).
How do these results compare to the exact error?
What subsample size do you think you need to get a good sense of what the error is?
Explain.
Report these results in your report.

Measure the runtime of these error evaluations, comparing the exact (full dataset) error evaluating with your subsampled error evaluation method (at sizes 100 and 1000). Which one seems to be faster? Can you explain why? Report these results in your report.

What to submit:

- An implementation of the functions in main.py.
- A lab report containing:
- A one-paragraph summary of what you observed during this programming assignment.
- The content that you were instructed to include in the report in the directions above.

Setup:

- Run
`pip3 install -r requirements.txt`

to install the required python packages