<jupyter_start><jupyter_text>Linear Regression Let's implement a simple linear regression using torch<jupyter_code>import torch<jupyter_output><empty_output><jupyter_text>To start with, we'll create some random data and store it as torch tensors. For now you can think of a torch tensor as equivalent to a numpy array, but we will soon see that tensors have some additional functionality that we can exploit.<jupyter_code>n = 100
p = 5
x = torch.randn(n, p) # predictors (100 observations, 5 features)
y = torch.randn(n, 1) # outcomes (100 observations)<jupyter_output><empty_output><jupyter_text>The Model We're looking for parameters (coefficients) $\beta$ so that $$\begin{array}{rcl}y_i & \approx & \beta_0 + x_i\beta_{1:p}  \\&= &\beta_0 + x_{i1}\beta_1 +x_{i2}\beta_2 \dots x_{ip}\beta_p\end{array}$$ Although you may not have seen it represented this way before, we can also write this model as a picture: This will come in handy later when we get to more complex models, but for now you can think of the linear model in whatever way seems natural to you.  The Loss To quantify what we mean by a "good" approximation, we'll use the mean-squared-error loss. So for a given guess $\beta$ we'll give it the "grade" $$L(y,\hat y) = \frac{1}{n}\sum_i (y_i - \hat y_i)^2$$where $\hat y_i  = x_i\beta_{1:p} + \beta_0$ and $n$ is the number of observations (rows) in the data. We're looking for the $\beta$ that gives us the best (lowest) grade. This combination of model (linear) and loss (mean-squared-error) is called linear regression.  The Optimization Algorithm There are an infinite number of possible values for $\beta$ that we could try, so it would take us forever to try them all and see which gives the best loss. To get around this problem, we need some kind of optimization algorithm that is better than brute-force search. The algorithm we will use here is an extremely useful approach called gradient descent, which you should already be familiar with from our [previous exploration]().To start, we'll initialize the coefficients $\beta$ with random numbers. That's our first guess. We'll update these random numbers using gradient descent to iteratively find better values that make the loss smaller.<jupyter_code>β = torch.randn(p+1, 1) # 5 coefficients (one per feature) plus one intercept
β.requires_grad_() # tell torch that β is going to have to save the gradient of something with respect to itself at some point<jupyter_output><empty_output><jupyter_text>Autograd Previously when using gradient descent we would have to analytically derive expressions for the gradient of the loss relative to each model parameter, implement these as functions in code, and call upon them at each gradient descent iteration. With pytorch, however, we can simply compute the current value of the loss and pytorch will automatically calculate all the necessary derivatives for us. Let's have a look.<jupyter_code>ŷ = torch.matmul(x, β[1:]) + β[0] # ŷ = xβ + β0
L = torch.sum((y-ŷ)**2)/n # L = Σ(yᵢ-ŷᵢ)²/n
L
L.backward() # compute the gradients of the loss with respect to any tensors that went into the loss with requires_grad=true
β.grad<jupyter_output><empty_output><jupyter_text>what we see here is a vector containing all of the derivatives we want. The first element is $\frac{\delta L}{\delta \beta_0}$, the second element is $\frac{\delta L}{\delta \beta_1}$, and so on. Note that this object is part of $\beta$ and not $L$. Let's calculate the gradient manually and make sure it matches up. We have:$$\begin{array}{rcl}\hat y_i &=& x_i\beta + \beta_0 \\L(y,\hat y) &=& \frac{1}{n}\sum_i (y_i - \hat y_i)^2\end{array}$$ So the derivative for $\beta_j$ with $j\ne 0$ is$$\begin{array}{rcl}\frac{\partial L}{\partial \beta_j} &=& \frac{1}{n}\sum_i \frac{\partial L}{\partial y_i} \frac{\partial y_i}{\partial \beta_j}\\&=&\frac{1}{n}\sum_i-2(y_i-\hat y_i)x_{ij}\\&=&-\frac{2}{n}x_j^T(y-\hat y)\end{array}$$and for $\beta_0$ is $$\begin{array}{rcl}\frac{\partial L}{\partial \beta_0}&=&\frac{1}{n}\sum_i-2(y_i-\hat y_i) \\&=&-\frac{2}{n}1^T(y-\hat y)\end{array}$$which means we can calculate the whole gradient as $-\frac{2}{n}[1,x]^T(y-[1,x]\beta)$<jupyter_code>import numpy as np
def beta_grad(β):
    x_with_1s = torch.Tensor(np.hstack((np.ones((n,1)), x)))
    return -2*torch.matmul(x_with_1s.transpose(0,1), y-torch.matmul(x_with_1s,β))/n
