{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"#
regularization visualized & explained⚙️
\n\n
\n\n**What you can expect from this notebook:** In this notebook **lasso(L1)**-, **ridge(L2)**- and **elasticnet**-regularization are explained based on some animated visualisations. After that section, a basic use case is covered by performing **regularized linear regression** on the house price dataset.\n\n
👉If you're just interested in the complete, with comments documented implementation of a linear model with regularization, feel free to click on show hidden code:
","metadata":{}},{"cell_type":"code","source":"class MultivariateLinearModel:\n def __init__(self):\n self.feature_func_ = False\n self.weights = None\n \n def fit(self, X, y, optimizer, epochs=100, learning_rate=0.01, batch_size=32, penalty='none', alpha=1, elastic_net_ratio=0.5, verbose=True):\n # initialize weight vecotor\n self.weights = np.zeros(len(self.feature_func)*len(X[0])) \n # initialize optimizer\n opt = optimizer(learning_rate)\n # iterate over epochs\n for e_i in range(epochs):\n # split in mini batches\n num_batches = int(np.ceil(len(X) / batch_size))\n batches = [np.array_split(X, num_batches), np.array_split(y, num_batches)]\n # iterate over batches\n for b_i in range(num_batches):\n # calculate gradient estimate\n y_pred = self.predict(batches[0][b_i])\n y_true = batches[1][b_i]\n features = np.concatenate([np.column_stack(func(batches[0][b_i])) for func in self.feature_func])\n error = y_true - y_pred\n gradient_point_estimates = -2 * error * features\n if penalty == 'none':\n penalty_gradient = 0\n elif penalty == 'l1':\n penalty_gradient = alpha * np.sign(self.weights)\n elif penalty == 'l2':\n penalty_gradient = alpha * self.weights\n elif penalty == 'elastic_net':\n penalty_gradient = elastic_net_ratio * alpha * np.sign(self.weights) + elastic_net_ratio * alpha * self.weights\n else:\n raise ValueError('please choose a penalty among none, l1, l2 and elastic_net')\n gradient_mean_estimate = (gradient_point_estimates + penalty_gradient).mean(axis=1)\n # weight update using optimizer\n self.weights = opt.update_weights(self.weights, gradient_mean_estimate)\n \n if verbose:\n mse = ((self.predict(X) - y) ** 2).mean()\n print(f'epoch {e_i + 1}, mse: {mse}')\n\n def predict(self, X):\n features = np.concatenate([np.column_stack(func(X)) for func in self.feature_func])\n return features.T @ self.weights\n \n def score(self, X, y):\n r_2 = 1 - (np.square(y - self.predict(X)).sum() / np.square(y - y.mean()).sum())\n mse = np.square(y - self.predict(X)).mean()\n print(f'model scored R2 of {r_2} and mse {mse}')\n return r_2, mse\n \n @property\n def feature_func(self):\n if type(self.feature_func_) != list:\n raise ValueError('please define the feature function before moving on')\n return self.feature_func_\n \n @feature_func.setter\n def feature_func(self, func_array): \n self.feature_func_ = [np.vectorize(func) for func in func_array]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-12T10:11:39.662092Z","iopub.execute_input":"2022-08-12T10:11:39.663177Z","iopub.status.idle":"2022-08-12T10:11:39.711165Z","shell.execute_reply.started":"2022-08-12T10:11:39.663066Z","shell.execute_reply":"2022-08-12T10:11:39.710208Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"****\n\n# 1 | Intuition\n\nAt its core, regularization is a way to **control the bias/variance ratio** of a model by **penalizing complexity** by a set factor.\n
\n
img src: staff.fnwi.uva.nl
\n\n
\n\n**How regularization can penalize model complexity** can be demonstrated well with polynomial regression:\n* the model gets more complex and **more likely to overfit the more x^n terms are added**\n * -> the **model variance increases** and it's likely to fit the noise instead of just the signal\n* how much (or even if) an x^n term contributes to the model **complexity is determined by its coefficient**\n * -> a polynomial of 20th order is equal to a polynomial of second order if every x^n term after x^2 has a coefficient of 0\n* **So by penalizing high (and many) coefficients, model complexity can be regularized.**\n\nThis concept can be generalized further: everything in a model weighted by some parameter (p.e. other variables) has the potential of increasing that model's variance and can be regularized by controlling the magnitude of its coefficient. \n\n
\n\n**How regularization operates** generally:\n1. Compute a **penalty that scales with more/larger parameters**\n2. **Scale that penalty** by some factor to determine 'how strongly' it influences the model parameters\n3. **Add this scaled penalty to an existing loss** function\n * this way reducing the loss also **requires keeping the parameter magnitudes low**\n \n
\n\n**Example formula** (MSE with penalty):\n\n$\\Large MSE_{regularized}=\\frac{1}{n}\\sum\\limits_{i=1}^N(y_i-\\hat{y}_i) + \\alpha \\cdot penalty$","metadata":{}},{"cell_type":"markdown","source":"****\n\n# 2 | Types of penalties\n\n**Formula overview:**\n\n| penalty: |
L1 (Lasso)
|
L2 (Ridge)
|
ElasticNet
|\n|---|---|---|---|\n| function: | $\\Large \\alpha\\lVert \\theta\\rVert_1 = \\alpha \\sum\\limits_iabs(\\theta_i)$
| $\\Large \\alpha\\lVert \\theta\\rVert_2^2 = \\alpha \\theta^T\\theta = \\alpha \\sum\\limits_i\\theta_i^2$
| $\\Large \\alpha(r\\lVert \\theta\\rVert_1 + (1-r)\\lVert \\theta\\rVert_2^2)$ |\n| gradient: | $\\Large \\alpha \\> sign(\\theta)$
| $\\Large 2\\alpha\\theta$
| $\\Large \\alpha \\> r \\> sign(\\theta) + 2\\alpha (1-r)\\theta$ | \n\n\n**Notes:**\n* **alpha is a scaling parameter** that controls the **Impact** of the penalty when added to a loss function\n* In **L1 regularization**, the **absolute value** of every parameter is summed up \n* In **L2 regularization**, the **squared value** of every parameter is summed up\n* **ElasticNet regularization** combines both **L1 and L2 regularization** with an **ratio r**\n\n
\n\n*To see how these formulas look visually and learn what differentiates them from each other, scroll down to the next section*\n","metadata":{}},{"cell_type":"markdown","source":"****\n\n# 3 | Animations & derived properties/details of different penalties\n","metadata":{}},{"cell_type":"markdown","source":"**Animation setup (necessary details will be covered, so no need to inspect the code):**","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom matplotlib import rc\nfrom matplotlib.animation import FuncAnimation\nfrom sklearn.experimental import enable_iterative_imputer\nfrom sklearn.impute import IterativeImputer\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import train_test_split\nfrom IPython.display import HTML\npd.options.mode.chained_assignment = None \n\ndf = pd.read_csv('../input/house-prices-advanced-regression-techniques/train.csv')\n\n############################### preprocess data: #################################\ndef start_pipeline(df):\n return df.copy()\n\ndef fix_missing_values(df):\n df.drop(columns=df.columns[df.isna().sum() > df.count()/4], inplace=True)\n \n # numerical\n num = df.iloc[:, (df.dtypes != object).values]\n imputer = IterativeImputer(max_iter=50, random_state=42)\n num_values = imputer.fit_transform(num)\n num[num.columns] = num_values\n \n # categorical\n cat = df.iloc[:, (df.dtypes == object).values]\n cat.fillna(df.mode().iloc[0, :], inplace=True)\n \n return pd.concat([num, cat], axis=1)\n\ndef encode_categorical(df):\n for c in df.columns[(df.dtypes == object).values]:\n df[c] = df[c].astype('category').cat.codes\n return df\n\ndef scale(df):\n scaler = StandardScaler()\n df[df.columns] = scaler.fit_transform(df)\n return df\n\npreprocessed = df.pipe(start_pipeline).pipe(fix_missing_values).pipe(encode_categorical).pipe(scale)\n\nx_ = preprocessed['OverallQual'].sample(100).to_numpy()\nx = np.repeat(np.repeat(np.expand_dims(np.expand_dims(x_, 0), 0), 200, axis=0), 200, axis=1)\ny_ = preprocessed['SalePrice'].sample(100).to_numpy() \ny = np.repeat(np.repeat(np.expand_dims(np.expand_dims(y_ + x_ * 3 + 0.7 * x_ ** 2, 0), 0), 200, axis=0), 200, axis=1)\ny_ = y_ + 3.5 * x_\n\nrc('animation', html='jshtml') \n\n############################### 1D, all penalties animation: #################################\ndef make_1d_animation():\n fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(32, 15))\n\n ax1.set_xlim([-1, 5])\n ax2.set_xlim([-1, 5])\n ax1.set_ylim([0, 20])\n ax2.set_ylim([-20, 20])\n\n L1, = ax1.plot([], [], color='purple', label='none')\n L2 = ax1.scatter([], [], color='purple')\n\n l3, = ax1.plot([], [], color='red', label='L1')\n l4, = ax1.plot([], [], color='red', alpha=0.6, linestyle=':')\n l5 = ax1.scatter([], [], color='red')\n\n l6, = ax1.plot([], [], color='green', label='L2')\n l7 = ax1.scatter([], [], color='green')\n l8, = ax1.plot([], [], color='green', alpha=0.6, linestyle=':')\n\n ax1.axvline(0, linestyle='--', color='black')\n ax1.axhline(0, linestyle='--', color='black')\n ax1.legend(title='penalty:', loc='upper right')\n ax1.set_title('loss function')\n\n l9, = ax2.plot([], [], color='purple', label='none')\n l10 = ax2.scatter([], [], color='purple')\n\n l11, = ax2.plot([], [], color='red', label='L1')\n l12 = ax2.scatter([], [], color='red')\n l13, = ax2.plot([], [], color='red', alpha=0.6, linestyle=':')\n\n l14, = ax2.plot([], [], color='green', label='L2')\n l15 = ax2.scatter([], [], color='green')\n l16, = ax2.plot([], [], color='green', alpha=0.6, linestyle=':')\n\n ax2.axvline(0, linestyle='--', color='black')\n ax2.axhline(0, linestyle='--', color='black')\n ax2.legend(title='penalty:', loc='upper right')\n ax2.set_title('loss function gradient')\n fig.tight_layout()\n\n def init():\n L1, L2, l3, l4, l5, l6, l7, l8, l9, l10, l11, l12, l13, l14, l15, l16 = visualise_alpha_impact(alpha=0)\n return L1, L2, l3, l4, l5, l6, l7, l8, l9, l10, l11, l12, l13, l14, l15, l16\n\n def visualise_alpha_impact(alpha):\n ############################### computation: #################################\n ß1 = np.linspace(-1, 5, 300)\n\n loss = np.square(y_.reshape(-1, 1) - ß1 * x_.reshape(-1, 1)).mean(axis=0)\n gradient = (-2 * x_.reshape(-1, 1) * (y_.reshape(-1, 1) - ß1 * x_.reshape(-1, 1))).mean(axis=0)\n\n loss_l1 = np.square(y_.reshape(-1, 1) - ß1 * x_.reshape(-1, 1)).mean(axis=0) + alpha * np.abs(ß1)\n gradient_l1 = (-2 * x_.reshape(-1, 1) * (y_.reshape(-1, 1) - ß1 * x_.reshape(-1, 1))).mean(axis=0) + alpha * np.sign(ß1)\n l1 = alpha * np.abs(ß1)\n l1_dev = alpha * np.sign(ß1)\n\n loss_l2 = np.square(y_.reshape(-1, 1) - ß1 * x_.reshape(-1, 1)).mean(axis=0) + alpha / 2 * np.square(ß1)\n gradient_l2 = (-2 * x_.reshape(-1, 1) * (y_.reshape(-1, 1) - ß1 * x_.reshape(-1, 1))).mean(axis=0) + alpha * ß1\n l2 = alpha / 2 * np.square(ß1)\n l2_dev = alpha * ß1\n\n ############################### visualisation: #################################\n L1.set_data(ß1, loss) \n L2.set_offsets([[ß1[np.argmin(loss)]] + [loss[np.argmin(loss)]]]) \n\n l3.set_data(ß1, loss_l1) \n l4.set_data(ß1, l1) \n l5.set_offsets([[ß1[np.argmin(loss_l1)]] + [loss_l1[np.argmin(loss_l1)]]]) \n\n l6.set_data(ß1, loss_l2) \n l7.set_offsets([[ß1[np.argmin(loss_l2)]] + [loss_l2[np.argmin(loss_l2)]]]) \n l8.set_data(ß1, l2) \n\n l9.set_data(ß1, gradient)\n l10.set_offsets([[ß1[np.argmin(np.abs(gradient))]] + [0]])\n\n l11.set_data(ß1, gradient_l1) \n l12.set_offsets([[ß1[np.argmin(np.abs(gradient_l1))]] + [0]]) \n l13.set_data(ß1, l1_dev) \n\n l14.set_data(ß1, gradient_l2) \n l15.set_offsets([[ß1[np.argmin(np.abs(gradient_l2))]] + [0]])\n l16.set_data(ß1, l2_dev)\n\n return L1, L2, l3, l4, l5, l6, l7, l8, l9, l10, l11, l12, l13, l14, l15, l16\n\n ani = FuncAnimation(fig, visualise_alpha_impact, frames=np.linspace(0.1, 10, 99), blit=True, init_func=init)\n return ani\n\n############################### 2D animation: #################################\ndef make_2d_animation(penalty_):\n model = np.vectorize(lambda ß0, ß1, ß2, x: ß0 + ß1 * x + ß2 * x ** 2) \n\n ß1 = np.linspace(-5, 5, 200)\n ß2 = np.linspace(-5, 5, 200)\n ß1ß1, ß2ß2 = np.meshgrid(ß1, ß2)\n\n penalty = penalty_(ß1ß1, ß2ß2) \n loss = np.mean((y - model(0, np.expand_dims(ß1ß1, -1), np.expand_dims(ß2ß2, -1), x)) ** 2, axis=2)\n\n fig2, (ax1, ax2) = plt.subplots(1, 2, figsize=(22, 12))\n\n ax1.set_xlim([-5, 5])\n ax2.set_xlim([-5, 5])\n ax1.set_ylim([-5, 5])\n ax2.set_ylim([-5, 5])\n\n ax1.contourf(ß1, ß2, loss, cmap=plt.get_cmap('bone_r'))\n ax1.contour(ß1, ß2, penalty, cmap=plt.get_cmap('winter_r'))\n ax1.axvline(0, linestyle='--', color='black')\n ax1.axhline(0, linestyle='--', color='black')\n ax1.set_title('loss(filled) vs. penalty(contour)')\n \n cmap = plt.get_cmap('bone_r')\n a1 = ax2.contourf([0,0], [0,0], [[0,0], [0,0]], cmap=cmap)\n ax2.axvline(0, linestyle='--', color='black')\n ax2.axhline(0, linestyle='--', color='black')\n ax2.set_title('loss + penalty')\n\n fig2.tight_layout()\n\n def init2():\n a1 = visualise_alpha_impact_L1(alpha=0)\n return a1\n\n def visualise_alpha_impact_L1(alpha, loss=loss, penalty=penalty, a1=a1, cmap=cmap):\n ############################### computation: #################################\n loss_w_penalty = loss + alpha * penalty \n\n ############################### visualisation: ###############################\n for c in a1.collections:\n try:\n c.remove() \n except ValueError:\n pass\n a1 = ax2.contourf(ß1, ß2, loss_w_penalty, cmap=cmap)\n\n return a1\n\n ani = FuncAnimation(fig2, visualise_alpha_impact_L1, frames=np.linspace(0.1, 10, 99), repeat=False, init_func=init2)\n return ani","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-12T14:19:50.142290Z","iopub.execute_input":"2022-08-12T14:19:50.143746Z","iopub.status.idle":"2022-08-12T14:19:51.903825Z","shell.execute_reply.started":"2022-08-12T14:19:50.143698Z","shell.execute_reply":"2022-08-12T14:19:51.902726Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### L1 and L2 regularization compared (1D example):","metadata":{}},{"cell_type":"code","source":"%matplotlib agg\n\nd1_ani = make_1d_animation()\nd1_ani.save('regularization_impact_alpha_1D.gif', writer='imagemagick')\n\nHTML(d1_ani.to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2022-08-12T10:17:34.579251Z","iopub.execute_input":"2022-08-12T10:17:34.579699Z","iopub.status.idle":"2022-08-12T10:18:00.365599Z","shell.execute_reply.started":"2022-08-12T10:17:34.579661Z","shell.execute_reply":"2022-08-12T10:18:00.363622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Notes:**\n* every frame corresponds to an increase of **alpha by 0.1**\n* The **translucent graphs** in the background correspond to the loss/derivative of the **penalty terms only**\n* The **graphs in the foreground** correspond to the loss/derivative of the **MSE loss regularized by the penalty terms**\n* The **points** marked on the graphs are the **optimal parameters** obtained from that loss function\n\n**Things to take away from the animation:**\n* L1 **can set parameters to 0** due to the proportional to alpha increasing y-range at x=0 of its derivative -> **promotes sparsity**\n * **L1 can be thought of** as **affecting the intercept** of the loss functions gradient (while keeping its slope the same)\n* L2 **shrinks parameters but never sets them exactly to 0** due to it, increasing alpha, just increases the slope of the derivative (compared to the loss function derivative without penalty) but does not change the y-intercept\n * **L2 can be thought of** as **affecting the slope** of the loss functions gradient (while keeping its intercept the same)\n* increasing alpha of **L1** regularization has a negative linear proportional effect on the magnitudes of the parameters -> increasing alpha reduces the parameter magnitudes linearly\n* increasing alpha of **L2** regularization has a negative quadratic effect on the magnitudes of the parameters -> increasing alpha reduces the parameter magnitudes less the higher alpha gets","metadata":{}},{"cell_type":"markdown","source":"### L1 regularization (2D example):","metadata":{}},{"cell_type":"code","source":"%matplotlib agg\n\nl1_penalty_ = np.vectorize(lambda ß1, ß2: np.abs(ß1) + np.abs(ß2))\n\nd2_ani_l1 = make_2d_animation(l1_penalty_)\nd2_ani_l1.save('L1_regularization_impact_alpha_2D.gif', writer='imagemagick')\n\nHTML(d2_ani_l1.to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2022-08-12T12:50:48.221288Z","iopub.execute_input":"2022-08-12T12:50:48.221761Z","iopub.status.idle":"2022-08-12T12:51:06.929278Z","shell.execute_reply.started":"2022-08-12T12:50:48.221723Z","shell.execute_reply":"2022-08-12T12:51:06.927888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Notes:**\n* darker colours correspond to higher values\n\n**Things to take away from the animation:**\n* with multiple parameters, L1 regularization tends to eliminate some features while keeping the most important\n * -> can be used as a feature selection technique","metadata":{}},{"cell_type":"markdown","source":"### L2 regularization (2D example):","metadata":{}},{"cell_type":"code","source":"%matplotlib agg\n\nl2_penalty_ = np.vectorize(lambda ß1, ß2: np.square(ß1) + np.square(ß2))\n\nd2_ani_l2 = make_2d_animation(l2_penalty_)\nd2_ani_l2.save('L2_regularization_impact_alpha_2D.gif', writer='imagemagick')\nHTML(d2_ani_l2.to_jshtml())","metadata":{"execution":{"iopub.status.busy":"2022-08-12T12:49:47.620333Z","iopub.execute_input":"2022-08-12T12:49:47.620790Z","iopub.status.idle":"2022-08-12T12:50:06.228303Z","shell.execute_reply.started":"2022-08-12T12:49:47.620744Z","shell.execute_reply":"2022-08-12T12:50:06.227284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Things to take away from the animation:**\n* in contrast to L1 regularization, L2 regularization doesn't operate as a feature selection method but rather as a feature 'weighting/scaling' method","metadata":{}},{"cell_type":"markdown","source":"****\n\n# 4 | Python implementation & Regression on house price dataset\n\n### libraries used:","metadata":{}},{"cell_type":"code","source":"import pandas as pd # for data handeling\nimport numpy as np # for linear algebra\n# preprocessing:\nfrom sklearn.experimental import enable_iterative_imputer\nfrom sklearn.impute import IterativeImputer\nfrom sklearn.preprocessing import StandardScaler\nfrom copy import copy # make deep copies of objects -> needed for cv implementation","metadata":{"execution":{"iopub.status.busy":"2022-08-12T14:24:16.335361Z","iopub.execute_input":"2022-08-12T14:24:16.335856Z","iopub.status.idle":"2022-08-12T14:24:16.342115Z","shell.execute_reply.started":"2022-08-12T14:24:16.335818Z","shell.execute_reply":"2022-08-12T14:24:16.341249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### mathematical recap:\n\n**loss functions:**\n\n$\\large SSE_{L1} = \\lVert y-\\hat{y}\\rVert_2^2 + \\alpha\\lVert \\theta\\rVert_1= \\lVert y-\\phi_X^T\\theta\\rVert_2^2 + \\alpha\\lVert \\theta\\rVert_1= (y-\\phi_X^T\\theta)^T(y-\\phi_X^T\\theta)=y^Ty-y^T\\phi_X^T\\theta-\\phi_X\\theta^Ty+\\phi_X\\theta^T\\phi_X^T\\theta + \\alpha \\vec{1}^T \\left\\{ \\begin{array}{ c l }\\theta_i & \\quad \\textrm{if } \\theta_i >= 0 \\\\ -\\theta_i & \\quad \\textrm{if } \\theta_i < 0 \\end{array}\\right.$\n\n$\\large SSE_{L2} = \\lVert y-\\hat{y}\\rVert_2^2 + \\alpha\\lVert \\theta\\rVert_2^2 = \\lVert y-\\phi_X^T\\theta\\rVert_2^2 = (y-\\phi_X^T\\theta)^T(y-\\phi_X^T\\theta) + \\alpha \\theta^T \\theta=y^Ty-y^T\\phi_X^T\\theta-\\phi_X\\theta^Ty+\\phi_X\\theta^T\\phi_X^T\\theta + \\alpha \\theta^T \\theta$\n\n**derivatives:**\n\n$\\large \\frac{\\partial \\>\\>y^Ty-y^T\\phi_X^T\\theta-\\phi_X\\theta^Ty+\\phi_X\\theta^T\\phi_X^T\\theta + \\alpha \\lVert \\theta\\rVert_1}{\\partial \\>\\>\\theta} = -2\\phi_Xy + 2\\phi_X\\phi_X^T\\theta + \\alpha \\> sign(\\theta)$
\n\n$\\large \\frac{\\partial \\>\\>y^Ty-y^T\\phi_X^T\\theta-\\phi_X\\theta^Ty+\\phi_X\\theta^T\\phi_X^T\\theta + \\alpha \\theta^T \\theta}{\\partial \\>\\>\\theta} = -2\\phi_Xy + 2\\phi_X\\phi_X^T\\theta + 2\\alpha \\theta = -\\phi_Xy + (\\phi_X\\phi_X^T + \\alpha I) \\theta$
\n\n### implementation:\nLet's implement a linear model with all covered types of regularization. For this implementation, I'll build on the linear model derived in [this notebook](https://www.kaggle.com/code/vincentbrunner/ml-from-scratch-custom-linear-models). ","metadata":{}},{"cell_type":"markdown","source":"As optimizer we make usage of the Adam algorithm. Adam is a gradient descent based optimizer I explained here: [ml from scratch: neural network and GD-optimizers](https://www.kaggle.com/code/vincentbrunner/ml-from-scratch-neural-network-and-gd-optimizers). It's a bit overkill for this optimization problem, but I had fun playing around with the implementation a bit.","metadata":{}},{"cell_type":"code","source":"class Adam:\n def __init__(self, learning_rate, ß1=0.9, ß2=0.99):\n self.lr = learning_rate\n self.ß1 = ß1\n self.ß2 = ß2\n self.v = 0\n self.s = 0\n self.t = 1\n self.e = 1e-10\n \n def update_weights(self, weights, gradient):\n # exp. moving average term for gradient -> 'momentum buffer'\n self.v = self.ß1 * self.v + (1 - self.ß1) * gradient\n # exp. moving average term for squared gradient -> addaptive learning rate\n self.s = self.ß2 * self.s + (1 - self.ß2) * gradient ** 2\n # bias correction for both terms\n v_corrected = self.v / (1 - self.ß1 ** self.t)\n s_corrected = self.s / (1 - self.ß2 ** self.t)\n self.t += 1\n # final weight update\n return weights - self.lr * (v_corrected / (np.sqrt(s_corrected) + self.e))","metadata":{"execution":{"iopub.status.busy":"2022-08-12T13:47:09.523089Z","iopub.execute_input":"2022-08-12T13:47:09.524088Z","iopub.status.idle":"2022-08-12T13:47:09.533599Z","shell.execute_reply.started":"2022-08-12T13:47:09.524046Z","shell.execute_reply":"2022-08-12T13:47:09.532271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class MultivariateLinearModel:\n def __init__(self, penalty='none', alpha=1, elastic_net_ratio=0.5):\n self.feature_func_ = False\n self.weights = None\n self.penalty = penalty\n self.alpha = alpha\n self.elastic_net_ratio = elastic_net_ratio\n \n def fit(self, X, y, optimizer=Adam, epochs=500, learning_rate=1, batch_size=16, verbose=False):\n # initialize weight vecotor -> one coefficient per term per variable\n self.weights = np.zeros(len(self.feature_func)*len(X[0])) \n # initialize optimizer\n opt = optimizer(learning_rate)\n # iterate over epochs\n for e_i in range(epochs):\n # split in mini batches\n num_batches = int(np.ceil(len(X) / batch_size))\n batches = [np.array_split(X, num_batches), np.array_split(y, num_batches)]\n # iterate over batches\n for b_i in range(num_batches):\n # calculate gradient estimate of loss\n y_pred = self.predict(batches[0][b_i])\n y_true = batches[1][b_i]\n features = np.concatenate([np.column_stack(func(batches[0][b_i])) for func in self.feature_func])\n error = y_true - y_pred\n gradient_point_estimates = -2 * error * features\n # calculate gradient estimate of penalty\n if self.penalty == 'none':\n penalty_gradient = 0.0\n elif self.penalty == 'l1':\n penalty_gradient = self.alpha * np.sign(self.weights)\n elif self.penalty == 'l2':\n penalty_gradient = self.alpha * self.weights\n elif self.penalty == 'elastic_net':\n penalty_gradient = self.elastic_net_ratio * self.alpha * np.sign(self.weights) + self.elastic_net_ratio * self.alpha * self.weights\n else:\n raise ValueError('please choose a penalty among none, l1, l2 and elastic_net')\n # average point estimates over mini batch\n gradient_mean_estimate = (gradient_point_estimates).mean(axis=1) + penalty_gradient\n # weight update using optimizer\n self.weights = opt.update_weights(self.weights, gradient_mean_estimate)\n \n if verbose:\n if (e_i + 1) % 10 == 0:\n mse = ((self.predict(X) - y) ** 2).mean()\n print(f'epoch {e_i + 1}, mse: {mse}')\n\n def predict(self, X):\n # passing X through the feature function to obtain the calculated terms \n features = np.concatenate([np.column_stack(func(X)) for func in self.feature_func])\n # linear combination of terms with their coefficients\n return features.T @ self.weights\n \n def score(self, X, y, verbose=False):\n # coefficient of determination\n r_2 = 1 - (np.square(y - self.predict(X)).sum() / np.square(y - y.mean()).sum())\n # mean squared error\n mse = np.square(y - self.predict(X)).mean()\n if verbose == True:\n print(f'model scored R2 of {r_2} and mse {mse}')\n return r_2\n \n # using a property to deal with the assignment of the feature function\n @property\n def feature_func(self):\n # checking if the feature function is allready set\n if type(self.feature_func_) != list:\n raise ValueError('please define the feature function before moving on')\n return self.feature_func_\n \n @feature_func.setter\n def feature_func(self, func_array): \n # np.vectorize allows to pass a numpy array trough the function\n self.feature_func_ = [np.vectorize(func) for func in func_array]","metadata":{"execution":{"iopub.status.busy":"2022-08-12T14:13:27.581011Z","iopub.execute_input":"2022-08-12T14:13:27.581425Z","iopub.status.idle":"2022-08-12T14:13:27.602264Z","shell.execute_reply.started":"2022-08-12T14:13:27.581386Z","shell.execute_reply":"2022-08-12T14:13:27.600772Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### basic data preprocessing pipeline:\n1. split into train/test set\n2. drop columns with too many missing values\n3. imputing the rest \n4. encoding categorical data\n5. **scaling data**: this is especially important with regularization cause features on a larger scale would heavily influence the penalty","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv('../input/house-prices-advanced-regression-techniques/train.csv')\n\n############################### preprocess data: #################################\ndef start_pipeline(df, test_split=0.2):\n new_df = df.copy()\n df['is_train'] = 1\n df['is_train'][df.index.isin(df.sample(frac=test_split).index)] = 0\n return df\n\ndef fix_missing_values(df):\n df.drop(columns=df.columns[df.isna().sum() > df.count()/4], inplace=True)\n \n # numerical\n num = df.iloc[:, (df.dtypes != object).values]\n imputer = IterativeImputer(max_iter=50, random_state=42)\n imputer.fit(num)\n num_values = imputer.transform(num)\n num[num.columns] = num_values\n \n # categorical\n cat = df.iloc[:, (df.dtypes == object).values]\n cat.fillna(df.mode().iloc[0, :], inplace=True)\n \n return pd.concat([num, cat], axis=1)\n\ndef encode_categorical(df):\n for c in df.columns[(df.dtypes == object).values]:\n df[c] = df[c].astype('category').cat.codes\n return df\n\ndef scale(df):\n scaler = StandardScaler()\n scaler.fit(df[df['is_train'] == 1].iloc[:, ~df.columns.isin(['is_train', 'SalePrice'])])\n df.iloc[:, ~df.columns.isin(['is_train', 'SalePrice'])] = scaler.transform(df.iloc[:, ~df.columns.isin(['is_train', 'SalePrice'])])\n return df\n\npreprocessed = df.pipe(start_pipeline).pipe(fix_missing_values).pipe(encode_categorical).pipe(scale)\npreprocessed.head()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-12T14:27:56.021443Z","iopub.execute_input":"2022-08-12T14:27:56.021902Z","iopub.status.idle":"2022-08-12T14:27:57.917096Z","shell.execute_reply.started":"2022-08-12T14:27:56.021868Z","shell.execute_reply":"2022-08-12T14:27:57.916245Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# splitting into features and labels (X & y)\nX_train = preprocessed[preprocessed['is_train'] == 1].iloc[:, preprocessed.columns != 'SalePrice'].to_numpy()\ny_train = preprocessed[preprocessed['is_train'] == 1].iloc[:, preprocessed.columns == 'SalePrice'].to_numpy().reshape(-1)\nX_test = preprocessed[preprocessed['is_train'] == 0].iloc[:, preprocessed.columns != 'SalePrice'].to_numpy()\ny_test = preprocessed[preprocessed['is_train'] == 0].iloc[:, preprocessed.columns == 'SalePrice'].to_numpy().reshape(-1)","metadata":{"execution":{"iopub.status.busy":"2022-08-12T14:28:21.223684Z","iopub.execute_input":"2022-08-12T14:28:21.224383Z","iopub.status.idle":"2022-08-12T14:28:21.238539Z","shell.execute_reply.started":"2022-08-12T14:28:21.224339Z","shell.execute_reply":"2022-08-12T14:28:21.237453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### House price prediction using linear regression with L1 regularization\n* choosing best alpha utilising crossvalidation (covered [here](https://www.kaggle.com/code/vincentbrunner/ml-from-scratch-nested-cross-validation))\n\n**the crossvalidation implementation from that notebook:**","metadata":{}},{"cell_type":"code","source":"# turn into class to work with sklearn model pipeline and allow for easy nesting of cross-validation loops\nclass KfoldCrossval:\n def __init__(self, k, models, name=\"k_fold_crossval\"):\n self.k = k\n self.models = models\n self.name = name\n \n def evaluate(self, X, y):\n self.fit(X, y)\n return self.error_scores, self.trained_models\n \n def fit(self, X, y, verbose=True):\n # initialise empty lists to store scores and trained models\n self.error_scores = []\n self.trained_models = []\n #self.fold_scores = []\n # split dataset into folds\n folds_X = np.array_split(X, self.k)\n folds_y = np.array_split(y, self.k)\n # iterate over folds\n for test_fold_i in range(self.k):\n # log\n if verbose:\n print(f\"{self.name} fold {test_fold_i+1}: ---------------------------------------------\")\n # specify train/test data\n X_train = folds_X[:test_fold_i] + folds_X[test_fold_i+1:]\n y_train = folds_y[:test_fold_i] + folds_y[test_fold_i+1:]\n X_test = folds_X[test_fold_i]\n y_test = folds_y[test_fold_i]\n # train on everything but the test_fold\n fold_scores = []\n fold_models = []\n for model in self.models:\n # make deep copy of model\n model = copy(model)\n # fit and test\n model.fit(np.concatenate(X_train), np.concatenate(y_train))\n fold_models.append(model)\n score = model.score(X_test, y_test)\n fold_scores.append(score)\n # log\n if verbose:\n print(f\"{model} trained, score = {score}\")\n # save fold scores and models\n self.error_scores.append(np.array(fold_scores))\n self.trained_models.append(np.array(fold_models, dtype=object))\n \n # just neccessary if inner loop:\n # find best model\n scores = np.stack(self.error_scores).mean(axis=0)\n self.best_model = copy(np.array(self.models)[scores == scores.max()][0])\n # fit best model on whole training data\n self.best_model.fit(X, y)\n \n def score(self, X_test, y_test):\n return self.best_model.score(X_test, y_test)","metadata":{"execution":{"iopub.status.busy":"2022-08-12T13:48:02.800302Z","iopub.execute_input":"2022-08-12T13:48:02.800801Z","iopub.status.idle":"2022-08-12T13:48:02.817458Z","shell.execute_reply.started":"2022-08-12T13:48:02.800757Z","shell.execute_reply":"2022-08-12T13:48:02.816170Z"},"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# list of possible alphas\nalphas = np.logspace(0, 3, 5)\n\n# list of resulting models with polinomial\nfeature_func = [lambda x: x, lambda x: 1] #, lambda x: x ** 2, lambda x: x ** 3, lambda x: x ** 4 ...\nmodels = [MultivariateLinearModel('l1', a) for a in alphas]\nfor model in models:\n model.feature_func = feature_func\n \n# gridsearch\nclf = KfoldCrossval(4, models)\nscores, trained_models = clf.evaluate(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2022-08-12T14:11:19.514035Z","iopub.execute_input":"2022-08-12T14:11:19.514486Z","iopub.status.idle":"2022-08-12T14:12:32.585178Z","shell.execute_reply.started":"2022-08-12T14:11:19.514433Z","shell.execute_reply":"2022-08-12T14:12:32.584203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# print out best value for alpha\nprint(f'best tested value for alpha: {clf.best_model.alpha}')","metadata":{"execution":{"iopub.status.busy":"2022-08-12T14:38:33.964825Z","iopub.execute_input":"2022-08-12T14:38:33.965287Z","iopub.status.idle":"2022-08-12T14:38:33.971001Z","shell.execute_reply.started":"2022-08-12T14:38:33.965248Z","shell.execute_reply":"2022-08-12T14:38:33.969872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**final performance evaluation on test set:**","metadata":{}},{"cell_type":"code","source":"# make predictions:\npredictions = clf.best_model.predict(X_test)\n\n# calculate rmse:\nrmse = np.sqrt(np.square(y_test - predictions)).mean()\n\n# calculate R2:\nr2 = 1 - np.square(y_test - predictions).sum()/np.square(y_test - y_test.mean()).sum()\n\nprint(f\"testing data: root mean squared error = {rmse}\\nR-squared = {r2}\")\n\n# make predictions:\npredictions_train = clf.best_model.predict(X_train)\n\n# calculate rmse:\nrmse_train = np.sqrt(np.square(y_train - predictions_train)).mean()\n\n# calculate R2:\nr2_train = 1 - np.square(y_train - predictions_train).sum()/np.square(y_train - y_train.mean()).sum()\n\nprint(f\"training data: root mean squared error = {rmse_train}\\nR-squared = {r2_train}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-12T14:12:58.812046Z","iopub.execute_input":"2022-08-12T14:12:58.812488Z","iopub.status.idle":"2022-08-12T14:12:58.937032Z","shell.execute_reply.started":"2022-08-12T14:12:58.812438Z","shell.execute_reply":"2022-08-12T14:12:58.934965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**That's all for this notebook, have a great day and happy learning!👋**\n\n*If you liked that kind of notebook, I've made quite a handful of similar ones that can be found here: [ml from scratch: table of contents](https://www.kaggle.com/discussions/general/334676)*","metadata":{}}]}