# Imports


```python
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
```


```python
plt.style.use(["ggplot"])
```

# Create Data

<h5> Generate some data with:
\begin{equation} \theta_0= 4 \end{equation} 
\begin{equation} \theta_1= 3 \end{equation} 

Add some Gaussian noise to the data

Here we have $$y=\theta_0 + \theta_1 *X+ random\ noise$$


```python
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)
```

Let's plot our data to check the relation between X and Y

## plot generated 2d data points


```python
plt.plot(X, y, "b.")
plt.xlabel("$x$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
ax = plt.axis([0, 2, 0, 15])
```

#  Analytical way of Linear Regression


```python
# ?? np.c_
```


```python
X_b = np.c_[np.ones((100, 1)), X]
X_b
```




    array([[1.        , 0.61161301],
           [1.        , 1.72430602],
           [1.        , 1.85534309],
           [1.        , 1.76073836],
           [1.        , 0.0572145 ],
           [1.        , 1.84703953],
           [1.        , 0.23746452],
           [1.        , 1.66383103],
           [1.        , 0.2041836 ],
           [1.        , 1.51621319],
           [1.        , 1.40221905],
           [1.        , 1.68842929],
           [1.        , 0.23088307],
           [1.        , 1.96224888],
           [1.        , 1.38270468],
           [1.        , 0.45367799],
           [1.        , 1.51284161],
           [1.        , 1.57549713],
           [1.        , 0.6783873 ],
           [1.        , 1.89365556],
           [1.        , 0.21850822],
           [1.        , 1.85162873],
           [1.        , 0.76138162],
           [1.        , 1.21935501],
           [1.        , 1.52439915],
           [1.        , 0.80818048],
           [1.        , 0.70451694],
           [1.        , 0.88903564],
           [1.        , 0.32331054],
           [1.        , 1.64909732],
           [1.        , 1.7124657 ],
           [1.        , 0.91713771],
           [1.        , 0.57914975],
           [1.        , 0.29721128],
           [1.        , 1.10640219],
           [1.        , 1.90267708],
           [1.        , 0.52570261],
           [1.        , 0.76683326],
           [1.        , 1.60819434],
           [1.        , 1.78457635],
           [1.        , 1.17548247],
           [1.        , 0.50681977],
           [1.        , 0.04839011],
           [1.        , 1.71238211],
           [1.        , 1.17437814],
           [1.        , 1.35064057],
           [1.        , 1.46373702],
           [1.        , 0.42525494],
           [1.        , 0.82589216],
           [1.        , 1.09337039],
           [1.        , 1.53121379],
           [1.        , 1.74522245],
           [1.        , 0.33265673],
           [1.        , 0.70568204],
           [1.        , 0.55856574],
           [1.        , 1.92412633],
           [1.        , 0.38815032],
           [1.        , 1.33906861],
           [1.        , 0.82516803],
           [1.        , 1.5456871 ],
           [1.        , 0.71522878],
           [1.        , 0.08590715],
           [1.        , 0.58955981],
           [1.        , 0.61951299],
           [1.        , 0.97332899],
           [1.        , 0.44926547],
           [1.        , 1.94176703],
           [1.        , 0.18386009],
           [1.        , 0.15639567],
           [1.        , 1.76878058],
           [1.        , 0.4371604 ],
           [1.        , 0.90612205],
           [1.        , 0.20664235],
           [1.        , 1.58972011],
           [1.        , 1.62836902],
           [1.        , 1.75144256],
           [1.        , 0.40163699],
           [1.        , 1.67312172],
           [1.        , 1.81273512],
           [1.        , 0.44030981],
           [1.        , 1.60587802],
           [1.        , 1.79572153],
           [1.        , 0.16385152],
           [1.        , 0.00740951],
           [1.        , 1.09224804],
           [1.        , 1.45951573],
           [1.        , 1.27301069],
           [1.        , 0.49603731],
           [1.        , 0.01205083],
           [1.        , 1.0129211 ],
           [1.        , 1.57774042],
           [1.        , 0.83671178],
           [1.        , 0.35653958],
           [1.        , 1.00719932],
           [1.        , 1.10596479],
           [1.        , 1.11281235],
           [1.        , 1.19792334],
           [1.        , 0.46265133],
           [1.        , 1.81943243],
           [1.        , 0.77343146]])



Normal Equation    

$\hat{\boldsymbol\beta} = (\mathbf{X}^\mathsf{T}\mathbf{X})^{-1} \mathbf{X}^\mathsf{T} \mathbf{y}$


```python
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
theta_best
```




    array([[4.2682626 ],
           [2.78901042]])



<h5>This is close to our real thetas 4 and 3. It cannot be accurate due to the noise I have introduced in data


```python
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
y_predict = X_new_b.dot(theta_best)
y_predict
```




    array([[4.2682626 ],
           [9.84628344]])



<h5>Let's plot prediction line with calculated:theta


```python
plt.plot(X_new, y_predict, "r-")
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
```

# Gradient Descent

## Cost Function & Gradients

<h4> The equation for calculating cost function and gradients are as shown below. Please note the cost function is for Linear regression. For other algorithms the cost function will be different and the gradients would have  to be derived from the cost functions



<b>Cost</b>
\begin{equation}
J(\theta) = 1/2m \sum_{i=1}^{m} (h(\theta)^{(i)} - y^{(i)})^2 
\end{equation}

<b>Gradient</b>

\begin{equation}
\frac{\partial J(\theta)}{\partial \theta_j} = 1/m\sum_{i=1}^{m}(h(\theta^{(i)} - y^{(i)}).X_j^{(i)}
\end{equation}

<b>Gradients</b>
\begin{equation}
\theta_0: = \theta_0 -\alpha . (1/m .\sum_{i=1}^{m}(h(\theta^{(i)} - y^{(i)}).X_0^{(i)})
\end{equation}
\begin{equation}
\theta_1: = \theta_1 -\alpha . (1/m .\sum_{i=1}^{m}(h(\theta^{(i)} - y^{(i)}).X_1^{(i)})
\end{equation}
\begin{equation}
\theta_2: = \theta_2 -\alpha . (1/m .\sum_{i=1}^{m}(h(\theta^{(i)} - y^{(i)}).X_2^{(i)})
\end{equation}

\begin{equation}
\theta_j: = \theta_j -\alpha . (1/m .\sum_{i=1}^{m}(h(\theta^{(i)} - y^{(i)}).X_j^{(i)})
\end{equation}


```python
def cal_cost(theta, X, y):
    """
    
    Calculates the cost for given X and Y. The following shows and example of a single dimensional X
    theta = Vector of thetas 
    X     = Row of X's np.zeros((2,j))
    y     = Actual y's np.zeros((2,1))
    
    where:
        j is the no of features
    """

    m = len(y)

    predictions = X.dot(theta)
    cost = (1 / 2 * m) * np.sum(np.square(predictions - y))
    return cost
```


```python
def gradient_descent(X, y, theta, learning_rate=0.01, iterations=100):
    """
    X    = Matrix of X with added bias units
    y    = Vector of Y
    theta=Vector of thetas np.random.randn(j,1)
    learning_rate 
    iterations = no of iterations
    
    Returns the final theta vector and array of cost history over no of iterations
    """
    m = len(y)
    cost_history = np.zeros(iterations)
    theta_history = np.zeros((iterations, 2))
    for it in range(iterations):

        prediction = np.dot(X, theta)

        theta = theta - (1 / m) * learning_rate * (X.T.dot((prediction - y)))
        theta_history[it, :] = theta.T
        cost_history[it] = cal_cost(theta, X, y)

    return theta, cost_history, theta_history
```

<h3> Let's start with 1000 iterations and a learning rate of 0.01. Start with theta from a Gaussian distribution


```python
lr = 1e-2
n_iter = 1000

theta = np.random.randn(2, 1)

X_b = np.c_[np.ones((len(X), 1)), X]
theta, cost_history, theta_history = gradient_descent(X_b, y, theta, lr, n_iter)


print(
    "Theta0:          {:0.3f},\nTheta1:          {:0.3f}".format(
        theta[0][0], theta[1][0]
    )
)
print("Final cost/MSE:  {:0.3f}".format(cost_history[-1]))
```

    Theta0:          3.969,
    Theta1:          3.033
    Final cost/MSE:  3869.407


<h3> Let's plot the cost history over iterations


```python
fig, ax = plt.subplots(figsize=(12, 8))

ax.set_ylabel("J(Theta)")
ax.set_xlabel("Iterations")
_ = ax.plot(range(n_iter), cost_history, "b.")
```

<h3> After around 150 iterations the cost is flat so the remaining iterations  are not needed or will not result in any further optimization. Let us zoom in till iteration 200 and see the curve


```python
fig, ax = plt.subplots(figsize=(10, 8))
_ = ax.plot(range(200), cost_history[:200], "b.")
```

<b>It is worth while to note that the cost drops faster initially and then the gain in cost reduction is not as much

### It would be great to see the effect of different learning rates and iterations together

### Let us  build a function which can show the effects together and also show how gradient decent actually is working


```python
def plot_GD(n_iter, lr, ax, ax1=None):
    """
     n_iter = no of iterations
     lr = Learning Rate
     ax = Axis to plot the Gradient Descent
     ax1 = Axis to plot cost_history vs Iterations plot

     """
    _ = ax.plot(X, y, "b.")
    theta = np.random.randn(2, 1)

    tr = 0.1
    cost_history = np.zeros(n_iter)
    for i in range(n_iter):
        pred_prev = X_b.dot(theta)
        theta, h, _ = gradient_descent(X_b, y, theta, lr, 1)
        pred = X_b.dot(theta)

        cost_history[i] = h[0]

        if i % 25 == 0:
            _ = ax.plot(X, pred, "r-", alpha=tr)
            if tr < 0.8:
                tr = tr + 0.2
    if not ax1 == None:
        _ = ax1.plot(range(n_iter), cost_history, "b.")
```

### Plot the graphs for different iterations and learning rates combination


```python
fig = plt.figure(figsize=(30, 25), dpi=200)
fig.subplots_adjust(hspace=0.4, wspace=0.4)

it_lr = [(2000, 0.001), (500, 0.01), (200, 0.05), (100, 0.1)]
count = 0
for n_iter, lr in it_lr:
    count += 1

    ax = fig.add_subplot(4, 2, count)
    count += 1

    ax1 = fig.add_subplot(4, 2, count)

    ax.set_title("lr:{}".format(lr))
    ax1.set_title("Iterations:{}".format(n_iter))
    plot_GD(n_iter, lr, ax, ax1)
```

<b> See how useful it is to visualize the effect of learning rates and iterations on gradient descent. The red lines show how the gradient descent starts and then slowly gets closer to the final value

## You can always plot Indiviual graphs to zoom in


```python
_, ax = plt.subplots(figsize=(14, 10))
plot_GD(100, 0.1, ax)
```

# Stochastic Gradient Descent


```python
def stocashtic_gradient_descent(X, y, theta, learning_rate=0.01, iterations=10):
    """
    X    = Matrix of X with added bias units
    y    = Vector of Y
    theta=Vector of thetas np.random.randn(j,1)
    learning_rate 
    iterations = no of iterations
    
    Returns the final theta vector and array of cost history over no of iterations
    """
    m = len(y)
    cost_history = np.zeros(iterations)

    for it in range(iterations):
        cost = 0.0
        for i in range(m):
            rand_ind = np.random.randint(0, m)
            X_i = X[rand_ind, :].reshape(1, X.shape[1])
            y_i = y[rand_ind].reshape(1, 1)
            prediction = np.dot(X_i, theta)

            theta = theta - (1 / m) * learning_rate * (X_i.T.dot((prediction - y_i)))
            cost += cal_cost(theta, X_i, y_i)
        cost_history[it] = cost

    return theta, cost_history
```


```python
lr = 0.5
n_iter = 50

theta = np.random.randn(2, 1)

X_b = np.c_[np.ones((len(X), 1)), X]
theta, cost_history = stocashtic_gradient_descent(X_b, y, theta, lr, n_iter)


print(
    "Theta0:          {:0.3f},\nTheta1:          {:0.3f}".format(
        theta[0][0], theta[1][0]
    )
)
print("Final cost/MSE:  {:0.3f}".format(cost_history[-1]))
```

    Theta0:          4.160,
    Theta1:          2.860
    Final cost/MSE:  38.051



```python
fig, ax = plt.subplots(figsize=(10, 8))

ax.set_ylabel("{J(Theta)}", rotation=0)
ax.set_xlabel("{Iterations}")
theta = np.random.randn(2, 1)

_ = ax.plot(range(n_iter), cost_history, "b.")
```

# Mini Batch Gradient Descent


```python
def minibatch_gradient_descent(
    X, y, theta, learning_rate=0.01, iterations=10, batch_size=20
):
    """
    X    = Matrix of X without added bias units
    y    = Vector of Y
    theta=Vector of thetas np.random.randn(j,1)
    learning_rate 
    iterations = no of iterations
    
    Returns the final theta vector and array of cost history over no of iterations
    """
    m = len(y)
    cost_history = np.zeros(iterations)
    n_batches = int(m / batch_size)

    for it in range(iterations):
        cost = 0.0
        indices = np.random.permutation(m)
        X = X[indices]
        y = y[indices]
        for i in range(0, m, batch_size):
            X_i = X[i : i + batch_size]
            y_i = y[i : i + batch_size]

            X_i = np.c_[np.ones(len(X_i)), X_i]

            prediction = np.dot(X_i, theta)

            theta = theta - (1 / m) * learning_rate * (X_i.T.dot((prediction - y_i)))
            cost += cal_cost(theta, X_i, y_i)
        cost_history[it] = cost

    return theta, cost_history
```


```python
lr = 0.1
n_iter = 200

theta = np.random.randn(2, 1)


theta, cost_history = minibatch_gradient_descent(X, y, theta, lr, n_iter)


print(
    "Theta0:          {:0.3f},\nTheta1:          {:0.3f}".format(
        theta[0][0], theta[1][0]
    )
)
print("Final cost/MSE:  {:0.3f}".format(cost_history[-1]))
```

    Theta0:          4.201,
    Theta1:          2.844
    Final cost/MSE:  750.136



```python
fig, ax = plt.subplots(figsize=(10, 8))

ax.set_ylabel("{J(Theta)}", rotation=0)
ax.set_xlabel("{Iterations}")
theta = np.random.randn(2, 1)

_ = ax.plot(range(n_iter), cost_history, "b.")
```

9d10ba27000b3b9d1850bd00415108558bbb0ba1
ipynb
Jupyter Notebook
content/week-05/samples/calculus-02/Gradient Descent (NumPy).ipynb
GiorgiBeriashvili/school-of-ai
abd033fecf32c1222da097aa8420db6c69b357e6
Apache-2.0 MIT
content/week-05/samples/calculus-02/Gradient Descent (NumPy).ipynb
GiorgiBeriashvili/school-of-ai
abd033fecf32c1222da097aa8420db6c69b357e6
Apache-2.0 MIT
content/week-05/samples/calculus-02/Gradient Descent (NumPy).ipynb
GiorgiBeriashvili/school-of-ai
abd033fecf32c1222da097aa8420db6c69b357e6
Apache-2.0 MIT
Qwen/Qwen-72B
1. YES
2. YES