<jupyter_start><jupyter_text>Gradient descent to solve linear regressionWe'll implement stochastic gradient descent from scratch, and use it to solve an ordinary least-squares linear regression.Caveat: we don't actually use gradient descent, or the normal equation, to solve this problem. [Read why](https://towardsdatascience.com/why-gradient-descent-and-normal-equation-are-bad-for-linear-regression-928f8b32fa4f). The purpose of this notebook is to look at how gradient descent works.--- First, let's look at what **gradient descent** is trying to do.The basic idea is that we have some function we'd like to minimize. If you think of the function as a physical surface, we're going to try to walk (or ski!) down it to a minimum, hopefully the global minimum. The function we're trying to minimize is the **cost function**. This is whatever we want it to be, but in our example it's going to be the mean squared error... hence the name 'least squares solution'.If we think of the cost function as a function of the **parameters** (weights, or coefficients), then we can use the partial derivative of the function with respect to each parameter to adjust that weight... although actually it turns out we can avoid updating the parameters separately by treating them as a single vector.  Load some dataWe'll use our 'rocks' data:<jupyter_code>import pandas as pd

df = pd.read_csv('../data/rocks.csv')

df.head()<jupyter_output><empty_output><jupyter_text>Let's try to predict `Vs` from `Vp` and `Rho_n`:<jupyter_code>data = df[['Vp', 'Rho_n', 'Vs']].values<jupyter_output><empty_output><jupyter_text>Random dataIf you prefer, you can make a dataset from a multivariate normal distribution, so we can control the statistics of the features and target.To do this, run the following code:import numpy as npnum_samples = 100The desired mean values of the sample.mu = np.array([5.0, 5.0, 10.0])The desired covariance matrix.r = np.array([        [  2.00,  0.00, -2.00],        [  0.00,  2.00,  2.00],        [ -2.00,  2.00,  4.00]    ])Generate the random samples.data = np.random.multivariate_normal(mu, r, size=num_samples)  Make `X` and `y` We can visualize the data:<jupyter_code>import seaborn as sns
import pandas as pd

df = pd.DataFrame(data, columns=['x0', 'x1', 'y'])

sns.pairplot(df, )
X = data[:, 0:2]
y = data[:, 2]

print(X.shape, y.shape)
from sklearn.preprocessing import StandardScaler

X = StandardScaler().fit_transform(X)<jupyter_output><empty_output><jupyter_text>Solve with `sklearn`<jupyter_code>from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

model = LinearRegression().fit(X, y)

y_pred = model.predict(X)
print(f"MSE: {mean_squared_error(y, y_pred):.3f}\n")

(w0, w1), b = model.coef_, model.intercept_
print(f"w0 = {w0:.3f}\nw1 = {w1:.3f}\nb = {b:.3f}")
import numpy as np
import matplotlib.pyplot as plt

X_plt = np.arange(-3, 3, 0.1)
y_plt_sk = model.coef_[0]*X_plt + model.coef_[1]*0 + model.intercept_

plt.scatter(X[:, 0], y) 
plt.plot(X_plt, y_plt_sk, color='r')
plt.show()<jupyter_output><empty_output><jupyter_text>Closed form solutionThere's actually a linear algebraic way to solve this equation. In principle, this is much faster and easier to do than gradient descent. In practice, it might not be possible, e.g. on very large datasets, or with custom loss functions.One way to solve the problem is with ordinary least-squares and the so-called [**normal equation**](https://en.wikipedia.org/wiki/Linear_least_squaresDerivation_of_the_normal_equations):$$ \mathbf{y} = \mathbf{w}\mathbf{X} $$We want the weights $\mathbf{w}$, so:$$ \mathbf{w} = \mathbf{X}^{-1}\mathbf{y} $$The catch is that $\mathbf{X}^{-1}$ may not be possible to compute, so we have to use [the pseudoinverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse) of `X`.This approach is easier if we add a column of 1's to `X`, and concatenate the weights `w` with the bias `b`, so that we have a single data matrix and a single weight vector:<jupyter_code>N = X.shape[0]
X_ = np.c_[X, np.ones((N, 1))]<jupyter_output><empty_output><jupyter_text>That [`np.c_`](https://numpy.org/doc/stable/reference/generated/numpy.c_.html) implements a concatenation trick. You can also use `np.hstack`.Now we can implement the 'naive' version of the normal equation:<jupyter_code>np.set_printoptions(suppress=True)
import numpy.linalg as la

X_inv = la.inv(X_.T @ X_) @ X_.T  # The 'pseudoinverse' of X.

X_inv @ y<jupyter_output><empty_output><jupyter_text>We solved it! Pretty cool. But it won't always work...  Gradient descentLet's solve the problem using gradient descent.In order to optimize the weights using gradient descent, we need to know the **gradient** of the **cost function**. The cost $C$ over $N$ samples is given by:$$ C = \frac{1}{N} \sum_{i=0}^{N} (\mathbf{y}_i - \hat{\mathbf{y}}_i)^2 $$where the estimate $\hat{\mathbf{y}}$ is given by the matrix product of the weights $\theta$ and the 'augmented' data $\bar{\mathbf{X}}$ (the data plus an extra column of ones):$$ \hat{\mathbf{y}} = \mathbf{w} \mathbf{X} + b = \theta \bar{\mathbf{X}} $$Using this 'combined' weight matrix is easier than dealing with the weights and bias term separately. Note that $\theta$ is a vector, although it's displayed with an italic character here. Chain ruleThere are 2 functions here: the expression for $C$ and the one for $\hat{\mathbf{y}}$. The chain rule tells us how to differentiate a function of a function:$$ \frac{\mathrm{d}}{\mathrm{d}x}[f(g(x))] = f^\prime(g(x))\ \times\ g^\prime(x) $$Therefore the gradient $g$ is given by:$$ g = \frac{\mathrm{d}C}{\mathrm{d}\theta} = \frac{\mathrm{d}C}{\mathrm{d}\hat{\mathbf{y}}}\ \times\ \frac{\mathrm{d}\hat{\mathbf{y}}}{\mathrm{d}\theta} $$So our problem now looks like this:$$ g = \frac{2}{N}\sum_{i=0}^{N} (\mathbf{y}_i - \hat{\mathbf{y}}_i)\ \times\ -\mathbf{X}_i $$so if the error $\mathbf{e} = (\mathbf{y} - \hat{\mathbf{y}})$ then:$$ g = -\frac{2}{N}\sum_{i=0}^{N} \mathbf{e}_i \mathbf{X}_i $$<jupyter_code>alpha = 0.001
epochs = 5_000

# Initialize.
cost_history = []
theta = np.zeros(3)

for i in range(epochs):
    # Make estimate with current weights.
    y_hat = theta @ X_.T
    
    # Compute error and cost.
    err = y - y_hat
    cost = np.sum(err**2) / N
    cost_history.append(cost)
    if i % (epochs // 10) == 0:
        print(f"{i:>8d} {cost:.3f}")
        
    # Compute gradients (one for each feature).
    grad = (-2/N) * np.sum(err[:, None] * X_, axis=0)
    
    # Adjust weights.
    theta -= alpha * grad

print(f"{'end':>8s} {cost:.3f}")
plt.plot(cost_history)
plt.yscale('log')<jupyter_output><empty_output><jupyter_text>Let's make sure it's similar to what we got before, which was<jupyter_code># With sklearn...
w0, w1, b
# With gradient descent...
theta<jupyter_output><empty_output><jupyter_text>Stochastic gradient descentHere, we're computing the cost over the entire dataset. In practice, this might not be possible, especially on high-dimensional datasets. We can reduce the computational burden by computing the cost over smaller 'batches'. This is known as 'stochastic gradient descent'. We'll also shuffle the data before each epoch (one epoch is one pass through the entire dataset). This helps to avoid getting stuck in a local minimum, because every time we take a step on the cost surface, we're on a slightly different version of the surface.<jupyter_code>alpha = 0.001
epochs = 1000
batch_size = 50

splits = np.arange(batch_size, N, batch_size)
idx = np.arange(N)

# Initialize.
cost_history = []
theta = np.random.random(3)

for i in range(epochs):
    
    # Shuffle the data.
    np.random.shuffle(idx)
    X_shuffle = X_[idx]
    y_shuffle = y[idx]
    X_batches = np.split(X_shuffle, splits)
    y_batches = np.split(y_shuffle, splits)

    for j, (X_batch, y_batch) in enumerate(zip(X_batches, y_batches)):

        # Make estimate with current weights.
        y_hat = theta @ X_batch.T

        # Compute error and (once per epoch) cost.
        err = y_batch - y_hat
        
        if j == 0:
            cost = np.sum(err**2) / batch_size
            cost_history.append(cost)

        # Compute gradients (one for each feature).
        grad = (-2/batch_size) * np.sum(err[:, None] * X_batch, axis=0)

        # Adjust weights.
        theta -= alpha * grad

    if i % (epochs // 10) == 0:
        print(f"{i:>8d} {cost:.3f}")

print(f"{'end':>8s} {cost:.3f}")
plt.plot(cost_history)
plt.yscale('log')<jupyter_output><empty_output><jupyter_text>Notice that this converges much faster than using the entire dataset.<jupyter_code>theta<jupyter_output><empty_output><jupyter_text>This compares well with what we got from `sklearn`:<jupyter_code>w0, w1, b<jupyter_output><empty_output><jupyter_text>And let's check the new parameters against the previous solution:<jupyter_code>w0_sgd, w1_sgd, b_sgd = theta

y_plt = w0_sgd*X_plt + w1_sgd*0 + b_sgd

plt.scatter(X[:, 0], y) 
plt.plot(X_plt, y_plt_sk, color='r')  # regression line
plt.plot(X_plt, y_plt, '--', color='c')  # regression line
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13357519788918207, 'fraction_numerical': 0.014951627088830254, 'mean_word_length': 4.332356389214537, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 53, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '68116', 'license': 'Apache-2.0', 'n_tokens_mistral': 3195, 'n_tokens_neox': 3026, 'n_words': 1011, 'path': 'notebooks/Demo__Gradient_descent.ipynb', 'repo_name': 'seismic-shift/geoml', 'size': '9096'}