beta_grad(β)<jupyter_output><empty_output><jupyter_text>which is exactly the same as the result in `β.grad`! How does pytorch do this? It turns out that every pytorch tensor records not only its own value, but also what functions were called to produce it and which tensors went into those functions. That's why we use torch tensors instead of numpy arrays. The `torchvis` package lets us see this:<jupyter_code>from torchviz import make_dot
make_dot(L, {'β':β})<jupyter_output><empty_output><jupyter_text>To calculate the gradient of $L$ with respect to $\beta$, pytorch takes the gradients of each of these functions in turn (which are simple and hardcoded into pytorch), evaluates them at their current value using the input tensors, and multiplies them together to arrive at the answer you would get via the chain rule. You can learn more about this process [here](https://towardsdatascience.com/pytorch-autograd-understanding-the-heart-of-pytorchs-magic-2686cd94ec95) and elsewhere.  Doing Gradient Descent To update the weights, we need to subtract the gradient (times a small learning rate) from the current value of the weights. We do this from inside a `no_grad():` "context" so that $\beta$ doesn't store the history of the update (try the update without the `with torch.no_grad():` and see what happens). We also clear out `β.grad`.<jupyter_code>with torch.no_grad():
    β -= 10e-5 * β.grad # β = β - 10e-5 * β.grad
    β.grad.zero_()
β<jupyter_output><empty_output><jupyter_text>As you can see, our new value of $\beta$ is slightly different than what we started with because we've taken a single gradient step. **ZEROING GRADIENTS**If we don't zero the gradient, the next time we calculate the gradient of something with respect to $\beta$, the new gradient will be added to whatever was stored there instead of overwriting it. That's just the way torch was made to work. ---**EXERCISE**Investigate for yourself (either by testing code or googling) why the parameter update and gradient zerioing should be performed within a `torch.no_grad()` context. ---  Looping Let's repeat what we have so far but now add a little loop to train our model for 500 gradient descent iterations instead of going slowly though a single iteration:<jupyter_code>x = torch.randn(n, p) # predictors (10 observations, 5 features)
y = torch.randn(n, 1) # outcomes (10 observations)

loss_record = [] # to kep track of the loss over the iterations
β = torch.randn(p+1,1) # 5 coefficients (one per feature) plus one intercept
β.requires_grad_() # tell torch that β is going to have to save the gradient of something with respect to itself at some point
for e in range(500):
    ŷ = torch.matmul(x, β[1:]) + β[0] # ŷ = xβ + β0 (calculate predictions)
    L = torch.sum((y-ŷ)**2)/n # L = Σ(yᵢ-ŷᵢ)²/n (use predictions to calculate loss)
    
    L.backward() # compute gradients (in this case δL/δβ, δL/δW)
    loss_record.append(L.item())
    
    with torch.no_grad(): # take the gradient descent step 
        β -= 10e-3 * β.grad
        β.grad.zero_()<jupyter_output><empty_output><jupyter_text>And we can see how the loss changes over the iterations:<jupyter_code>import altair as alt
import pandas as pd

loss_df = pd.DataFrame({'i':range(500), 'loss':loss_record})
alt.Chart(loss_df, height=100).mark_line().encode(x='i', y='loss')<jupyter_output><empty_output><jupyter_text>We see training loss goes down. And here's the value of $\beta$ after 500 iterations:<jupyter_code>β<jupyter_output><empty_output><jupyter_text>Remember our original data is totally random and there is no relationship between the predictors and outcomes. So the "right" answer for what $\beta$ should be in this case is $\beta=0$. As we see above, all the values are near 0, so our algorithm appears to be converging to the right answer. To predict for a new observation, all we have to do is multiply by $\beta$ and add the intercept:<jupyter_code>new_x = torch.randn(10, 5) # 10 new observations
torch.matmul(new_x, β[1:]) + β[0] # predictions<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup