<jupyter_start><jupyter_text><NAME> > Write a code in Python whose input is a training dataset {(x1, y1), . . . , (xN , yN )} and its output is the weight vector θ in the linear regression model y = θ'φ(x), for a given nonlinear mapping φ(·). Implement two cases: i) using the closed-form solution, ii) using a stochastic gradient descent on mini-batches of size m.<jupyter_code>import numpy as np
from numpy.linalg import inv
import scipy.io
import math
import pandas as pd
import time

# For reproducibility
np.random.seed(1000)

# Load the dataset
df = scipy.io.loadmat("dataset1.mat")
df2 = scipy.io.loadmat("dataset2.mat")
# This function generates data based on the formula we input 
# For this time we will generate a square function added with a noise
f = lambda  x: x**2
# Closed function regression
def closed_function_regression(train_X, train_Y):
  X = np.array(train_X)
  y = np.array(train_Y)
  
        
  Xt = X.T
  product = np.dot(Xt, X)
  theInverse = inv(product)
  theta = np.dot(np.dot(theInverse, Xt), y)
        
  return theta  
# Source of this function: http://adventuresinmachinelearning.com/stochastic-gradient-descent/
def shuffle_data(train_X, train_Y, batch_size):
    X = train_X
    y = train_Y
    
    random_indices = np.random.choice(len(X), len(y), replace=False)
    
    X_shuffled = X[random_indices,:]
    y_shuffled = y[random_indices]
    mini_batches = [(X_shuffled[i:i+batch_size,:], y_shuffled[i:i+batch_size]) for
                   i in range(0, len(y), batch_size)]
    
    return mini_batches
def mini_gradient_descent(train_X, train_Y, learning_rate, num_iter, batch_size):
    # Prepare the data
    X = np.array(train_X)
    y = np.array(train_Y)
    
    # Randomize the data and split into samples
    df = shuffle_data(X, y, batch_size)
    
    # Get the dimension
    x_col = X.shape[1]
    
    # Set the initial theta and other initial variables
    theta = np.zeros((x_col, 1))
    
    start_i = 0
    
    # Create the loop
    for i in range(start_i + 1, num_iter + 1):
        # For each randomized-mini-batch partition
        for j in range(0, len(df)):
            X = df[j][0]
            y = df[j][1]
            
            y_hat = np.dot(X, theta)
            
            # Break if overflow occurs
            theta = theta - learning_rate * np.dot(X.T, y_hat - y)
            
            # FIXME: Break if NaN occurs
            # if not (float('-inf') < float(temp[0]) < float('inf')): 
            #    return theta
                
            # theta =  temp
            
        
    return theta<jupyter_output><empty_output><jupyter_text>> Consider n-degree polynomials, φ(·) =  1 x x^2 · · · x^n . Work with ‘dataset1’. Run the code on the training data to compute θ for n ∈ {2, 3, 5}. Evaluate the regression error on both training and the test data. Report θ, training error and test error for both implementation (closed-form vs gradient descent). What is the effect of the size of the mini-batch on the speed and testing error of the solution.<jupyter_code>X_trn = df['X_trn']
Y_trn = df['Y_trn']
X_test = df['X_tst']
Y_test = df['Y_tst']


closed_function_regression(X_trn, Y_trn)
def create_polinomial_array(element, n):
    pol = np.empty([1, n+1])
    # Debug
    # print(pol)
    for i in range(0, n+1):
        # Debug
        # print("Element: \n", element)
        pol[0][i] = element ** i
        
    return pol[0]


# Here we set the degree of polynomial        
degree = 5
# degree + 1 for 1's in the first column
X_trn_new = np.ndarray((len(X_trn), degree + 1))

# print(X_trn_new)

i = 0
X = np.array(X_trn)
for x in X:
    X_trn_new[i] = create_polinomial_array(x, degree)
    i += 1
        

