<jupyter_start><jupyter_text>Linear RegressionIn this I python notebook we will focus on $\textit{Linear Regression}$ task. We will show how to prove it mathematically and how to create it in numpy (only with basic matrix operations).  ProblemFor a given set $D = \{(x_n, y_n)\}_{n=1}^N$ predict the $\textit{y}$ for the new $\textbf{x}$. To do it we will have to find function that best approximates the training set.  Example Flat size of $x_n$ has the price $y_n$. We have the following data:      Flat size $[m^2]$    Price $[k \$]$         20    100         30    180           55    320            80    450   We want to predict price $\textit{y}$ for the new flat $\textbf{x}$ based on its size.      Flat size $[m^2]$    Price $[k \$]$         45    ?         70    ?<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import random<jupyter_output><empty_output><jupyter_text>Design matrix $\Phi$Design matrix $\Phi$ is matrix of values of explanatory variables of a set of objects. Depending on the definition it takes different forms f.e:$\phi(x) = \begin{pmatrix}1 & x & x^2 \\\end{pmatrix}^T$$\phi(x) = \begin{pmatrix}1 & e^x & e^{2x} \\\end{pmatrix}^T$Where $x$ is the vector of observations $x_n$ Example: $ x = \begin{bmatrix}           1 \\ 2 \\ 3          \end{bmatrix} $$\phi(x) = \begin{pmatrix}1 & x & x^2 \\\end{pmatrix}^T \quad \; \; \Phi=\begin{bmatrix}       1 & 1 & 1 \\       1 & 2 & 4 \\       1 & 3 & 9 \\     \end{bmatrix}$ $\phi(x) = \begin{pmatrix}1 & e^x & e^{2x} \\\end {pmatrix}^T \quad \Phi=\begin{bmatrix}       1 & e^1 & e^2 \\       1 & e^2 & e^4 \\       1 & e^3 & e^6 \\     \end{bmatrix}$     In our example I would use only first form of polynomial matrix<jupyter_code>def design_matrix(x_train, M):
    '''
    :param x_train: vector of input values Nx1
    :param M: degree of polynomial 0,1,2,...
    :return: Design Matrix X_train for polynomial degree of M
    '''
    matrix = np.array([x_train ** i for i in range(M + 1)]).T
    
    return matrix[0]

a = np.array([1,2,3])
a = np.reshape(a, (3,1))

print("Input vector: \n", a)
print()
print ("Design matrix: \n", design_matrix(a, 4))
print ("Design shape: \n", design_matrix(a, 4).shape)<jupyter_output>Input vector: 
 [[1]
 [2]
 [3]]

Design matrix: 
 [[ 1  1  1  1  1]
 [ 1  2  4  8 16]
 [ 1  3  9 27 81]]
Design shape: 
 (3, 5)<jupyter_text>Model$\textbf{w} \in R^{Dx1} $ is a vector of parameters$\Phi \in R^{NxD}$ is a design matrix$\overline{y} = \Phi \textbf{w} \in R^{Nx1}$ is a vector of model predictionsWe need to create method to compute the polynomial value. This method takes set $X = \{x_n\}_{n=1}^N$ and set $w = \{w_m\}_{m=1}^{M+1}$. Then computs polynomial values for imputs $x_n$Example:$X = \{1, 2, 3, -4\}$Polynomial is quadratic polynomial (Second degree polynomial): $y = c + bx + ax^2$$\textbf{w} = \{1, 3, 2\}$$c = 1$$b = 3$$a = 2$$Result = \{1+3*1+2*1^2, \quad 1+3*2+2*2^2 , \quad 1+3*3+2*3^2 , \quad 1+3*(-4)+2*(-4)^2\} = \{6, 15, 28, 21\}$<jupyter_code>def polynomial(X, w):
    '''
    :param x: vector of input values Nx1
    :param w: vector of prameters (M+1)x1 
    :return: vector of polynomial values for the input x
    '''
    dm = [w[i] * X ** i for i in range(np.shape(w)[0])]
    return np.sum(dm, axis=0)

X = np.array([1,2,3,-4])
w = np.array([1, 3, 2])
p = polynomial(X, w)
print("X:", X)
print()
print("Parameters:", w)
print()
print("Result:", p)<jupyter_output>X: [ 1  2  3 -4]

Parameters: [1 3 2]

Result: [ 6 15 28 21]<jupyter_text>Mean Squared ErrorTo measure how good our model is we would use the cost function called Mean Squered Error$E = \frac{1}{2N}\lVert\overline{\textbf{y}} - \textbf{y}\rVert_2^2$$\lVert \overline{\textbf{y}} - \textbf{y}\rVert_2$ is called Norm 2 In linear algebra, functional analysis, and related areas of mathematics, a norm is a function that assigns a strictly positive length or size to each vector in a vector space$\lVert\textbf{x}\rVert_2 = \sqrt[2]{\sum_{n=1}^{N} x_n^2}$<jupyter_code>def mean_squared_error(y_pred, y_real):
    '''
    :param x: vector of input values Nx1
    :param y: vector of output values Nx1
    :param w: vector of prameters (M+1)x1 
    :return: mean_squared_error between y and predicted_y, where y are actual y values and predicted_y are computed 
    form a polynomial with w parameters
    '''
    N = y_pred.shape[0]

    error = np.sum((y_pred - y_real)**2)/(2*N)

    return error
y_pred = p
y_real = [7, 14, 22, 30]

print("Y predicted", y_pred)
print()
print("Y real", y_real)
print()
print("Mean Squared Error", mean_squared_error(y_pred, y_real))<jupyter_output>Y predicted [ 6 15 28 21]

Y real [7, 14, 22, 30]

Mean Squared Error 14.875<jupyter_text>Minimalize cost functionTo find best parameters $\textbf{w}$ we have to minimize the cost function. As we know from Calculus to find the function miniumum we have to calculate the derivative (called gradient because we have to find best vector $\textbf{w}$ not a single scalar), set it equal to zero and solve for $\textbf{w}$ . $\textbf{N}$ is const. so in the task of minimizing function we can skip it$E(w) = \frac{1}{2}\lVert\widehat{\textbf{y}} - \textbf{y}\rVert_2^2$$\overline{\textbf{y}} = \Phi \textbf{w}$$E(w) = \frac{1}{2}\lVert \Phi \textbf{w} - \textbf{y}\rVert_2^2$Calculate derivative$\nabla_w E(w) = \nabla_w (\frac{1}{2} \lVert\Phi \textbf{w} - \textbf{y}\rVert_2^2) = -2\frac{1}{2} \Phi^T( \textbf{y} - \Phi \textbf{w}) =  -\Phi^T( \textbf{y} - \Phi \textbf{w})$Set it equal to zero$\nabla_w E(w) = 0$$\Phi^T( \textbf{y} - \Phi \textbf{w}) = 0$$\Phi^Ty - \Phi^T\Phi \textbf{w} = 0 $$\Phi^Ty = \Phi^T\Phi \textbf{w} $$(\Phi^T\Phi)^{-1} \Phi^T\textbf{y} = \mathbb{1}\textbf{w}$$\mathbb{1}$ is Identity Matrix$ \widehat{\textbf{w}}= (\Phi^T\Phi)^{-1} \Phi^T\textbf{y} $ExampleLet us assume, that this is our real function:$y = (x+3)(x-2) = -6 + x + x^2$And in our data we have three points:$A = (1, -4), \space  B = (-3, 0), \space C = (-1, -6)$ which ideally overlap with our function y. First of all we need to divde our data, and prepare$\textbf{x} = \begin{bmatrix}1 \\ -3 \\ -1 \end{bmatrix} $ $\space \textbf{y} = \begin{bmatrix}-4 \\ 0 \\ 6 \end{bmatrix}$ and the real parameters of function $\textbf{w}  = \begin{bmatrix} -6 \\ 1 \\ 1 \end{bmatrix} $We can also construct our designed matrix and transpose it:$\Phi = \begin{bmatrix} 1 & 1 & 1 \\ 1 & -3 & 9 \\ 1 & -1 & 1 \end{bmatrix}$$\Phi^T = \begin{bmatrix} 1 & 1 & 1 \\ 1 & -3 & -1 \\ 1 & 9 & 1 \end{bmatrix}$Then we have to compute the parameters with the formula shown above:$\widehat{\textbf{w}}= (\Phi^T\Phi)^{-1} \Phi^T\textbf{y}$$\widehat{\textbf{w}} = \left ( \begin{bmatrix} 1 & 1 & 1 \\ 1 & -3 & -1 \\ 1 & 9 & 1 \end{bmatrix} \times \begin{bmatrix} 1 & 1 & 1 \\ 1 & -3 & 9 \\ 1 & -1 & 1 \end{bmatrix} \right )^{-1}  \times \begin{bmatrix} 1 & 1 & 1 \\ 1 & -3 & -1 \\ 1 & 9 & 1 \end{bmatrix} \times \begin{bmatrix} -4 \\ 0 \\ 6  \end{bmatrix} =  \begin{bmatrix} -6 \\ 1 \\ 1  \end{bmatrix} $  So our estimation $\widehat{\textbf{w}}$  founded with this method is equal to the real $\textbf{w}$ and the error goes to zero.<jupyter_code>def least_squares(x_train, y_train, M):
    '''
   :param x: vector of input values Nx1
   :param y: vector of output values Nx1
   :param M: degree of polynomial
   :return: tuple (w,err), where w are best polynomial parameters and err mean squuared error
   '''
    dm = design_matrix(x_train, M)

    w = np.linalg.inv(dm.T @ dm) @ dm.T @ y_train
    
    y_pred = polynomial(x_train, w)

    err = mean_squared_error(y_pred, y_train)

    err = np.asscalar(err)

    return (w, err)
x_train = np.array([[1], [-3], [-1]])
y_train = np.array([[-4], [0], [-6]])

print("X coordinates of our training data: \n", x_train)
print()
print("Y coordinates of our training data: \n", y_train)
w = least_squares(x_train, y_train, 2)
print()
print("Estimated parameters w: \n", w[0])<jupyter_output>X coordinates of our training data: 
 [[ 1]
 [-3]
 [-1]]

Y coordinates of our training data: 
 [[-4]
 [ 0]
 [-6]]

Estimated parameters w: 
 [[-6.]
 [ 1.]
 [ 1.]]<jupyter_text>Select best modelNow we need to select model, that fits x_val in the best way. To find best polinomial degree, on each degree we will: 1) Using least_squeres method find best polynomial parameters $ \textbf{w}$ for current degree2) Using founded parameters $ \textbf{w}$, count error on validation set3) If error on validation set is smaller then the smallest error, then we change old values to the new parameters as it is the best solution founded by now<jupyter_code>def model_selection(x_train, y_train, x_val, y_val, M_values):
    '''
   :param x_train: vector of input training values Nx1
   :param y_train: vector of output training values Nx1
   :param x_val: vector of input validation values Nx1
   :param y_val: vector of output vallidation values Nx1
   :param M_values: list of polynomial degrees 
   :return: function returns tuple (w,train_err,val_err), w is vector of model prameters, 
   which fits model in best way on validation set
   train_err and val_err are mean squered errors on training and validation sets
   '''

    val_err = np.inf
    for m in M_values:
        current_w, current_train_err = least_squares(x_train, y_train, m)
        y_pred = polynomial(x_val, current_w)
        current_val_err = mean_squared_error(y_pred, y_val)

        if current_val_err < val_err:
            val_err = current_val_err
            w = current_w
            train_err = current_train_err

    val_err = np.asscalar(val_err)

    return (w, train_err, val_err)<jupyter_output><empty_output><jupyter_text>Finding best parameters visualization<jupyter_code>def run():
    def target_output(x):
        return np.sin(2 * np.pi * x)

    def plot_model(subplot, x_train, y_train, x, y_obj, y_model, x_val=None, y_val=None, train_err=None,
               val_err=None):
        x_min = np.min([np.min(x_train), np.min(x)])
        x_max = np.max([np.max(x_train), np.max(x)])
        y_min = -1.5
        y_max = 1.5
        int_x = x_max - x_min
        x_beg = x_min - int_x / 14.0
        x_end = x_max + int_x / 14.0
        x_ticks = [x_min, x_max]
        int_y = y_max - y_min
        y_ticks = [y_min, y_min + 0.5 * int_y, y_max]

        sub.set_xlim(x_beg, x_end)
        sub.set_ylim(1.1 * y_min, 1.1 * y_max)
        sub.set_xticks(x_ticks)
        sub.set_yticks(y_ticks)
        sub.plot(x_train, y_train, 'o', markerfacecolor='none', markeredgecolor='blue', markersize=8,
                 markeredgewidth=2)
        sub.plot(x, y_obj, '-g', linewidth=2)
        sub.plot(x, y_model, '-r', linewidth=2)
        if x_val is not None and y_val is not None:
            sub.plot(x_val, y_val, 'o', markerfacecolor='none', markeredgecolor='red', markersize=8,
                     markeredgewidth=2)
        if train_err is not None and val_err is not None:
            sub.text(0, -1.3, 'Train error: {0:.5f}\nVal error:    {1:.5f}'.format(train_err, val_err),
                     bbox={'facecolor': 'none', 'pad': 10})

    
    
    def create_data():
        np.random.seed(1)
        x_train_8 = np.random.sample((8,1))*1.5
        x_train_50 = np.random.rand(50,1)
        x_val_10 = np.random.rand(10,1)
        x_val_20 = np.random.rand(20,1)
        y_train_8 = target_output(x_train_8) + np.random.randn(8,1)/3
        y_train_50 = target_output(x_train_50) + np.random.randn(50,1)/3
        y_val_10 = target_output(x_val_10) + np.random.randn(10,1)/3
        y_val_20 = target_output(x_val_20) + np.random.randn(20,1)/3
        
        created_data = {'x_train_8': x_train_8, 'x_train_50': x_train_50, 'x_val_10': x_val_10, 'x_val_20': x_val_20, 
                        'y_train_8': y_train_8, 'y_train_50': y_train_50, 'y_val_10': y_val_10, 'y_val_20': y_val_20}
        return created_data
    data = create_data()#pickle.load(open('LinearRegression/data.pkl', mode='rb'))
    
    x_plot = np.arange(0, 1.01, 0.01)
    y_obj = target_output(x_plot)

    # Dopasowanie wielomianow metoda najmniejszych kwadratow
    print('\n---Mean ---')
    print('-------------- Training examples amount N=8. --------------')
    fig = plt.figure(figsize=(12, 6), num='Least squuares task for N=8')

    for i in range(8):
        w, err = least_squares(data['x_train_8'], data['y_train_8'], i)
        y_model = polynomial(x_plot, w)
        sub = fig.add_subplot(2, 4, i + 1)
        plot_model(sub, data['x_train_8'], data['y_train_8'], x_plot, y_obj, y_model)
        sub.set_title("M = {}".format(i))

    plt.show()

    print('\n--- Fitting polynomials with least squares method ---')
    print('-------------- Training examples amount N=50. --------------')
    fig = plt.figure(figsize=(12, 6), num='Least squares task for N=50')

    for i in range(8):
        w, err = least_squares(data['x_train_50'], data['y_train_50'], i)
        y_model = polynomial(x_plot, w)
        sub = fig.add_subplot(2, 4, i + 1)
        plot_model(sub, data['x_train_50'], data['y_train_50'], x_plot, y_obj, y_model)
        sub.set_title("M = {}".format(i))

    plt.show()

    # Model selection
    print('\n--- Model Selection ---')
    print('---------------- Polynomial models degree M=0,...,7 ----------------')
    print('- Data training points number N=50. Validation data points number N=20 -')

    M_values = range(0, 7)
    w, train_err, val_err = model_selection(data['x_train_50'], data['y_train_50'],
                                            data['x_val_20'], data['y_val_20'], M_values)
    M = np.shape(w)[0] - 1
    y_model = polynomial(x_plot, w)

    fig = plt.figure(figsize=(6, 5), num='Selection model for M')
    sub = fig.add_subplot(1, 1, 1)
    sub.set_title('The best M={}'.format(M))
    plot_model(sub, data['x_train_50'], data['y_train_50'], x_plot, y_obj, y_model,
               data['x_val_20'], data['y_val_20'], train_err, val_err)

    plt.show()
    
run()<jupyter_output>---Mean ---
-------------- Training examples amount N=8. --------------<jupyter_text>RegularizationThe goal of regularization is to avoid overfitting, in other words we are trying to avoid models that fit extremely well to the training data (data used to build the model), but fit poorly to testing data (data used to test how good the model is). When we fit a model we try to minimize a cost function (a measure of error). We can get really good and bring the cost function to zero. However, when we generalize the model (e.g. use it on new data that has not been seen before), the model could not work very well. This is known as overfitting. To avoid overfitting, we include in our cost function a penalty term, known as regularization, such that we want the error to be small but the sum of the absolute value of the parameters to also be small$\lambda$ is a hyperparameter called lambda that we need to choose on our own. To find best model it is common to chceck errors for a few lambda valuesRegularization using squared norm 2 is called L2 Regularization$E(w) = \frac{1}{2}\lVert\overline{\textbf{y}} - \textbf{y}\rVert_2^2 + \frac{\lambda}{2} \lVert\textbf{w}\rVert^2_2$$\overline{\textbf{y}} = \Phi \textbf{w}$$E(w) = \frac{1}{2}\lVert \Phi \textbf{w} - \textbf{y}\rVert_2^2 + \frac{\lambda}{2} \lVert\textbf{w}\rVert^2_2$Calculate derivative$\nabla_w E(w) = \nabla_w (\frac{1}{2} \lVert \Phi \textbf{w} - \textbf{y}\rVert_2^2 + \frac{\lambda}{2}\lVert\textbf{w}\rVert^2_2) = -2\frac{1}{2} \Phi^T( \textbf{y} - \Phi \textbf{w}) + \frac{\lambda}{2} 2\textbf{w}=-\Phi^T( \textbf{y} - \Phi\textbf{w}) + \lambda\textbf{w}$Set it equal to zero$\nabla_w E(w) = 0$$-\Phi^T( \textbf{y} - \Phi \textbf{w}) + \lambda\textbf{w}= 0$$\Phi^T y - \Phi^T \Phi \textbf{w} - \lambda\textbf{w}= 0 $$\Phi^Ty = \Phi^T\Phi \textbf{w} + \lambda\textbf{w} $$\Phi^Ty = (\Phi^T\Phi + \lambda \mathbb{1}) \textbf{w} $$\mathbb{1}$ is Identity Matrix$ \widehat{\textbf{w}}= (\Phi^T\Phi + \lambda \mathbb{1})^{-1} \Phi^T\textbf{y}$<jupyter_code>def regularized_least_squares(x_train, y_train, M, regularization_lambda):
    '''
   :param x_train: train input data Nx1
   :param y_train: train input target Nx1
   :param M: polynomial degree
   :param regularization_lambda: regularization parameter
   :return: function returns tuple (w, err), where w are parameters of polynomial fit according to the regularization l2 criterium
   and err is a mean squared error.
   '''

    dm = design_matrix(x_train, M)

    w = np.linalg.inv(dm.T @ dm + regularization_lambda * np.eye(M+1)) @ dm.T @ y_train
    
    y_pred = polynomial(x_train, w)

    err = mean_squared_error(y_pred, y_train)

    err = np.asscalar(err)

    return (w, err)<jupyter_output><empty_output><jupyter_text>Select best model with regularizationIn function below we assume that we know the polynomial degree, so we will look only for best lambda value. We will iterate over lambdas list and choose the best one.<jupyter_code>def regularized_model_selection(x_train, y_train, x_val, y_val, M, lambda_values):
    '''
    :param x_train: vector of input training values Nx1
   :param y_train: vector of output training values Nx1
   :param x_val: vector of input validation values Nx1
   :param y_val: vector of output vallidation values Nx1
   :param D: degree of polynomial
   :param lambda_values: list of diffrent lambda values
   :return: function return tuple (w,train_err,val_err,regularization_lambda), where w are parameters of polynomial which generalizes 
    in the best way -that gives least error on validation set.
    Polynomial is adjusted to the best lambda parameter - lambda_value is the best found lambda.
   '''
    val_err = np.inf

    for l in lambda_values:
        current_w, current_train_err = regularized_least_squares(x_train, y_train, M, l)
        y_pred = polynomial(x_val, current_w)
        current_val_err = mean_squared_error(y_pred, y_val)

        if current_val_err < val_err:
            val_err = current_val_err
            w = current_w
            train_err = current_train_err
            lambda_value = l

    val_err = np.asscalar(val_err)

    return (w, train_err, val_err, lambda_value)
def run_regularized():
    
    def target_output(x):
        return np.sin(2 * np.pi * x)

    def plot_model(subplot, x_train, y_train, x, y_obj, y_model, x_val=None, y_val=None, train_err=None,
               val_err=None):
        x_min = np.min([np.min(x_train), np.min(x)])
        x_max = np.max([np.max(x_train), np.max(x)])
        y_min = -1.5
        y_max = 1.5
        int_x = x_max - x_min
        x_beg = x_min - int_x / 14.0
        x_end = x_max + int_x / 14.0
        x_ticks = [x_min, x_max]
        int_y = y_max - y_min
        y_ticks = [y_min, y_min + 0.5 * int_y, y_max]

        sub.set_xlim(x_beg, x_end)
        sub.set_ylim(1.1 * y_min, 1.1 * y_max)
        sub.set_xticks(x_ticks)
        sub.set_yticks(y_ticks)
        sub.plot(x_train, y_train, 'o', markerfacecolor='none', markeredgecolor='blue', markersize=8,
                 markeredgewidth=2)
        sub.plot(x, y_obj, '-g', linewidth=2)
        sub.plot(x, y_model, '-r', linewidth=2)
        if x_val is not None and y_val is not None:
            sub.plot(x_val, y_val, 'o', markerfacecolor='none', markeredgecolor='red', markersize=8,
                     markeredgewidth=2)
        if train_err is not None and val_err is not None:
            sub.text(0, -1.3, 'Train error: {0:.5f}\nVal error:    {1:.5f}'.format(train_err, val_err),
                     bbox={'facecolor': 'none', 'pad': 10})
    
    
    def create_data():
        np.random.seed(1)
        x_train_8 = np.random.sample((8,1))*1.5
        x_train_50 = np.random.rand(50,1)
        x_val_10 = np.random.rand(10,1)
        x_val_20 = np.random.rand(20,1)
        y_train_8 = target_output(x_train_8) + np.random.randn(8,1)/3
        y_train_50 = target_output(x_train_50) + np.random.randn(50,1)/3
        y_val_10 = target_output(x_val_10) + np.random.randn(10,1)/3
        y_val_20 = target_output(x_val_20) + np.random.randn(20,1)/3
        
        created_data = {'x_train_8': x_train_8, 'x_train_50': x_train_50, 'x_val_10': x_val_10, 'x_val_20': x_val_20, 
                        'y_train_8': y_train_8, 'y_train_50': y_train_50, 'y_val_10': y_val_10, 'y_val_20': y_val_20}
        return created_data
    data = create_data()
    
    x_plot = np.arange(0, 1.01, 0.01)
    y_obj = target_output(x_plot)
    
    print('\n--- Model selection for linear task of least squares error with regularization ---')
    print('-- Degree M=7. Data training points N=50. Data validation points N=20 --')

    M = 3
    lambdas = [0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100, 300]
    w, train_err, val_err, best_lambda = regularized_model_selection(data['x_train_50'],
                                                                     data['y_train_50'],
                                                                     data['x_val_20'],
                                                                     data['y_val_20'],
                                                                     M, lambdas)
    y_model = polynomial(x_plot, w)

    fig = plt.figure(figsize=(6, 5), num='Model selection for regularization parameter')
    sub = fig.add_subplot(1, 1, 1)
    sub.set_title('M={}    Best $\lambda$={}'.format(M, best_lambda))
    plot_model(sub, data['x_train_50'], data['y_train_50'], x_plot, y_obj, y_model,
               data['x_val_20'], data['y_val_20'],
               train_err, val_err)

    plt.show()

run_regularized()<jupyter_output>--- Model selection for linear task of least squares error with regularization ---
-- Degree M=7. Data training points N=50. Data validation points N=20 --<jupyter_text>Gradient DescentTo say the truth - method to compute parameters explained above isn't computionally effective. We can even say, that the computing cost is terrible, aspeccialy for larger data matrixes, and things goes even worse and worse when the number of data dimensions increases. Inversing a matrix has $n^2\cdot\log{n}$ complexity.There is a way more servicable solution to compute optimal parameters. The method we are going to describe is called Gradient Descent.How does it work?At first we initialize random parameters, and multiply them properly with inputs. Then we are computing the error between those random predictions and target values. As we got the error, we can compute the derivative of error of the parameters at the exact point we are. Why are we doing this? Derivative give us information about the way in which function changes. With information like that we can substract the value of derivative from parameters, as there we got two possible options:1. The function grows, so the derivative is positive, and when we are substracting it we are "going back" to the local minimum.2. The function decreases, so the derivative is negative, and when we are substracting it we are "going forward" to the local minimum. It looks beautiful, when we are working with low-dimensional data, as for example imagine yourself that you have $y=x^2$ function, and you are at point (2, 4). Derivative is positive, so after substracting it ( of course after multiplication it by a small parameter called learning rate ) we are moving closer to zero - the minimum.However, when we have much more dimensions - it is hard to create a cost function that has only one minimum. In that case we are exposed for stucking in some of the local minimum, but there are a few ways of avoiding this, which are not described in this notebook. Learning rateThe learning rate parameter is a hyperparameter that defines "how fast" our parameters are going to be learned, as we are multiplying value of derivatie by it. We need to be cautios when choosing it, as too large learning rate can make us to miss the local minimum. We wouldn't be able to go to "zero" in the $y = x^2$ function we can jump from point (2, 4) to (-2, 4) and back, when it is going to be too large. On the other side choosing a very small learning rate can make our program to find optimal parameters even slower than using the formula which involves matrix inversing. TL;DR - Gradient Descent in steps1. Define Cost function Define the cost function to measure the diffrence between predictions and target outputs. We assume that we are working with first degree polynomial, to make derivative calculations easier to present. ( Linear function y = wx +b )$$Error = \frac{1}{N}\sum_{i=1}^N (y_i - \overline{y}_i)^2 = \frac{1}{N}\sum_{i=1}^N (y_i - (x_iw+b))^2 $$Where N is an observations amount2. Compute the derivative$$\frac{\delta Error}{\delta w} = \frac{2}{N}\sum_{i=1}^N -x_i(y_i -(m  x_i +b ))  $$$$\frac{\delta Error}{\delta b} = \frac{2}{N}\sum_{i=1}^N -(y_i -(m  x_i +b ))  $$3. Update current parameters$$ w:= w- learning\_rate \cdot \frac{\delta Error}{\delta w}   $$ $$ b:= b- learning\_rate \cdot \frac{\delta Error}{\delta b}   $$ 4. Repeat until it fits good enough<jupyter_code>def gradient_descent(x, y_expected, epochs=1000, learning_rate=0.00001):
    '''   
    :param x: vector of input training values Nx1
    :param y_expected: vector of output training values Nx1
    :param epochs: amount of training epochs
    :param learning_rate: rate of learning speed
    '''
    w = random.uniform(0, 1)
    b = random.uniform(0, 1)

    N = x.shape[0]
    cost = []
    for i in range(epochs):
        y_current = w * x + b
        cost.append(np.sum((y_expected - y_current) ** 2) / N)
        w_gradient = - 2 / N * np.sum(2 * (y_expected - y_current) * x)
        b_gradient = - 2 / N * np.sum(2 * (y_expected - y_current))

        w = w - (learning_rate * w_gradient)
        b = b - (learning_rate * b_gradient)

    return w, b, cost
x = np.array([[i for i in range(10)]]).T
y = np.array([[5 * i + 10 + random.uniform(-.1, .1) for i in range(10)]]).T
print("W expected: {} and b expected: {}".format(5, 10))
print()
w, b, _ = gradient_descent(x,y)
print("W calculated: {:.2f} and b calculated: {:.2f}".format(w, b))<jupyter_output>W expected: 5 and b expected: 10

W calculated: 4.47 and b calculated: 1.47<jupyter_text>Real world exampleWe are going to predict student score according to the amount of hours they studied.<jupyter_code>def real_world_example():
    points = np.genfromtxt("LinearRegression/data.txt", delimiter=',')
    print('Some example hours and scores from exam:')
    print(points[:5])
    x, y = points[0], points[1]
    w, b, costs = gradient_descent(x, y, epochs=30000, learning_rate=1e-8)
    plt.plot(costs)
    plt.title('Cost')
    plt.xlabel('Epoch number')
    plt.ylabel('Cost value')
    plt.grid()
    plt.show()
    x_plot = np.arange(20, 70, 0.5)
    y_plot = x_plot*w+b
    
    hours = points[:,0]
    scores = points[:,1]
    plt.plot(hours, scores, '.', x_plot, y_plot, '-')
    plt.xlabel('Hours of study')
    plt.ylabel('Scores on exam')
    plt.legend(['Score points', 'Linear regression'])
    plt.grid()
    
    plt.show()

real_world_example()<jupyter_output>Some example hours and scores from exam:
[[32.50234527 31.70700585]
 [53.42680403 68.77759598]
 [61.53035803 62.5623823 ]
 [47.47563963 71.54663223]
 [59.81320787 87.23092513]]
{'fraction_non_alphanumeric': 0.13877235592846102, 'fraction_numerical': 0.03378229189666593, 'mean_word_length': 3.048115596603605, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '154438', 'license': 'MIT', 'n_tokens_mistral': 9877, 'n_tokens_neox': 9219, 'n_words': 2967, 'path': 'LinearRegression.ipynb', 'repo_name': 'tugot17/Linear-Regression-From-Scratch', 'size': '27174'}