mini_gradient_descent(X_trn, Y_trn, 0.001, 100, 10)
# Analysing the regression error and speed
# Hyper-parameter are the alpha (learning rate), the batch size and the number of iterations
def analyze(learning_rate, n_iteration, batch_size):
    degrees = [2, 3, 5]
    
    # Time for performance
    grad_desc_performance = []
    
    # For train data
    error_closed_form_train = []
    error_grad_desc_train = []
    
    # For test data
    error_closed_form_test = []
    error_grad_desc_test = []
    
    X = np.array(X_trn)
    X_tst = np.array(X_test)
    
    for degree in degrees:
        X_trn_new = np.ndarray((len(X_trn), degree + 1))
        X_tst_new = np.ndarray((len(X_tst), degree + 1))
        
        i = 0
        for x in X:
            X_trn_new[i] = create_polinomial_array(x, degree)
            i += 1
        
        i = 0
        for x in X_tst:
            X_tst_new[i] = create_polinomial_array(x, degree)
            i += 1
        
        tetha_closed_form = closed_function_regression(X_trn_new, Y_trn)
        t = time.process_time()
        tetha_grad_desc = mini_gradient_descent(X_trn_new, Y_trn, learning_rate, n_iteration, batch_size)
        elapsed_time = time.process_time() - t
        print("Theta for degree {} done in: {}\n".format(degree, elapsed_time))
        
        
        y_hat_closed_form_train = np.dot(X_trn_new, tetha_closed_form)
        y_hat_grad_desc_train = np.dot(X_trn_new, tetha_grad_desc)
        
        y_hat_closed_form_test = np.dot(X_tst_new, tetha_closed_form)
        y_hat_grad_desc_test = np.dot(X_tst_new, tetha_grad_desc)
        
        error_closed_form_train.append(np.sqrt(((Y_trn - y_hat_closed_form_train) ** 2).mean()))
        error_grad_desc_train.append(np.sqrt(((Y_trn - y_hat_grad_desc_train) ** 2).mean()))
        
        error_closed_form_test.append(np.sqrt(((Y_test - y_hat_closed_form_test) ** 2).mean()))
        error_grad_desc_test.append(np.sqrt(((Y_test - y_hat_grad_desc_test) ** 2).mean()))
        
    d = {
        'degree': [2, 3, 5],
        'close_form_error_train': error_closed_form_train,
        'grad_desc_error_train': error_grad_desc_train,
        'close_form_error_test': error_closed_form_test,
        'grad_desc_error_test': error_grad_desc_test
    }
    errors = pd.DataFrame(data = d)
    return errors



# Batch size 10
result1 = analyze(0.00001, 10, 10)
result1
# Batch size 20
result1 = analyze(0.00001, 6, 20)
result1<jupyter_output>Theta for degree 2 done in: 0.001145999999998537

Theta for degree 3 done in: 0.0009610000000037644

Theta for degree 5 done in: 0.0008180000000024279<jupyter_text>With the hyperparameter of `number of iteration` is 10 and `learning rate` is 0.00001 and we split the data by 10, so we get 12 partition. We get smaller RMSE compared to the closed form for our training data.The error we get varies depends on our hyperparameters we mention above. In case of polinomial of 5, we get better result if we increase the batch size and adjust the number of iteration accordingly.With batch size of 10, we get the time decreasing from degree 2, 3 to 5 but for batch size of 20, the time is similar, except for degree 2. Therefore, we can conclude there is no significant influence of the batch time to the performance in this case.We have known bug that if we set the learning rate or the number of iteration too high, we will overflow exception (known bugs are marked as ` FIXME` in this report). > Apply Ridge regression to the dataset to compute θ for a given λ. Implement two cases: 2using a closed-form solution and using a stochastic gradient descent method with mini-batches of size m. Use K-fold cross validation on the training dataset to obtain the best regularization λ and apply the optimal θ to compute the regression error on test samples. Report the optimal λ, θ, test and training set errors for K ∈ {2,10,N}, where N is the number of samples. In all cases try n ∈ {2, 3, 5}. How does the test error change as a function of λ and n?<jupyter_code>class RidgeRegression(object):
    def __init__(self, lmbda=0.1):
        self.lmbda = lmbda
    
    def fit(self, X, y):
        C = X.T.dot(X) + self.lmbda * np.eye(X.shape[1])
        self.w = np.linalg.inv(C).dot(X.T.dot(y))
    
    def set_params(self, lmbda=0.1):
        self.lmbda = lmbda
        return self
        
    def get_weight(self):
        return self.w
    
    # def get_error(self):
    #     y_hat = np.dot(self.X, self.w)
    #     return np.sqrt(((self.Y - y_hat) ** 2).mean())
class Minimisation(object):
    def __init__(self, X, y, model):
        self.model = model
        # Prepare the data
        self.X = np.array(X)
        self.y = np.array(y)
        
    def closed_form(self):
        self.model.fit(self.X, self.y)
        return self.model.get_weight()
    
    def mini_batch_stochastic(self, learning_rate, num_iter, batch_size):
        # Randomize the data and split into samples
        X = self.X
        y = self.y
        
        df = shuffle_data(X, y, batch_size)
    
        # Get the dimension
        x_dim = X.shape[1]
    
        # Set the initial tetha and other initial variables
        theta = np.zeros((x_dim, 1))
    
        start_i = 0
    
        # Create the loop
        for i in range(start_i + 1, num_iter + 1):
            # For each randomized-mini-batch partition
            for j in range(0, len(df)):
                X = df[j][0]
                y = df[j][1]
                
                self.model.fit(self.X, self.y)
                theta = theta - learning_rate * self.model.get_weight() 
        
        return theta
        
    def shuffle_data(self, batch_size):
        X = self.X
        y = self.y
    
        random_indices = np.random.choice(len(X), len(y), replace=False)
    
        X_shuffled = X[random_indices,:]
        y_shuffled = y[random_indices]
        mini_batches = [(X_shuffled[i:i+batch_size,:], y_shuffled[i:i+batch_size]) for
                   i in range(0, len(y), batch_size)]
        return mini_batches

# Load train data
X_trn2 = df2['X_trn']
Y_trn2 = df2['Y_trn']


# Load test data
X_test2 = df2['X_tst']
Y_test2 = df2['Y_tst']
def k_fold(X, y, lmbda, fold):
    # Split the data into fold part
    # Prepare the data
    X = np.array(X)
    y = np.array(y)
    ridge = RidgeRegression()
    
    # Randomize the data and split into samples
    batch_size = int(len(X) / fold)
    df = shuffle_data(X, y, batch_size)
    
    test_result = []
    
    for i in range(0, len(df)):
        # Temporary data for k-fold
        temp = df.copy()
        
        # To store partitions for training
        train = []
        
        # Save and pop an element for test
        test = temp.pop(i)
        
        # Concat the remaining array
        X_train = np.empty(shape=[X.shape[0], X.shape[1]])
        Y_train = np.empty(shape=[y.shape[0], y.shape[1]])
        
        for j in range(0, len(df)):
            _X = df[j][0]
            _y = df[j][1]
            X_train = np.concatenate((X_train, _X))
            Y_train = np.concatenate((Y_train, _y))
        
        
        ridge.set_params(lmbda)
        # test[0] is X in test data
        # test[1] is Y in test data
        
                            
        ridge.fit(X_train, Y_train)
        #print(X_train, Y_train)
        w = ridge.get_weight()
        
        if ~(math.isnan(w[0])):
            y_hat = np.dot(test[0], w)
            test_result.append(np.sqrt(((test[1] - y_hat) ** 2).mean()))
    
    return np.nanmean(test_result)
def k_fold_analysis(degree, fold):
    X = np.array(X_trn2)
    k_fold_result = []
    
    X_trn_new = np.ndarray((len(X_trn2), degree + 1))
    i = 0
    for x in X:
        X_trn_new[i] = create_polinomial_array(x, degree)
        i += 1
        
    
    for j in range(-10, 2):
        k_fold_result.append(k_fold(X_trn_new, Y_trn2, 10 ** j, fold))
    
    return k_fold_result

degrees = [2, 3, 5]
folds = [2, 10, 20]

for degree in degrees:
    for fold in folds:
        temp = k_fold_analysis(degree, fold)
        print("Polynomial degree: {}, K-Fold: {}\n".format(degree, fold))
        print("Minimum error at index: {} value {} with lambda: 10^{}\n".format(
            np.argmin(temp), np.min(temp), (np.argmin(temp) - 10)))<jupyter_output>/usr/local/lib/python3.6/site-packages/ipykernel_launcher.py:48: RuntimeWarning: Mean of empty slice<jupyter_text>From the above result, we see that the best lambda is depending upon the number of parameters in X--polynomial with degree 5 has different lambda that result the smallest error compared to other polynomial degree. With the values above, let's try to compute the weight value using closed form and stochastic. Note: The value can change since we shuffle the data.<jupyter_code>error_closed_form_test = []
error_grad_desc_test = []

X = np.array(X_trn2)
X_tst = np.array(X_test2)
# For degree 2
degree = 2
X_trn_new = np.ndarray((len(X_trn2), degree + 1))
X_tst_new = np.ndarray((len(X_test2), degree + 1))
        
i = 0
for x in X:
    X_trn_new[i] = create_polinomial_array(x, degree)
    i += 1
    
i = 0
for x in X_tst:
    X_tst_new[i] = create_polinomial_array(x, degree)
    i += 1

ridge = RidgeRegression()
ridge.set_params(10**-9)
m = Minimisation(X_trn_new, Y_trn2, ridge)

w1 = m.closed_form()
y_hat1 = np.dot(X_tst_new, w1)
e1 = np.sqrt(((Y_test2 - y_hat1) ** 2).mean())
error_closed_form_test.append(e1)

w2 = m.mini_batch_stochastic(0.00001, 10, 10)
y_hat2 = np.dot(X_tst_new, w2)
e2 = np.sqrt(((Y_test2 -  y_hat2) ** 2).mean())
error_grad_desc_test.append(e2)
d = {
        'degree': [2],
        'close_form_error_test': error_closed_form_test,
        'grad_desc_error_test': error_grad_desc_test
    }

# Print the error in close form and gradient descent using the theta we get from the close form
d<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10455013989103225, 'fraction_numerical': 0.021204535414519218, 'mean_word_length': 2.5660278288264635, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '222868', 'license': 'MIT', 'n_tokens_mistral': 4931, 'n_tokens_neox': 4656, 'n_words': 1346, 'path': 'LinearRegression.ipynb', 'repo_name': 'prdx/MLFS', 'size': '13610'}