\UseRawInputEncoding
\documentclass[a4paper,12pt]{article}
\usepackage{graphicx} % Required for inserting images
\usepackage[margin=1.5cm]{caption}
\usepackage[T1]{fontenc}  % \ae\o\aa
\usepackage[utf8]{inputenc}
\usepackage{babel}          % ordelingsregler
\usepackage{amsmath, mathtools}
\usepackage{enumitem}
\setlist[itemize]{itemsep=0.3pt}
\usepackage{xcolor}
\usepackage{textgreek}
\usepackage{subfigure}
\usepackage{tabularx}
\usepackage{float}
\usepackage{hyperref}
\usepackage{subcaption} %needed for subfigures
\usepackage{placeins}
\usepackage{comment}
\usepackage{algorithm}
\usepackage{algpseudocode}
\usepackage{hyperref}

\usepackage{listings}
\usepackage{xcolor}

\definecolor{codegreen}{rgb}{0,0.6,0}
\definecolor{codegray}{rgb}{0.5,0.5,0.5}
\definecolor{codepurple}{rgb}{0.58,0,0.82}
\definecolor{backcolour}{rgb}{0.95,0.95,0.92}

\lstdefinestyle{code}{
    language=Python,
    backgroundcolor=\color{backcolour},   
    commentstyle=\color{codegreen},
    keywordstyle=\color{magenta},
    numberstyle=\tiny\color{codegray},
    stringstyle=\color{codepurple},
    basicstyle=\ttfamily\footnotesize,
    breakatwhitespace=false,         
    breaklines=true,                 
    captionpos=b,                    
    keepspaces=true,                 
    numbers=left,                    
    numbersep=5pt,                  
    showspaces=false,                
    showstringspaces=false,
    showtabs=false,                  
    tabsize=2
}

\lstset{style=code}

\title{Project 1 ( FYS-STK4155 )}
\author{Yoan Tutunarov \& Vladislav Foss \& Sebastian Koranda}
\date{September 2025}

\begin{document}

\maketitle

\begin{abstract}
This project investigates three regression methods—Ordinary Least Squares (OLS), Ridge, and LASSO—applied to polynomial fitting of the Runge function. We implement both closed-form solutions and gradient descent approaches, and compare their performance, noting small deviations that are expected due to numerical optimization. To analyze model complexity, we study the bias–variance trade-off using bootstrap resampling, which provides estimates of bias, variance, and noise. The results show that OLS overfits at higher polynomial degrees, Ridge reduces variance by shrinking coefficients, and LASSO promotes sparsity by eliminating irrelevant terms. Overall, the study illustrates how regularization and resampling techniques improve predictive performance and provide deeper insight into the generalization properties of regression models.
\end{abstract}

\section*{Acknowledgments and AI Disclaimer}
We would like to thank our supervisor for valuable feedback and discussions throughout the project.
This project has also made limited use of AI tools, such as ChatGPT, as a support mechanism for improving clarity, refining explanations, and troubleshooting code. The primary analysis, problem solving, and discussion were carried out independently by the project members. AI assistance was used only to enhance language in the report, verify understanding, and used for polishing figure names, ensuring that all team members fully comprehend the underlying concepts and methods.

\section*{Introduction}
In this project, we study different regression methods with a particular focus on Ordinary Least Squares (OLS), Ridge regression, and Lasso regression. These methods form the foundation of many statistical and machine learning approaches, and they provide useful insights into the trade-off between model complexity and predictive performance.
As a test case, we apply these methods to the Runge function, a one-dimensional function that is well known to illustrate the challenges of polynomial interpolation, especially at higher polynomial degrees (Runge’s phenomenon). By fitting polynomials of varying complexity to this function, we investigate how different regression techniques handle overfitting, variance, and bias.
The project is based on weekly exercises, lecture material, and additional references. All numerical experiments are implemented in Python, with the code provided in a separate file and linked through GitHub. This report serves as documentation of our implementation, results, and analysis.

\section*{Methods}
\subsection*{Statistical methods}
\subsection*{Mean Squared Error (MSE)}

Mean Squared Error (MSE) is one of the most fundamental measures of model performance in statistics and machine learning. It quantifies how well a predictive model approximates the true data by measuring the average squared difference between the actual (target) values and the predicted values \cite{GeeksMSE}. 

Formally, for a dataset with $n$ samples, the MSE is defined as
\[
\text{MSE} = \frac{1}{n}\sum_{i=1}^n \left(y_i - \hat{y}_i\right)^2,
\]
where
\begin{itemize}
    \item $y_i$ denotes the true (actual) value of the response variable for observation $i$,
    \item $\hat{y}_i$ denotes the predicted value from the model for the same observation,
    \item $n$ is the total number of observations.
\end{itemize}
The term $\left(y_i - \hat{y}_i\right)$ is called the residual, which represents the error between the actual and predicted values for observation $i$. Squaring the residual ensures that negative and positive deviations are treated equally and emphasizes larger errors. Averaging over all $n$ samples provides a single metric that reflects the model's overall predictive accuracy. 

A lower MSE indicates that the model predictions are closer to the actual values, while a higher MSE suggests larger discrepancies. Therefore, in regression problems, the objective is typically to minimize the MSE, ideally approaching zero.





\subsubsection*{$\mathbf{R^2}$ score}
The \(R^2\) score (often called the coefficient of determination) is a standard metric in regression analysis that quantifies how well the variation in the independent variables explains the variation in the dependent variable. In simpler terms, \(R^2\) measures the proportion of the variance in the target \(y\) that is captured by the regression model, relative to a baseline model that always predicts the mean of \(y\) \cite{GeeksRSquared}.

Formally, let \(y_i\) be the actual observed values, \(\hat y_i\) the predicted values, and \(\bar y = \frac{1}{n}\sum_{i=1}^n y_i\) the sample mean. Then
\[
\mathrm{SSE} = \sum_{i=1}^n (y_i - \hat y_i)^2, \quad
\mathrm{SST} = \sum_{i=1}^n (y_i - \bar y)^2,
\]
and
\[
R^2 = 1 - \frac{\mathrm{SSE}}{\mathrm{SST}} 
= 1 - \frac{\sum_{i=1}^n (y_i - \hat y_i)^2}{\sum_{i=1}^n (y_i - \bar y)^2}.
\]

Here:
\begin{itemize}
  \item \(\mathrm{SSE}\) (“Sum of Squared Errors”) is the sum of squared deviations between predictions and actual values.
  \item \(\mathrm{SST}\) (“Total Sum of Squares”) is the sum of squared deviations of actual values around their mean.
  \item \(R^2\) ranges (in theory) from \(-\infty\) to 1. An \(R^2 = 1\) indicates perfect prediction (all points lie exactly on the regression line). An \(R^2 = 0\) indicates that the model is no better than the constant mean predictor. Negative \(R^2\) occurs when the model performs worse than simply predicting the mean.
\end{itemize}
A higher \(R^2\) implies a better-fitting model in that it reduces residual error relative to the variance of the data itself.


\subsection*{Regression}
Regression is a statistical method used to study the relationship between a dependent variable (outcome) and one or more independent variables (predictors). It estimates how changes in predictors are associated with changes in the outcome. In its simplest form, linear regression fits a straight line through data to best explain this relationship. More complex forms include multiple regression (many predictors) and nonlinear regression (curved relationships). Regression is widely used for prediction, trend analysis, and decision-making in fields like economics, finance, and machine learning \cite{Beers2025}.


\subsubsection*{Ordinary Least Squaress (OLS)}
Ordinary Least Squares (OLS) is a method used to estimate the coefficients of a linear regression model. It finds the values of the coefficients that minimize the total squared distance between the observed data points and the values predicted by the model. Mathematically, OLS solves:

\[
\min_{\beta} \; (Y - X\beta)'(Y - X\beta)
\]

where \(Y\) is the vector of observed values, \(X\) is the matrix of predictors, and \(\beta\) is the vector of coefficients. OLS relies on certain conditions: the errors should be independent, have constant variance, and ideally follow a normal distribution. The method can struggle when predictors are highly correlated or when the dataset is too small compared to the number of variables. Once estimated, OLS coefficients let us assess model quality (e.g., using \(R^2\)) and make predictions~\cite{XLSTAT2025}.
In Task A, we implemented our own OLS function to compute the optimal $\hat{\theta}$. From the lecture notes in Week 35, we know that the solution to the regression problem is given as $y = X\theta$, where $X$ is the design matrix. In this assignment we used a polynomial of degree 15, meaning that $\hat{\theta}$ is a one-dimensional array with 15 parameters.  

The solution for OLS is obtained by differentiating the cost function and solving for its minimum:
\begin{align*}
    \frac{\partial C}{\partial \theta} = 0 &= X^\top (y - X\theta), \\
    X^\top y &= X^\top X \theta,
\end{align*}
and given that $X^\top X$ is invertible, the solution becomes
$$
\hat{\theta} = (X^\top X)^{-1}X^\top y.
$$

Here $y$ is given by Runge's function $y = \frac{1}{1+25x^2}$ with added noise. The noise was generated using NumPy’s \texttt{np.random.normal()}, with parameters specified in the project description.  

We then constructed the design matrix $X$ and split the data into training and test sets using scikit-learn’s built-in \texttt{train\_test\_split()} function.

\subsubsection*{Ridge regression}
Ridge Regression (also called L2 regularization) is a variation of linear regression that introduces a penalty on large coefficient values in order to reduce model variance and improve stability in the presence of multicollinearity. It adjusts the ordinary least squares estimator by adding a term \(\lambda I\) (where \(\lambda \ge 0\) is the regularization parameter and \(I\) is the identity matrix) such that the coefficient estimates are given by

\[
\hat{\beta}_\lambda = \bigl(X^\top X + \lambda I\bigr)^{-1} X^\top y
\]

This estimator shrinks the coefficients toward zero, but unlike some other methods, does not set any coefficients exactly to zero. The choice of \(\lambda\) controls the bias-variance tradeoff: a larger \(\lambda\) increases bias but lowers variance and may improve prediction for new data. Care must be taken: if \(\lambda\) is too large, the model may underfit; if too small, it may behave similarly to OLS and still overfit~\cite{GFGRidge2025}.


\subsubsection*{LASSO Regression}
LASSO (Least Absolute Shrinkage and Selection Operator) regression is an extension of linear regression that introduces both regularization and variable selection. It modifies the ordinary least squares (OLS) objective by adding an $L_1$ penalty term, which constrains the absolute values of the model coefficients. The LASSO objective function is defined as

\[
\min_{\boldsymbol{\beta}} \left\{ 
\frac{1}{2n} \sum_{i=1}^{n} \left( y_i - \mathbf{x}_i^\top \boldsymbol{\beta} \right)^2 
+ \lambda \sum_{j=1}^{p} |\beta_j| 
\right\},
\]
where $\lambda \ge 0$ is a regularization parameter controlling the strength of the penalty.

The $L_1$ regularization term encourages sparsity in the coefficient vector $\boldsymbol{\beta}$ by shrinking some coefficients exactly to zero, thereby performing feature selection. This simplifies the model and helps manage the bias–variance trade-off: as $\lambda$ increases, model variance decreases while bias increases. However, if $\lambda$ is too large, the model may underfit and fail to capture important data patterns.

Unlike OLS and Ridge regression, LASSO does not admit a closed-form analytical solution because the $L_1$ norm introduces a non-differentiable point at zero. Consequently, iterative optimization algorithms such as coordinate descent or gradient descent are required to estimate the parameters \cite{geeksforgeeks_lasso}.



\subsection*{Gradient Descent Variants}

In machine learning, we typically have a dataset $X$ and a model defined by parameters $\theta$ (sometimes also denoted $\beta$). 
The objective is to determine the parameter values that minimize a chosen cost function $C(\theta)$. 
Gradient descent is an iterative optimization algorithm designed to solve this problem.

The main idea is the following: if we have a function $F(x)$, with $x = (x_1, x_2, \dots, x_n)$, then the gradient $\nabla F(x)$ points in the direction of the steepest ascent of the function. 
By moving in the opposite direction of the gradient, the algorithm approaches a minimum of the cost function after taking enough steps. 
The general update rule for gradient descent is
\begin{align*}
    \mathbf{x}_{k+1} = \mathbf{x}_k - \eta_k \nabla F(\mathbf{x}_k),
\end{align*}
where the parameter $\eta$ is the \emph{learning rate}. 
This can be thought of as the step size taken in the chosen direction. 
If $\eta$ is too small, convergence towards the minimum is very slow, while if $\eta$ is too large, the algorithm may diverge. 
Choosing a suitable learning rate is therefore of great importance.

In our project, the goal is to minimize the Mean Squared Error (MSE) cost function for OLS, and MSE with an additional penalty term for Ridge regression. 
Minimization is performed by iteratively adjusting the model parameters $\theta$, thereby reducing the error between the predicted values and the true outcomes \cite{IBMGradientDescent}. 
We start by initializing all coefficients to zero. At each iteration, we compute the gradient of the cost function, which tells us how the error changes with respect to each parameter. 
The gradient indicates the slope of the loss surface: a positive gradient means the parameter should be decreased, while a negative gradient means it should be increased. 
Using a fixed learning rate $\eta$, we then update the coefficients according to
\[
\theta \; \leftarrow \; \theta - \eta \nabla_\theta C(\theta).
\]

For OLS, the gradient is 
\[
\nabla_\theta C_{\text{OLS}}(\theta) = \frac{2}{n} X^\top (X\theta - y),
\]
while for Ridge it becomes
\[
\nabla_\theta C_{\text{Ridge}}(\theta) = \frac{2}{n} X^\top (X\theta - y) + 2\lambda \theta,
\]
where $\lambda$ is the regularization parameter \cite{CompPhysicsWeek36}.
 

\subsubsection*{Stochastic Gradient Descent (SGD)}

In the most basic version of the gradient descent, we usually include the whole dataset to compute the the gradient descent at each stop. This is in contrast to \textbf{SGD}, where we would calculte the gradient on a subset of data called mini-batches. If our whole dataset is \textbf{n} datapoints, and we want to have a mini-batch of size \textbf{M}, then the number of batches would be $\mathbf{\frac{n}{M}}$.

Compared to the full batch (whole dataset) gradient descent, SGD typically will converge faster per iteration, but will however experience larger variance due to noise and the random nature of the method. A trade off with the SGD is the size of each batch: larger batches will be computationally more expensive but less noisy, and smaller batches will be very cost efficient, however will have a lot more noise in the trajectory they are calculating \cite{Week37SGD}.


The gradient step of the SGD is calculated as:
\begin{align*}
    \theta_{t+1} = \theta_{t} - \eta_t \nabla C_{B_k}(\theta_t),
\end{align*}

Here $k$ represents the randomly chosen mini-batch from the dataset at each iteration. 
One complete pass over all mini-batches is referred to as an \emph{epoch}~\cite{Week37SGD}. 
Calculating the gradient on subsets of the data reduces computational cost compared to full-batch gradient descent, 
and the inherent randomness can also help the algorithm escape shallow local minima and saddle points.
\\
\\
Since the learning rate plays a critical role in the performance of SGD, 
a variety of adaptive learning rate algorithms have been developed. 
These methods adjust the step size dynamically based on information from previous updates. 
In this project, we focus on four  variants: \textbf{Momentum}, \textbf{AdaGrad}, \textbf{RMSProp}, and \textbf{Adam}, 
each of which introduces specific improvements to the basic GD and SGD functions.

\subsubsection*{Momentum based GD}

SGD generally today is implemented with momentum. The general idea behind momentum based GD / SGD is that if the landscape we are traversing with our gradient is a canyon, the vanilla gradients will have large oscillations instead of gliding towards the minima. The implementation of momentum build on the idea that if the direction of the gradient descent aligns with the previous point, then we gain "momentum", and if not we dampen it. This way we dampen the oscillations that we might experience. This is mathematically described as:
\begin{align*}
    \theta_{t+1} = \theta_t - v_t , \ \text{where} \\
    v_t = \gamma v_{t-1} + \eta_t \nabla E(\theta_t)
\end{align*}
The function $E(\theta)$ represents our function we are trying to find the minimum of \cite{Week37SGD}. Here we can see that we will get a smoothed path where we accelerate if the gradient is still pointing in the same direction. We also introduce a new parameter $\gamma$, which represents the momentum, and we usually choose between $0 \leq \gamma \leq 1$ \cite{Week37SGD}.



\subsubsection*{Adagrad}

The adaptive learn rate of Adagrad drops the parameter $\gamma$ represented in momentum completely and has a different approach. For this algorithm, we first define a vector:
\begin{align*}
    r_t = r_{t-1} + (\nabla C(\theta))^2
\end{align*}
Calculating the AdaGrad at step t will then become 
\begin{align*}
    \theta_{t+1} = \theta_t -\eta H_t^{-1/2}\nabla C(\theta)
\end{align*}
, where we have that the matrix $H_t$ is defined as $H_t = \text{diag}(r_t)$, where the entries in the matrtix will diagonally be entries of $r_t$ over time. The update of $\theta$ will be:
\begin{align*}
    \theta_{t+1 , j} = \theta_{t , j} - \frac{\eta}{\sqrt{\epsilon + r_{t , j}}}\nabla C(\theta_{t , j})
\end{align*}
Here we introduce the parameter $\epsilon$ as a small constant to avoid division by zero if $r_{t , j}$ is 0. From the mathematical expressions, we can see that there we have removed the momentum ($\gamma$) as implemented in momentum gradient descent described in a section above. We see that if a gradient is very large the step size we take ($\frac{\eta}{\sqrt{\epsilon + r_{t , j}}}$) will automatically become very small, meaning we have less large jumps, while if the gradient is very small the step size becomes larger and speeds up \cite{Week37SGD}.
\\
\\
Unfortunately, a key limitation of AdaGrad is that the accumulated squared gradients in $H_t$ grow without bound. 
This causes the algorithm to take very large steps in the beginning, but as the vector $r_{t , j}$ grows, the effective learning rate
$$
-\frac{\eta}{\sqrt{r_{t , j} + \epsilon}}
$$
becomes smaller and smaller \cite{Week37SGD}. 
As a result, the algorithm slows down excessively, and the step sizes eventually become infinitesimally small, effectively stopping further learning. 
This limitation is the main motivation behind the development of RMSProp and Adam, which modify AdaGrad to maintain a more stable learning rate~\cite{Week37SGD}.

\subsubsection*{RMSProp}
Subsequently, the root-mean-square propagation (RMSProp) algorithm was developed to address AdaGrad's limitation of diminishing learning rates that eventually lead to extremely small step sizes. 
The key idea is to use an exponentially decaying moving average of past squared gradients rather than accumulating all of them. 
This acts as a weighted average that gives more importance to recent gradients than to older ones, preventing the learning rate from decaying towards zero over many iterations~\cite{Optmizationgradient,Week37SGD}.

The exponentially decaying average is defined as
$$
v_t = \rho v_{t-1} + (1-\rho)(\nabla C(\theta_t))^2 ,
$$
where $\rho$ is the decay rate, typically chosen in the range $0.9 \leq \rho \leq 0.99$.  

The update rule for the parameters then becomes
$$
\theta_{t+1} = \theta_t - \frac{\eta}{\sqrt{v_t + \epsilon}} \nabla C(\theta_t),
$$
where $\epsilon$ is a small constant included to avoid division by zero. 

From these expressions, we see that recent gradients have greater influence than earlier ones, ensuring that the effective learning rate remains stable. 
This modification makes RMSProp more robust than AdaGrad and it has become a widely used optimizer in practice, as well as an important component of the Adam algorithm \cite{Optmizationgradient,Week37SGD}.


\subsubsection*{Adam}
The adaptive moment estimation (Adam) algorithm combines the ideas of RMSProp, with its adaptive learning rate, and momentum, which smooths the optimization path. 
Adam maintains running averages of both the first and the second moments of the gradients~\cite{Week37SGD}.

The first moment is defined as
$$
m_t = \beta_1 m_{t-1} + (1-\beta_1)\nabla C(\theta_t),
$$
and the second moment as
$$
v_t = \beta_2 v_{t-1} + (1-\beta_2)(\nabla C(\theta_t))^2,
$$
where the parameters are typically chosen as $\beta_1 = 0.9$ and $\beta_2 = 0.999$, 
with initial values $m_0 = 0$ and $v_0 = 0$~\cite{Week37SGD}.  

When we first calculate the first and second moments ($m_t, v_t$), both are biased towards zero in the early iterations. 
This is because the running averages favor the most recent gradients, and thus require time to build up from zero. 
To correct for this bias, Adam introduces normalized versions of the moments. The bias-corrected terms are given as~\cite{Week37SGD}:
$$
\hat{m}_t = \frac{m_t}{1-\beta_1^t}, \qquad 
\hat{v}_t = \frac{v_t}{1-\beta_2^t}.
$$

The parameter update rule then becomes
$$
\theta_{t+1} = \theta_t - \frac{\eta}{\sqrt{\hat{v}_t} + \epsilon}\hat{m}_t ,
$$
where $\eta$ is the learning rate and $\epsilon$ is a small constant to avoid division by zero.  

From these expressions we see that in the first few iterations, when $t$ is small, the denominators $(1-\beta_i^t)$ are close to zero, leading to significantly larger values of $\hat{m}_t$ and $\hat{v}_t$ compared to the uncorrected terms. 
This correction is critical for stability in the early iterations of the algorithm~\cite{Week37SGD}. 
As the number of iterations increases, the bias-corrected estimates converge towards the uncorrected values.  

By combining the exponential moving averages from RMSProp with the smoothing effect of momentum, Adam achieves fast convergence with relatively little tuning. 
It is also robust when dealing with sparse gradients, making it one of the most widely used optimization algorithms in modern machine learning~\cite{Week37SGD}.

\subsection*{Bias-Variance}
An important concept in machine learning and statistical analysis is the bias-variance trade-off. 
We can decompose the expected prediction error of a model into three components: bias, variance, and irreducible error (noise).
This means that the expected squared error between the true values and the model's predictions can be expressed as the sum of these three components.

Knowing this, and the definition of mean squared error (MSE), we have:
$$
\mathbb{E}\left[(\boldsymbol{y}-\boldsymbol{\tilde{y}})^2\right]=\mathrm{Bias}[\tilde{y}]+\mathrm{var}[\tilde{y}]+\sigma^2,
$$
where
$$
\mathrm{Bias}[\tilde{y}]=\mathbb{E}\left[\left(\boldsymbol{y}-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right]\right)^2\right],
$$
and
$$
\mathrm{var}[\tilde{y}]=\mathbb{E}\left[\left(\tilde{\boldsymbol{y}}-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right]\right)^2\right]=\frac{1}{n}\sum_i(\tilde{y}_i-\mathbb{E}\left[\boldsymbol{\tilde{y}}\right])^2.
$$

We can derive this decomposition as follows. Assume the data is generated according to the model:
$$
y_i = f(x_i) + \varepsilon_i,\qquad \mathbb{E}[\varepsilon_i]=0,\quad \mathrm{Var}(\varepsilon_i)=\sigma^2,
$$
and let $\tilde{y}_i$ be the predictions of the model
$$
\mathbb{E}\!\left[(y_i-\tilde{y}_i)^2\right]
= \mathbb{E}\!\left[(f(x_i)+\varepsilon_i-\tilde{y}_i)^2\right].
$$

Add and subtract $\mathbb{E}[\tilde{y}_i]$:
$$
\begin{aligned}
\mathbb{E}\!\left[(y_i-\tilde{y}_i)^2\right]
&= \mathbb{E}\!\left[\big(f(x_i)-\mathbb{E}[\tilde{y}_i] + \mathbb{E}[\tilde{y}_i]-\tilde{y}_i + \varepsilon_i\big)^2\right] \\
&= \underbrace{\big(f(x_i)-\mathbb{E}[\tilde{y}_i]\big)^2}_{\text{bias}^2}
+ \underbrace{\mathbb{E}\!\left[(\tilde{y}_i-\mathbb{E}[\tilde{y}_i])^2\right]}_{\mathrm{Var}(\tilde{y}_i)}
+ \underbrace{\mathbb{E}[\varepsilon_i^2]}_{\sigma^2} \\
&\quad +\; \underbrace{2\,(f(x_i)-\mathbb{E}[\tilde{y}_i])\,\mathbb{E}[\tilde{y}_i-\mathbb{E}[\tilde{y}_i]]}_{=\,0}
\;+\; \underbrace{2\,\mathbb{E}[(\tilde{y}_i-\mathbb{E}[\tilde{y}_i])\,\varepsilon_i]}_{=\,0}.
\end{aligned}
$$

The cross-terms vanish because $\mathbb{E}[\tilde{y}_i-\mathbb{E}[\tilde{y}_i]]=0$ and, at test time, $\tilde{y}_i$ is independent of $\varepsilon_i$ with $\mathbb{E}[\varepsilon_i]=0$. Therefore,
$$
\boxed{\;
\mathbb{E}\!\left[(y_i-\tilde{y}_i)^2\right]
= \big(f(x_i)-\mathbb{E}[\tilde{y}_i]\big)^2 \;+\; \mathrm{Var}(\tilde{y}_i) \;+\; \sigma^2
\;}
$$

Averaging over all points $i=1,\dots,n$ yields the dataset-level decomposition:
$$
\boxed{\;
\frac{1}{n}\sum_{i=1}^n \mathbb{E}\!\left[(y_i-\tilde{y}_i)^2\right]
= \frac{1}{n}\sum_{i=1}^n \big(f(x_i)-\mathbb{E}[\tilde{y}_i]\big)^2
\;+\; \frac{1}{n}\sum_{i=1}^n \mathrm{Var}(\tilde{y}_i) \;+\; \sigma^2
\;}
$$











\section*{Implementation}

\subsection*{Measurements for error}

\subsubsection*{Mean Squared Error (MSE)}
For calculating the Mean Squared Error (MSE), we make use of the \href {https://scikit-learn.org/stable/modules/generated/sklearn.metrics.mean_squared_error.html}{scikit-learn} library’s built-in function 
\texttt{mean\_squared\_error}. This ensures both efficiency in our implementations


\subsubsection*{$R^2$ score}
The $R^2$ score is implemented from its theoretical definition. 
The numerator and denominator are calculated separately, and then combined into the final return, as shown in Code~\ref{code:r2}.
This metric will be used to evaluate the performance of our models later.

\subsection*{Scaling}
On the topic of scaling, we usually scale our data before performing regression. 
This ensures that the different features $\theta$ are on comparable scales, so that no single feature 
dominates or is weighted less simply because of the magnitude of the feature. 
As a general practice in our implementations (see Listings~\ref{code:adaptive}, \ref{code:ols}, and \ref{code:ridge}), 
we include an explicit scaling step.

The standard way we scale the data is to call the \texttt{StandardScaler} 
function from the Scikit-learn library~\cite{sklearnStandardScaler}, 
and create it as an instance called \textit{scaler}. 
We then scale the training data and apply the same fitted scaler to the test data 
to avoid data contamination between the train and test split.



\subsection*{Regression Models}

  \subsubsection*{ OLS for the Runge function}

    We begin by generating a dataset of 100 equally spaced points 
    \(x \in [-1,1]\). For each point, the response value \(y\) is calculated 
    from the Runge function, with added random noise to simulate measurement errors.  
    To prepare the data for regression, we implemented the function 
    \texttt{polynomial\_features}. This function constructs the design matrix for a chosen polynomial degree \(p\). Each row in the design matrix contains the powers of \(x\), i.e. 
    \([1, x, x^2, \dots, x^p]\), with the option of including or excluding an intercept term.  
    The regression parameters are estimated with our function 
    \texttt{OLS\_parameters}, which implements the ordinary least squares (OLS) solution. 
    Model performance is measured using two standard metrics: the mean squared error (MSE) 
    and the coefficient of determination (\(R^2\)), both coded directly from their mathematical definitions 
    (\texttt{R2\_score}).  
    In the main script, the design matrix is created and the dataset is split into training (80\%) 
    and test (20\%) sets. The polynomial degree is varied from 1 to 15. For each degree, the following steps are carried out:   
    \begin{enumerate}
        \item Construct the design matrix with \texttt{polynomial\_features}.
        \item Split the dataset into training and test sets.
        \item Compute the OLS coefficients from the training data.
        \item Use the coefficients to predict values on the test set.
        \item Compute MSE and \(R^2\) for the predictions.
        \item Store the values for later comparison.
    \end{enumerate}
    Finally, we create an array of polynomial degrees and plot both MSE and \(R^2\) 
    as functions of the polynomial degree. These plots provide a visual overview of how 
    model complexity affects prediction accuracy.  
    This implementation directly reflects the theoretical framework from the lectures: 
    polynomial features control model complexity, the OLS solution provides a baseline estimator, 
    and MSE and \(R^2\) serve as performance measures. Gradient-based optimization methods 
    (introduced later in the course) could be applied as an alternative, but here we focus 
    on the analytical OLS approach.  



  
  \subsubsection*{Ridge}
Based on the theoretical definition, we implement Ridge regression as a function 
\texttt{ridge\_regression(X, y, $\lambda$)}. 
The parameters are estimated by
$$
\theta = (X^\top X + \lambda I)^{-1}X^\top y.
$$
To be able to study the effect of the regularization parameter $\lambda$, we selected multiple values of $\lambda$ using \texttt{train\_test\_split} to divide data into training and test sets. 
For each $\lambda$, the model was fitted on training data and evaluated on the test set using mean squared error (MSE) and $R^2$. 
The implementation is provided in Code Listing~\ref{code:ridge}. We then subsequently plot the \textbf{MSE} values along with the \textbf{$R^2$} values with the different values for $\lambda$.
\\
\\
We also implemented a way to gauge the difference between the \textbf{MSE} and $R^2$ values between the Ridge and OLS implementations, with differing $\lambda$ values. The way this as done was by computing the mean of the \textbf{MSE} values for each of the $\lambda$ values, and then plotting them against eachother, as we can see in the code \ref{code:mser2olsridge}


  \subsubsection*{Lasso}



  
\subsection*{Gradient Descent Variants}
After running the algorithm for a sufficient number of iterations, the coefficients converge close to the closed-form solutions we derived in Parts A and B. Comparing the two approaches shows that gradient descent gives nearly the same coefficients, with small numerical differences due to the finite number of iterations and the choice of learning rate. Ridge regression, as expected, shrinks the coefficients compared to OLS.

A key observation from experimenting with the learning rate $\eta$ is that it strongly influences convergence: if $\eta$ is too small, the algorithm converges very slowly; if it is too large, the updates overshoot and the algorithm can diverge. Choosing a moderate value allows stable and efficient convergence.

\cite{Optmizationgradient}.

\subsubsection*{Momentum GD}
We implement the momentum gradient descent as a standalone function to allow flexibility and reuse in other parts of the project. 
The algorithm is initialized with $\theta = 0$, and a variable \texttt{change} is introduced to store the previous update direction. 
At each iteration, the gradient is computed depending on whether we solve the OLS, Ridge, or Lasso case. For our momentum, we select a parameter of our choosing.
We also track the mean squared error (MSE) at each iteration to monitor convergence, as shown in code listing \ref{code:adaptive}.

\subsubsection*{AdaGrad}
We implement the Adagrad also as a standalone function to allow flexibility and the ability to reuse in other parts of the project. We initialize $\theta = 0$ and also $\mathbf{r} = 0$. For gradient calculation we can easily choose OLS or Ridge in our function. Then we square the gradient and add it to the previous value of $\mathbf{r}$. We calculate the square of the gradient, the adaptive step and the \textbf{MSE} value, and return the $\theta, \mathbf{MSE}$. The code for the Adagrad implementation is seen in \ref{code:adaptive}.


\subsubsection*{RMSProp}
RMSProp was also implemented as a standalone function as well, allowing for reussability and flexibility. We start by initializing the parameters $\theta , \mathbf{v} = 0$, followed by setting up a matrix for storing our \textbf{MSE} values and calculating the decaying average of gradients. We then calculate the parameter \textbf{v}, and update the $\theta$ value.The code for the RMSProp implementation can be seen in \ref{code:adaptive}.




\subsubsection*{ADAM}
We also implement ADAM as a standalone function, where we first initialize the parameters $\theta , \mathbf{v} , \mathbf{m} = 0$. We then loop over the iterations and calculate the 1st and 2nd moment, but also creating parameters $\mathbf{\hat{m}} , \mathbf{\hat{v}}$ for the bias corrections, and using these in the final update for $\theta$. The code for the ADAM implementation can be seen in \ref{code:adaptive}.




\subsubsection*{Stochastic Gradient Descend (SGD)}
We imeplented a Stochastic gradient descent (SGD) as a standalong function with different adaptive learning rates that we tested, along with \textbf{OLS} and \textbf{Ridge} regression methods. The implementation follows the principles laid out in the section about the SGD method, further above.
\\
\\
At each epoch, the data is shuffled randomly to make sure the batches are not correlated between the updates. Then we have the minibatches being created with the chosen length of the batches. For every minibatch, the gradient is computed based on the chosen style, and further, a learn rate optimizer is being used. The ones that are in question are the \textbf{Momentum} , \textbf{AdaGrad} , \textbf{RMSProp} and \textbf{ADAM}. The function then inside of these optimizers caluclated the $\theta$ update, and also calculates the \textbf{MSE} values for the minibatch. This way we can check convergence later in plotting.
\\
\\
The learn rate was set to $\eta=0.03$, with a momentum factor of $0.3$, as well as the $\lambda$ being able to be specified. The reasons for setting these parameters was that they during testing provided most stable plots and convergences.
The full implementation of the SGD is shown in code listing~\ref{code:SGD}







\subsection*{Resampling}
  \subsubsection*{Bootstrap setup}
    To perform bootstrap resampling, we call the \texttt{resample} utility function of scikit-learn iteratively on the training data for the number of bootstraps given as argument. 
    This function samples the dataset, with replacement, meaning it returns a dataset of the same size, consisting of elements randomly selected from the original dataset while allowing duplicates. 
    Essentially, we get variations of the original data with some points repeated and others omitted. Repeating this process multiple times generates a distribution of datasets that reflect the variability inherent in the original data.

    For each iteration, we fit the model using the resampled data and perform a prediction using each resulting fit.
    After completing all bootstrap iterations, we compute the mean and standard deviation of the predictions across all bootstrap samples. 
    This provides an estimate of the uncertainty in our model's predictions, as the variability across the bootstrap samples reflects the sensitivity of the model to changes in the training data.

    Bootstrap resampling is implemented as a method of the OLS class in Appendix B~\ref{code:OO-OLS}
  \subsubsection*{Cross-validation setup}
    For cross-validation, we use the \texttt{KFold} function from scikit-learn to split the data into \(k\) folds. For each fold, we fit the model on \(k-1\) folds and validate it on the remaining fold. We repeat this process for all folds and compute the average performance metrics (MSE and \(R^2\)) across all folds to assess the model's generalization ability.

    K-folds cross validation is also implemented for the OLS class~\ref{code:OO-OLS}



\section*{Results}

\subsection*{Part A: Ordinary Least Square (OLS) for the Runge function}
We fitted the noisy Runge function $f(x)=1/(1+25x^2)$ using OLS
with polynomial features up to degree~15. We tested both the scaled an unscaled versions of the feature matrix to evaluate the effect of scaling the data. The results are from the code listed in \ref{code:ols}. The results can be seen in figure~\ref{fig:part-a}.
\begin{figure}[h!]
    \centering
    \includegraphics[width=1\linewidth]{part-a-ols-ridge.png}
    \caption{OLS and Ridge regressions. Scaled and unscaled versions}
    \label{fig:part-a}
\end{figure}


\subsection*{Part B: Adding Ridge regression for the Runge function  }
We implemented the \textbf{Ridge} regression model and fitted it to the noisy Runge function 
$f(x) = 1 / (1 + 25x^2)$ using different values of the regularization parameter~$\lambda$. 
The tested values were
$$
\lambda = [0.0001,\ 0.01,\ 0.1,\ 0.5].
$$
The data were scaled prior to fitting, and the full implementation and plotting code 
are provided in Code Listing~\ref{code:ridge}. 
The resulting fits are shown in Figure~\ref{fig:ridge-1}.
\begin{figure}[h!]
    \centering
    \includegraphics[width=1\linewidth]{part-b-1.png}
    \caption{Ridge Analysis for different polynomial degrees}
    \label{fig:ridge-1}
\end{figure}
\\
\\
We also plotted the MSE values of OLS vs Ridge by taking the average of the values for the dataset, and then subtracting the differences between them for each of the $\lambda$ parameters. The results can be seen in figure \ref{fig:olsridgemse}, and the code for implementation can be found in code listing \ref{code:mser2olsridge}.
\begin{figure}[h!]
    \centering
    \includegraphics[width=1\linewidth]{part-b-2.png}
    \caption{OLS vs Ridge comparisson in MSE}
    \label{fig:olsridgemse}
\end{figure}








\subsection*{Part C: Writing our own Gradient descent code}
\begin{figure}[H]
    \centering 
    \includegraphics[width=\textwidth]{gd1.png} 
    \caption{GD Convergence} 
    \label{fig:GD_1} 
    \end{figure}

\begin{figure}[H]
    \centering 
    \includegraphics[width=\textwidth]{gd2.png} 
    \caption{GD: Model prediction vs True function} 
    \label{fig:GD_2} 
    \end{figure}


\subsection*{Part D: Including momentum and more advanced ways to update the learning the rate }
For this assignment, we applied the different adaptive gradient descent optimizers 
(\textbf{Momentum}, \textbf{AdaGrad}, \textbf{RMSProp}, and \textbf{Adam}) 
to both \textbf{OLS} and \textbf{Ridge} regression on the noisy Runge function 
$f(x) = 1 / (1 + 25x^2)$ using polynomial features of degree~6. In testing this was proven to give the best fit.  
The models were trained on the training data (80/20 split), and the fitted parameters were 
then applied to the entire dataset to visualize the final approximation to the true function.

The learning rate was set to $\eta = 0.03$ with a momentum factor of $0.4$, 
and the parameter $\lambda$ for Ridge regression was $\lambda = 1\times10^{-4}$.  
For each optimizer, we also recorded the mean squared error (MSE) at each iteration 
to study convergence.

The complete implementation is provided in Code Listing~\ref{code:adaptive}, 
and the resulting fits and convergence curves are shown in Figure~\ref{fig:adaptive-results}.
\begin{figure}[h!]
    \centering
    \includegraphics[width=\linewidth]{part-d.png}
    \caption{Optimizers fitted on the Runge Function}
    \label{fig:adaptive-results}
\end{figure}


\subsection*{Part E: Writing own code for Lasso regression}
\begin{figure}[H]
    \centering 
    \includegraphics[width=\textwidth]{lasso1.png} 
    \caption{Lasso Convergence} 
    \label{fig:lasso} 
    \end{figure}

\begin{figure}[H]
    \centering 
    \includegraphics[width=\textwidth]{lasso2.png} 
    \caption{Effect of regularization strength on Lasso} 
    \label{fig:lasso_2} 
    \end{figure}


\subsection*{Part F: Stochastic Gradient Descent}

We implemented \textbf{SGD} and tested the adaptive optimizers 
\textbf{Momentum}, \textbf{AdaGrad}, \textbf{RMSProp}, and \textbf{Adam} 
for both \textbf{OLS} and \textbf{Ridge} regression on the Runge function. 
The models were trained using mini-batches with an 80/20 train–test split, 
a polynomial degree of~10, and a learning rate of~$\eta = 0.03$. 
The results were compared with the baseline \textbf{vanilla gradient descent (GD)}.

For each optimizer, we recorded the mean squared error (MSE) at each iteration 
to study convergence behaviour. 
The implementation of the SGD procedure is provided in Code Listing~\ref{code:SGD}, 
and the resulting convergence curves are shown in Figure~\ref{fig:sgd}.

\begin{figure}[h!]
    \centering
    \includegraphics[width=\linewidth]{SGD-convergence.png}
    \caption{SGD convergence for OLS (left) and Ridge (right) regressions 
    compared with vanilla gradient descent.}
    \label{fig:sgd}
\end{figure}
\newpage

\subsection*{Part G: Bias-Variance Trade-off with Bootstrap Resampling}
\begin{figure}[h!]
    \centering
    \includegraphics[width=\linewidth]{bias_variance_trade_off.png}
    \caption{Total Error, Bias, and Variance, with 100 bootstrap resamplings for an Ordinary Least Squares (OLS) model on 100 data points.}
    \label{fig:bias-variance}
\end{figure}
\newpage

\subsection*{Part H: Cross-Validation}
\begin{figure}[h!]
    \centering
    \includegraphics[width=1\linewidth]{cross_validation.png}
    \caption{5-fold cross validation on our three linear regression models, optimized using fixed learning rate, full-batch gradient descent for 100 data points. Values of regularization hyperparameters specified in legend.}
    \label{fig:cross-validation}
\end{figure}

\section*{Discussion}
\subsection*{Part A}
As we can see from the results in  figure \ref{fig:part-a}, we observe that generally the \textbf{MSE} values decrease as and the $R^2$ values increase in the unscaled version with increasing polynomial degrees. This should not perplex us given that a more complicated model can acount for more of the variation in the dataset than a more simple one. We can see that the lowest values for the error are given with a degree 6 and onward in the unscaled version of the \textbf{OLS}.
\\
\\
When observing the scaled data (2 bottom images): we can observe that at first glance we see more of a variation, but if we notice the scale of the data we can actually see that the variation of the data is quite significantly lower than in the unscaled versions. The scaling of the data magnifies the differences between the \textbf{MSE} values, making our plot less smooth, however we can see that having the scaled data makes makes our variance much lower.
This demonstrates the effect of scaling on numerical stability and ties directly 
to the bias–variance trade-off: scaling reduces variance and stabilizes coefficient 
estimates without changing the underlying model complexity. 
\\
\\
Scaling ensures that all polynomial features contribute proportionally 
to the optimization process. Without scaling, higher-order terms dominate the regression 
due to their larger magnitudes, which leads to ill-conditioned matrices and unstable coefficient estimates. 
Standardizing the features (using \texttt{StandardScaler}) centers the data
which improves numerical stability and convergence of the optimization.
This stabilization reduces variance in the parameters but does not 
change the underlying model complexity as said above.

\subsection*{Part B}
As we can see from the results in figure \ref{fig:ridge-1} we can see the scaled data for the ridge regression fitted for several different $\lambda$. As discussed in the section for \textbf{Ridge} regression, a lower value of $\lambda$ will generally give us a lower \textbf{MSE} values due to penalization of the large coefficients being less. From the image we can see that the lowest error comes from the value $\lambda = 0.0001$. We can also see that in likeness to Part A, a more complex model will generally map out the parameters $\theta$ more optimally than a simpler model, however overfitting can happen.
\\
\\
To analyze further the effect of regularization, we computed the difference in MSE 
between Ridge and the corresponding OLS models across the dataset.  
The results, shown in Figure~\ref{fig:olsridgemse}, confirm that 
larger $\lambda$ values increase the overall error due to stronger coefficient shrinkage, 
while smaller $\lambda$ values preserve model flexibility and yield lower MSE.  
This behavior reflects the fundamental bias–variance trade-off:  
as $\lambda$ increases, variance decreases and bias increases.  
Hence, Ridge regression provides a mechanism to control model complexity and improve 
generalization compared to unregularized OLS. However as we see the scale of the $\lambda$ values above 0.01 do not give any significant effect.



\subsection*{Part C}
The graph in Figure~\ref{fig:GD_1} illustrates how Gradient Descent (GD) for OLS and Ridge regression behave differently in terms of Mean Squared Error (MSE) minimization. The optimization was performed with a fixed learning rate $\eta = 0.01$ over $10{,}000$ iterations (as shown on the $x$-axis).  

Both curves display a rapid initial decrease in MSE, indicating that the gradient descent algorithm quickly approaches a region close to the optimal solution.  

The blue curve represents GD for OLS. It continues to decrease slowly toward very small MSE values, reflecting the fact that OLS lacks a regularization term and thus aims to minimize the training error as much as possible. However, this tendency to fully minimize the error may cause the model to overfit the training data, resulting in high variance.  

In contrast, the orange curve represents GD for Ridge regression. Here, the MSE converges more quickly and stabilizes at a slightly higher value. This behavior arises from the inclusion of the regularization parameter $\lambda$, which penalizes large coefficients and effectively shrinks them toward zero. As a result, Ridge regression accepts a slightly higher bias (poorer fit on the training data) in exchange for lower variance and improved generalization performance on unseen data. 

\vspace{1cm} The graph in Figure~\ref{fig:GD_1} illustrates how Gradient Descent (GD) for OLS and Ridge regression behave differently in terms of Mean Squared Error (MSE) minimization. The optimization was performed with a fixed learning rate $\eta = 0.01$ over $10{,}000$ iterations (as shown on the $x$-axis).  

Both curves display a rapid initial decrease in MSE, indicating that the gradient descent algorithm quickly approaches a region close to the optimal solution.  

The blue curve represents GD for OLS. It continues to decrease slowly toward very small MSE values, reflecting the fact that OLS lacks a regularization term and thus aims to minimize the training error as much as possible. However, this tendency to fully minimize the error may cause the model to overfit the training data, resulting in high variance.  

In contrast, the orange curve represents GD for Ridge regression. Here, the MSE converges more quickly and stabilizes at a slightly higher value. This behavior arises from the inclusion of the regularization parameter $\lambda$, which penalizes large coefficients and effectively shrinks them toward zero. As a result, Ridge regression accepts a slightly higher bias (poorer fit on the training data) in exchange for lower variance and improved generalization performance on unseen data.  




\subsection*{Part D}
The implementation of our adaptive learning rate and momentum methods can be found in code listing \ref{code:adaptive}. From the results in figure \ref{fig:adaptive-results} we can observe clear differences in the different methods.All optimizers were trained using a fixed learning rate of $\eta = 0.03$ and momentum factor $0.4$, 
with $\lambda = 1\times10^{-4}$ for Ridge regression. The chosen polynomial degree was 6 since this proved to be the best one tested up to a degree of 15, where the program became slower and much more computationally expensive. Based on the results, the adaptive methods that mapped the Runge function best 
in order of performance were:
\begin{align*}
1) & \text{ADAM} \\
2) & \text{RMSProp} \\
3) & \text{Mometum} \\
4) & \text{AdaGrad} \\
5) & \text{Vanilla GD}
\end{align*}

The superiority of \textbf{Adam} is consistent with the theoretical expectations expressed in the section about gradient descent variants. By combining momentum and adaptive learning rate adjustments, the algorithm provides fast and stable convergence, as well as the best approximation. 
\\
\textbf{RMSProp} also performs well due to its exponentially weighted moving average 
of squared gradients, which helps maintain a stable learning rate. The algorithm is heavily influenced by the choice of hyperparameters chosen, and therefore might be better is it was tuned separatly.
\\
\textbf{Momentum} is an improvment above the vanilla gradient descent, as well as the AdaGrad, where it smoothes the path of the gradient descent and accelerates convergence if the direction is consistent.
\\
\textbf{AdaGrad}, while stable, performs worse due to its continuously decaying 
learning rate, which eventually becomes too small for further updates. The lack of performance is probably due to the iteration number chosen in this task (400 iterations).
\\
\\
The overall results highlight how the adaptive learning rate algorithms improve the convergence and numerical stability. We have chosen to overlay the actual function with those created by utilizing the algorithms so we can visually get a sense of how well they have performed. Adam and RMSProp provide the best performance, while AdaGrad provides the worst performance because of it's decaying step size limit. These findings are in line with the theoretical behaviors described in the ~\cite{Week37SGD}.

\subsection{Part E}

Figure~\ref{fig:lasso} shows the convergence of the LASSO regression model using gradient descent for $\lambda = 0.01$. The mean squared error (MSE) decreases steadily with each iteration, demonstrating that the optimization algorithm is functioning correctly and approaching a minimum. The convergence is smooth and monotonic, indicating that the chosen learning rate and number of iterations are appropriate for stable optimization. The final MSE stabilizes around 0.13, suggesting that the model reaches a satisfactory trade-off between fitting the data and applying regularization.

Figure~\ref{fig:lasso_2} illustrates how the test MSE varies with the regularization parameter $\lambda$. For small $\lambda$ values ($10^{-3}$–$10^{-2}$), the MSE remains low, implying that mild regularization helps control variance without significantly increasing bias. As $\lambda$ increases beyond $0.1$, the MSE rises sharply, indicating that excessive regularization causes underfitting by shrinking relevant coefficients too strongly towards zero. This behaviour aligns with the expected bias–variance trade-off inherent in LASSO regression.

Overall, the results confirm that LASSO effectively regularizes the model by penalizing large coefficients and promoting sparsity. However, for the smooth Runge function, strong LASSO regularization can remove useful higher-order polynomial terms, resulting in poorer generalization performance. The optimal performance is achieved for small $\lambda$ values, where the model balances bias and variance effectively.

\vspace{1cm}Compared to the gradient descent results for OLS and Ridge regression, the LASSO method exhibits slower convergence and higher final MSE values. This behavior arises from the non-smooth $L_1$ penalty, which makes the optimization problem more difficult and less stable near zero coefficients. While Ridge regression benefits from the smooth $L_2$ penalty that improves conditioning and stabilizes convergence, LASSO must rely on subgradient updates, resulting in a more gradual descent. Furthermore, the regularization path of LASSO shows that increasing $\lambda$ causes several coefficients to shrink exactly to zero, producing a sparse model. Although this sparsity can be beneficial for feature selection, it introduces additional bias and leads to underfitting for the smooth Runge function, where many polynomial terms contribute meaningfully. In contrast, Ridge achieves a better bias–variance balance, whereas LASSO prioritizes simplicity at the expense of accuracy.



\subsection*{Part F}

In this final part, we extended our gradient descent implementations by introducing 
\textbf{stochastic gradient descent (SGD)} using the same adaptive learning rate methods 
as in Parts C–E. The implementation is shown in Code Listing~\ref{code:SGD}, 
and the results are presented in Figure~\ref{fig:sgd}. As is the practice throughout this project: the train-test split was 80/20.
We compare the convergence behaviour of SGD with that of the full-batch gradient descent used previously.
In our \text{SGD} we introduce the mini batches and sampling of the dataset, which leads to higher variance in our dataset. Based on the results, the adaptive methods that had the least \textbf{MSE} values in order were:
\begin{align*}
    1) & \text{Vanilla GD} \\
    2) & \text{Momentum SGD} \\
    3) & \text{ADAM SGD} \\
    4) & \text{AdaGrad SGD} \\
    5) & \text{RMSProp SGD} \\
\end{align*}
Most interestingly, the vanilla GD achieved the lowest final MSE. This is perhaps because the full-batch gradient computes the exact gradient at each step, giving us a more stable and smoother conversion. We don't actually see any of the SGD implementations outperforming the vanilla GD with the varying adaptive learning rates. Perhaps this is a result of the parameters chosen, but also perhaps this is a flaw of our code.
\\
\\
Evidently from the results in Figure~\ref{fig:sgd}, the \textbf{Momentum} optimizer 
performed closest to the baseline \textbf{Vanilla GD}, achieving the best performance 
among the adaptive methods. \textbf{Adam} followed closely in third place, 
while \textbf{AdaGrad} ranked fourth with a noticeably higher variance in its MSE curve.  
This behaviour is consistent with theoretical expectations, as AdaGrad’s learning rate 
decays over time, eventually becoming too small for significant updates.
\\
\\
Most surprisingly, the \textbf{RMSProp} implementation performed the worst, 
displaying severe oscillations that far exceeded those of the other methods.  
This likely stems from the sensitivity of RMSProp to hyperparameter tuning—particularly 
the decay factor~$\rho$ and the learning rate~$\eta$.  
Even after testing multiple parameter configurations, the instability persisted, 
suggesting that the combination of stochastic mini-batches and our chosen parameters 
amplified oscillatory behaviour.
\\
\\
Overall, in our case, \textbf{SGD did not outperform vanilla GD} in terms of final convergence.  
While stochastic updates can improve per-iteration efficiency, 
they also introduce gradient noise that can hinder convergence to the true minimum, 
especially for smaller datasets or when the learning rate is not finely tuned.  
Momentum and Adam helped mitigate this noise, but the fully deterministic gradient 
descent ultimately produced the lowest final MSE and smoothest convergence.

\subsection*{Part G}
Figure~\ref{fig:bias-variance} illustrates the bias–variance trade-off for an Ordinary Least Squares (OLS) model fitted to the noisy Runge function using 100 bootstrap resamplings. The plot shows the total error, bias, and variance as functions of polynomial degree.
As the polynomial degree increases from 1 to 15, we observe the following trends:
\begin{itemize}
    \item \textbf{Bias}: The bias decreases monotonically with increasing polynomial degree. This is expected, as more complex models can better capture the underlying function, reducing systematic error.
    \item \textbf{Variance}: The variance increases with polynomial degree. Simpler models (low degree) are stable and have low variance, while complex models (high degree) are sensitive to fluctuations in the training data, leading to high variance.
    \item \textbf{Total Error}: The total error, which is the sum of bias and variance, exhibits a U-shaped curve. It decreases initially as bias reduction dominates, reaching a minimum around degree 6, before increasing again as variance becomes the dominant factor.
\end{itemize}
These results illustrate the classic bias–variance trade-off:
simple models (low degree) have high bias but low variance, while complex models (high degree) have low bias but high variance. The optimal polynomial degree balances these two sources of error, minimizing the total prediction error. In this case, a polynomial degree of around 6 achieves the best trade-off, yielding the lowest total error on average across the bootstrap samples.

\subsection*{Part H}
Figure~\ref{fig:cross-validation} presents the results of 5-fold cross-validation
for three linear regression models: Ordinary Least Squares (OLS), Ridge, and Lasso.
Each model was optimized using full-batch gradient descent with a fixed learning rate
over 50 data points sampled from the noisy Runge function. The regularization hyperparameters
for Ridge ($\lambda=0.01$) and Lasso ($\lambda=0.1$) were selected based on prior experiments.
The plot shows the mean squared error (MSE) on the test folds as a function of
polynomial degree from 1 to 15.
The results illustrate several key points:
\begin{itemize}
    \item \textbf{OLS}: The OLS model exhibits a U-shaped curve, with low MSE at low polynomial degrees due to high bias. As the degree increases, the MSE rises sharply due to overfitting and high variance.
    \item \textbf{Ridge}: The Ridge regression model shows a more stable MSE curve. The regularization helps control variance, allowing the model to maintain lower MSE at higher polynomial degrees compared to OLS. The MSE remains relatively flat from degree 5 to 15, indicating effective bias–variance trade-off.
    \item \textbf{Lasso}: The Lasso model also demonstrates improved stability over OLS, with a gradual increase in MSE as polynomial degree increases. However, the MSE is generally higher than Ridge, likely due to the sparsity induced by the $L_1$ penalty, which may eliminate useful features.
\end{itemize}
Overall, the cross-validation results confirm that regularization techniques like Ridge and Lasso
effectively mitigate overfitting compared to unregularized OLS. Ridge, in particular,
provides the best generalization performance across polynomial degrees by balancing bias and variance.

These effects are more pronounced with fewer data points, where the risk of overfitting is higher.

If we compare these results to those from the bootstrap analysis in Part G, we see consistent trends. Both methods highlight the bias–variance trade-off, with OLS suffering from high variance at complex polynomial degrees. Ridge and Lasso regularization reduce variance, leading to more stable performance across folds and bootstrap samples. However, Ridge tends to outperform Lasso in this scenario due to its ability to retain all features while controlling coefficient magnituboth resampling techniques reinforce the importance of regularization for achieving robust generalization in polynomial regression.

\section*{Conclusion}
This project provided valuable insight into the fundamental principles of machine learning and regression analysis. By using the Runge function as a test case, we demonstrated how increasing polynomial complexity can lead to overfitting when data points are evenly spaced. Although high-degree polynomials appear capable of capturing complex patterns, they tend to generalize poorly due to excessive sensitivity to training data.
Through this, we gained a deeper understanding of the bias–variance tradeoff, a central concept in machine learning that highlights the balance between a model’s accuracy and its stability. Regularization techniques such as Ridge and LASSO proved effective in mitigating overfitting by penalizing large coefficients and controlling model complexity.
Comparing different optimization methods, we observed that Ridge regression achieved the best generalization performance, while Adam provided the fastest and most stable convergence among the gradient-based algorithms. These findings emphasize the importance of choosing both an appropriate regularization scheme and optimization strategy for achieving reliable predictive models.
However, our analysis was limited to one-dimensional data and polynomial features. Future work could extend this framework to multivariate datasets or alternative basis functions, such as splines or kernel expansions, to better assess generalization in higher dimensions.
Overall, this project deepened our practical and theoretical understanding of how regression methods, regularization, and optimization interact to control overfitting and improve predictive performance—concepts that are foundational across modern machine learning.

\section*{Sources}
\begin{thebibliography}{9}

\bibitem{Beers2025}
Brian Beers,
\textit{Regression: Definition, Analysis, Calculation, and Example},
Investopedia, 2025.
Available: \url{https://www.investopedia.com/terms/r/regression.asp} [Accessed: 14-Sep-2025].

\bibitem{GeeksMSE}
GeeksforGeeks, 
``Mean Squared Error (MSE),'' 
\url{https://www.geeksforgeeks.org/maths/mean-squared-error/}, 
accessed October 2025.

\bibitem{XLSTAT2025}
XLSTAT,
\textit{Ordinary Least Squares Regression (OLS)},
Addinsoft, 2025.
Available: \url{https://www.xlstat.com/solutions/features/ordinary-least-squares-regression-ols}. [Accessed: 14-Sep-2025].

\bibitem{GFGRidge2025}
GeeksForGeeks, 
\textit{What is Ridge Regression?},
GeeksForGeeks, 2025. 
Available: \url{https://www.geeksforgeeks.org/machine-learning/what-is-ridge-regression/} [Accessed: 14-Sep-2025].

\bibitem{geeksforgeeks_lasso}
GeeksforGeeks, 
\textit{What is Lasso Regression?}, 
available at: \url{https://www.geeksforgeeks.org/machine-learning/what-is-lasso-regression/}, 
accessed October 2025.

\bibitem{IBMGradientDescent}
IBM Corporation, 
“Gradient Descent,” 
IBM THINK, 
\url{https://www.ibm.com/think/topics/gradient-descent}, 
accessed September 2025.

\bibitem{CompPhysicsWeek36}
CompPhysics,
\emph{MachineLearning: Lecture Notes, Week 36},
GitHub repository, 
\url{https://github.com/CompPhysics/MachineLearning/blob/master/doc/LectureNotes/week36.ipynb}, 
accessed September 2025.

\bibitem{Optmizationgradient}
M. Hjorth-Jensen, 
“ Optimization, the central part of any Machine Learning algorithm” 
University of Oslo, 
\url{https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/chapteroptimization.html}, 
accessed September 2025.

\bibitem{Week37SGD}
M. Hjorth-Jensen, 
“Stochastic Gradient Descent vs. Full Batch Gradient Descent,” 
in \emph{Machine Learning Lecture Notes, Week 37}, 
University of Oslo, 
\url{https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/week37.html#sgd-vs-full-batch-gd-convergence-speed-and-memory-comparison}, 
accessed September 2025.

\bibitem{GeeksRSquared}
GeeksforGeeks,
“R-squared (Coefficient of Determination),”
\url{https://www.geeksforgeeks.org/maths/r-squared/},
accessed October 2025.


\bibitem{sklearnStandardScaler}
Scikit-learn Developers, 
``StandardScaler — Scikit-learn Documentation,'' 
\url{https://scikit-learn.org/stable/auto_examples/preprocessing/plot_all_scaling.html#plot-all-scaling-standard-scaler-section}, 
accessed October 2025.


\end{thebibliography}



\appendix
\section{Code Listings}

\begin{lstlisting}[caption = Implementation of OLS regression, label = code:ols]
#choosing datapoints between 1 and -1 
n = 100
x = np.linspace (-1,1, n)


#From the Runge function: for every x, compute the correspponding y
y = 1/(1 + 25 * x**2)
#Adding noise (Gaussian) with mean 0 and variance 1: N(0,1)
noise = np.random.normal(0,1,n)
y_noisy = y + 0.1 * noise


#Creating feature matrix with shape [1,x,x^2,x^3,x^4 ...]
def polynomial_features(x, p, intercept = True):
    n = len(x)
    if intercept == True: # adding functionality for intercept
        X = np.zeros((n, p + 1))
        X[:, 0] = 1  # intercept column first
        for i in range(1, p + 1): # from 1 to p inclusive
            X[:, i] = x ** i
    else: #if intercept is False
            X = np.zeros((n, p)) #no intercept column
            for i in range(p): #from 0 to p-1
                X[:, i] = x ** (i + 1) 

    return X

def OLS_parameters(X, y):
    X_T_X = np.linalg.pinv(X.T @ X)
    beta_OLS = X_T_X @ X.T @ y
    return beta_OLS

def R2_score(actual , predicted):
     numerator = np.sum((actual - predicted) ** 2) #the numerator of the r2 formula
     actual_mean = np.mean(actual) #our mean value of the actual values
     denominator = np.sum((actual - actual_mean) ** 2) #denominator of the r2 formula
     R2 = 1 - (numerator / denominator) #calculating the r2
     return R2


import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler

mse_values_OLS_unscaled , R2_values_OLS_unscaled = np.zeros(15) ,np.zeros(15)  #storing the mse values for each polynomial degree - unscaled
mse_values_OLS_scaled , R2_values_OLS_scaled = np.zeros(15) , np.zeros(15) #storing the r2 values for each polynomial - Scaled

for i in range(1 , 16):
    #unscaled version
    X = polynomial_features(x , i)
    X_train , X_test , y_train , y_test = train_test_split(X , y_noisy , test_size = 0.2 , random_state=67) #splitting into train_test_split
    beta = OLS_parameters(X_train , y_train) #finding optimal beta
    y_pred = X_test @ beta #predicting using the training model
    mse_values_OLS_unscaled[i-1] = mean_squared_error(y_test , y_pred)
    R2_values_OLS_unscaled[i-1] = R2_score(y_test , y_pred)

    #scaled version
    '''
    Important note: We split, and then introduce the scaler because we do not want to contaminate the dataset of train with the information from test!!! 
    '''
    scaler = StandardScaler()
    X_train_scaled = scaler.fit_transform(X_train)
    X_test_scaled  = scaler.transform(X_test)

    beta_scaled = OLS_parameters(X_train_scaled , y_train)
    y_pred_scaled = X_test_scaled @ beta_scaled

    mse_values_OLS_scaled[i-1] = mean_squared_error(y_test , y_pred_scaled)
    R2_values_OLS_scaled[i-1] = R2_score(y_test , y_pred_scaled)





#now we create degrees array as well
degrees = np.arange(1 , 16)
plt.figure(figsize = (12 , 10))

# Unscaled MSE
plt.subplot(2 , 2 , 1)
plt.plot(degrees , mse_values_OLS_unscaled , marker = 'o', color='red')
plt.xlabel('Polynomial Degree')
plt.ylabel('Mean Squared Error')
plt.title('MSE vs Polynomial Degree (Unscaled)')
plt.grid()

# Unscaled $R^2$
plt.subplot(2 , 2 , 2)
plt.plot(degrees , R2_values_OLS_unscaled , marker = 'o', color='orange')
plt.xlabel('Polynomial Degree')
plt.ylabel('R^2 Score')
plt.title('R^2 vs Polynomial Degree (Unscaled)')
plt.grid()

# Scaled MSE
plt.subplot(2 , 2 , 3)
plt.plot(degrees , mse_values_OLS_scaled , marker = 's', color='blue')
plt.xlabel('Polynomial Degree')
plt.ylabel('Mean Squared Error')
plt.title('MSE vs Polynomial Degree (Scaled)')
plt.grid()

# Scaled R^2
plt.subplot(2 , 2 , 4)
plt.plot(degrees , R2_values_OLS_scaled , marker = 's', color='green')
plt.xlabel('Polynomial Degree')
plt.ylabel('R^2 Score')
plt.title('R^2 vs Polynomial Degree (Scaled)')
plt.grid()

plt.tight_layout()
plt.show()
\end{lstlisting}



\begin{lstlisting}[caption =$ R^2$ Implementation , label = code:r2]
def R2_score(actual , predicted):
     numerator = np.sum((actual - predicted) ** 2) #the numerator of the r2 formula
     actual_mean = np.mean(actual) #our mean value of the actual values
     denominator = np.sum((actual - actual_mean) ** 2) #denominator of the r2 formula
     R2 = 1 - (numerator / denominator) #calculating the r2
     return R2
\end{lstlisting}







\begin{lstlisting}[caption = Ridge Implementation , label = code:ridge]
def ridge_regression(X , y , lamb):
def ridge_regression(X , y , lamb):
    I = np.identity(X.shape[1]) # need to have the identity matrix
    return (np.linalg.inv(X.T @ X + lamb * I) @ X.T @ y) #the actual expression for optimal beta

#setting differing values for lambda
lambd = [0.0001 , 0.01 , 0.1 , 0.5]

mse_values_ridge = np.zeros((15 , len(lambd))) #defining the mse values matrix
R2_values_ridge = np.zeros((15, len(lambd))) #R2 matrix for storage

for l in range(len(lambd)): #we want to test for all lambdas
    for i in range(1 , 16): #we want to test for teh different polynomials
        X = polynomial_features(x , i)
        X_train , X_test , y_train , y_test = train_test_split(X , y_noisy , test_size = 0.2 , random_state= 60) #split the data
        
        scaler = StandardScaler()
        X_train = scaler.fit_transform(X_train)
        X_test = scaler.transform(X_test)

        beta_ridge = ridge_regression(X_train , y_train , lambd[l]) #do the ridge regression
        y_pred_ridge = X_test @ beta_ridge #do the prediction

        mse_values_ridge[i-1, l ] = mean_squared_error(y_test , y_pred_ridge) #calculate the mse values
        R2_values_ridge[i-1 , l] = R2_score(y_test , y_pred_ridge) #calculate the r2 values




plt.figure(figsize=(12, 5))

# ---- MSE plot ----
plt.subplot(1, 2, 1)
for l in range(len(lambd)):
    plt.plot(degrees, mse_values_ridge[:, l], marker='o', label=f"\lambda={lambd[l]}")
plt.xlabel('Polynomial Degree')
plt.ylabel('Mean Squared Error')
plt.title('MSE Ridge vs Polynomial Degree')
plt.legend()
plt.grid()

# ---- R^2 plot ----
plt.subplot(1, 2, 2)
for l in range(len(lambd)):
    plt.plot(degrees, R2_values_ridge[:, l], marker='o', label=f"\lambda={lambd[l]}")
plt.xlabel('Polynomial Degree')
plt.ylabel('R^2 Score')
plt.title('R^2 Ridge vs Polynomial Degree')
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()
\end{lstlisting}

\begin{lstlisting}[caption = Capturing difference between OLS and Ridge using MSE and $R^2$, label = code:mser2olsridge]
difference_mse = np.zeros((15 , 4))
difference_r2 = np.zeros((15 , 4))
for j in range(len(lambd)):
    for i in range(1 , 16):
        difference_mse[i-1, j] = mse_values_ridge[i-1, j] - mse_values_OLS_scaled[i-1]
        difference_r2[i-1 , j] = R2_values_ridge[i-1 , j] - R2_values_OLS_scaled[i-1]


mean_diff_mse = np.mean(difference_mse, axis=0)   # one value per \lambda
mean_diff_r2  = np.mean(difference_r2, axis=0)    # one value per \lambda

plt.figure(figsize=(10,4))

# ---- Mean difference in MSE ----
plt.subplot(1,2,1)
plt.bar([str(l) for l in lambd], mean_diff_mse, color='skyblue')
plt.axhline(0, color='black', linestyle='--')
plt.xlabel("\lambda")
plt.ylabel("Mean \Delta MSE (Ridge - OLS)")
plt.title("Average MSE Difference")

# ---- Mean difference in R^2 ----
plt.subplot(1,2,2)
plt.bar([str(l) for l in lambd], mean_diff_r2, color='salmon')
plt.axhline(0, color='black', linestyle='--')
plt.xlabel("\lambda")
plt.ylabel("Mean R^2 (Ridge - OLS)")
plt.title("Average R^2 Difference")

plt.tight_layout()
plt.show()
    
\end{lstlisting}


\begin{lstlisting}[caption = Gradient descent for OLS and Ridge]
eta = 0.01
num_iters = 10000
lam = 0.1

n_samples, n_features = X.shape

# Initialize parameters to zero
theta_gdOLS = np.zeros(n_features)
theta_gdRidge = np.zeros(n_features)

def ols_gradient(X, y, theta):
    return (2/n_samples) * X.T @ (X @ theta - y)

def ridge_gradient(X, y, theta, lam):
    return (2/n_samples) * X.T @ (X @ theta - y) + 2 * lam * theta

# Track MSE
ols_loss = []
ridge_loss = []

for t in range(num_iters):
    grad_OLS = ols_gradient(X, y, theta_gdOLS)
    grad_Ridge = ridge_gradient(X, y, theta_gdRidge, lam)

    theta_gdOLS   -= eta * grad_OLS
    theta_gdRidge -= eta * grad_Ridge

    ols_loss.append(np.mean((y - X @ theta_gdOLS) ** 2))
    ridge_loss.append(np.mean((y - X @ theta_gdRidge) ** 2))

print("Gradient Descent OLS coefficients:", theta_gdOLS)
print("Gradient Descent Ridge coefficients:", theta_gdRidge)

#Closed-form solution
def OLS_parameters(X, y):
    return np.linalg.inv(X.T @ X) @ X.T @ y

def Ridge_parameters(X, y, lam):
    I = np.eye(X.shape[1])
    return np.linalg.inv(X.T @ X + lam * I) @ X.T @ y

theta_closed_formOLS   = OLS_parameters(X, y)
theta_closed_formRidge = Ridge_parameters(X, y, lam)

print("Closed-form OLS:", theta_closed_formOLS)
print("Closed-form Ridge:", theta_closed_formRidge)

#plot GD Convergience
plt.figure(figsize=(10, 4))
plt.plot(range(num_iters), ols_loss, label='OLS (GD)')
plt.plot(range(num_iters), ridge_loss, label='Ridge (GD)')
plt.xlabel('Iterations')
plt.ylabel('MSE')
plt.title('Gradient Descent Convergence')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()


# Plot Model Predictions vs True Function
y_pred_OLS_gd       = X @ theta_gdOLS
y_pred_Ridge_gd     = X @ theta_gdRidge
y_pred_OLS_closed   = X @ theta_closed_formOLS
y_pred_Ridge_closed = X @ theta_closed_formRidge

plt.figure(figsize=(10, 6))
plt.scatter(x, y, color='steelblue', alpha=0.6, label='Noisy data')
plt.plot(x, y, color='black', linewidth=2, label='True Runge function')
plt.plot(x, y_pred_OLS_closed, '--', color='darkblue', label='Closed-form OLS')
plt.plot(x, y_pred_OLS_gd, ':', color='orange', label='GD OLS')
plt.plot(x, y_pred_Ridge_closed, '--', color='green', label='Closed-form Ridge')
plt.plot(x, y_pred_Ridge_gd, ':', color='red', label='GD Ridge')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Model Predictions vs True Function')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()
\end{lstlisting}


\begin{lstlisting}[caption = Different adaptive learning rates , label = code:adaptive]
import numpy as np
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
import matplotlib.pyplot as plt


n_samples = 200 #number of samples we will be using


def ols_gradient(X, y, theta):
    return (2/n_samples) * X.T @ (X @ theta - y)

def ridge_gradient(X, y, theta, lam):
    return (2/n_samples) * X.T @ (X @ theta - y) + 2 * lam  * theta





def momentum_gd(X , y , iterations, momentum , n_steps , func , lam , theta_init = None):
    if theta_init is None:
        theta = np.zeros(X.shape[1])
    else:
        theta = theta_init.copy() #we start by implementing the starting point of the theta
    
    change = np.zeros_like(theta) #we initialize the change

    mse_val = np.zeros(iterations) # We also initialize mse_val as the measure of how far away from our goal we are
    for i in range(iterations): #here we want to iterate until we reach our point
        if func == "OLS": #we split between OLS and ridge gradients
            grad = ols_gradient(X , y , theta)
        elif func == "Ridge":
            grad = ridge_gradient(X , y , theta , lam)

        change = momentum * change + n_steps * grad #the change is based on the momentum * change + direction we are headed
        theta -= change #the final result is theta is changed based on the previous change + the new one.

        mse_val[i] =  mean_squared_error(y , X @ theta) #we want to calculate the error away from the true values of the function
    return theta , mse_val





def ADAgrad(X , y , iterations , n_steps , func , lam  , eps = 1e-6 , theta_init = None): 
    if theta_init is None: #if we do not have an initial theta we set it to 0
        theta = np.zeros(X.shape[1])
    else:
        theta = theta_init.copy() #initialize the parameters we need from the function call
    r = np.zeros_like(theta)

    mse_val = np.zeros(iterations)
    for i in range(iterations):
        if func == "OLS":
            grad = ols_gradient(X , y , theta)
        elif func == "Ridge":
            grad = ridge_gradient(X , y , theta , lam)
        
        gtgt = grad ** 2 #here we are taking the square of the gradient
        r= r+ gtgt #we apply the adaptive step

        theta = theta - (n_steps / (np.sqrt(r) + eps)) * grad #we calculate the new theta here based on the adaptive step

        mse_val[i] = mean_squared_error(y , X @ theta)  #and we also calculate the MSE
    return theta , mse_val




def RMSProp(X , y , iterations , n_steps , func , lam  , p = 0.9 , eps = 1e-5, theta_init = None):
    if theta_init is None:
        theta = np.zeros(X.shape[1])
    else:
        theta = theta_init.copy() #initalize the first steps
    v = np.zeros_like(theta)

    mse_val = np.zeros(iterations)
    for i in range(iterations):
        if func == "OLS":
            grad = ols_gradient(X , y , theta)
        elif func == "Ridge":
            grad = ridge_gradient(X , y , theta , lam)

        v = p*v + (1-p) * grad**2 # Here is our step

        theta = theta - (n_steps/np.sqrt(v + eps)) * grad #the new theta calculated from the RMSprop algorythm

        mse_val[i] = mean_squared_error(y , X @ theta) #as always we calculate the MSE
    return theta , mse_val



def ADAM(X , y , iterations , n_steps , func , lam  , b1 = 0.9 , b2 = 0.999 , eps = 1e-5, theta_init = None):
    if theta_init is None:
        theta = np.zeros(X.shape[1])
    else:
        theta = theta_init.copy() # we initialize the different things we need
    v = np.zeros_like(theta)
    m = np.zeros_like(theta)

    mse_val = np.zeros(iterations)
    for i in range(1 , iterations + 1): #we cannot start with 0 in this loop because of the bias corrections need to raised to the power of the iteration.
        if func == "OLS":
            grad = ols_gradient(X , y , theta)
        elif func == "Ridge":
            grad = ridge_gradient(X , y , theta , lam)
        
        m = b1*m + (1-b1) * grad #we define our unbiased terms
        v = b2*v + (1-b2)* grad ** 2

        m_correct = m / (1-b1**i) #here we apply the bias correction terms
        v_correct = v / (1-b2**i)

        theta = theta - ( n_steps / (np.sqrt(v_correct) + eps) ) * m_correct #calculate the theta according to the algorythm

        mse_val[i-1] = mean_squared_error(y , X @ theta)  #calculkate the MSE as always
    return theta , mse_val












'''
In this block of the code we define the Runge function, create a design matrix and try to extract the optimal theta values using the 
gradient descents we have implemented. This will be plotted against the actual runge function to visually see how close we are, and we will
also implement a plot of the MSE vs iterations of the different methods
'''
x = np.linspace(-1, 1, 200)  #we define x
y_true = 1 / (1 + 25 * x**2) #the true runge function
lam = 1e-4
degree = 6 #we choose a polynomial of 4th degree


# setup
poly = PolynomialFeatures(degree)
X = poly.fit_transform(x.reshape(-1, 1))

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)   # mean/std learned here

X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_true, test_size=0.2) #We then split the data







'''
Here we define a function called run_optimizer that will take all the different functions and go thourgh them creating different gradient descent methods
and extracting the mse_val and theta values
'''


def run_optimizer_OLS(method, X, y , lambd):
    iterations = 400
    n_steps = 0.03
    momentum = 0.4
    lam = lambd
    if method == "GD":
        # Standard gradient descent (no momentum)
        return momentum_gd(X, y, iterations, momentum, n_steps, "OLS" , lam)
    elif method == "Momentum":
        return momentum_gd(X, y, iterations, momentum, n_steps, "OLS" , lam)
    elif method == "Adagrad":
        return ADAgrad(X, y, iterations, n_steps, "OLS" , lam)
    elif method == "RMSProp":
        return RMSProp(X, y, iterations, n_steps, "OLS" , lam)
    elif method == "Adam":
        return ADAM(X, y, iterations, n_steps, "OLS" , lam)
    else:
        raise ValueError("Unknown method")
    
def run_optimizer_ridge(method, X, y, lamb):
    iterations = 400
    n_steps = 0.03
    momentum = 0.4
    lam = lamb
    if method == "GD":
        return momentum_gd(X, y, iterations, momentum, n_steps, "Ridge", lam)
    elif method == "Momentum":
        return momentum_gd(X, y, iterations, momentum, n_steps, "Ridge", lam)
    elif method == "Adagrad":
        return ADAgrad(X, y, iterations, n_steps, "Ridge", lam)
    elif method == "RMSProp":
        return RMSProp(X, y, iterations, n_steps, "Ridge", lam)
    elif method == "Adam":
        return ADAM(X, y, iterations, n_steps, "Ridge", lam)
    else:
        raise ValueError("Unknown method")













if __name__ == "__main__" :

    plt.figure(figsize=(12, 10))

    # --- OLS fit ---
    plt.subplot(2, 2, 1)
    plt.plot(x, y_true, 'k', label="Runge function")
    for method in ["GD", "Momentum", "Adagrad", "RMSProp", "Adam"]:
        theta, mse_val = run_optimizer_OLS(method, X_train, y_train, lam)

        # Use the same poly + scaler from above
        X_pred = poly.transform(x.reshape(-1, 1))
        X_pred = scaler.transform(X_pred)   # <--- FIX: use transform, not fit_transform
        y_pred = X_pred @ theta

        plt.plot(x, y_pred, label=method)
    plt.legend()
    plt.title(f"OLS Runge approximation (degree={degree})")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.grid(True)




    # OLS convergence
    plt.subplot(2, 2, 2)
    for method in ["GD", "Momentum", "Adagrad", "RMSProp", "Adam"]:
        _, mse_val = run_optimizer_OLS(method, X_train, y_train, lam)
        plt.plot(mse_val, label=method)
    plt.legend()
    plt.title("OLS Convergence (MSE vs iterations)")
    plt.xlabel("Iteration")
    plt.ylabel("MSE")
    plt.grid(True)












    # Ridge fit
    plt.subplot(2, 2, 3)
    plt.plot(x, y_true, 'k', label="Runge function")
    for method in ["GD", "Momentum", "Adagrad", "RMSProp", "Adam"]:
        theta, mse_val = run_optimizer_ridge(method, X_train, y_train, lam)

        X_pred = poly.transform(x.reshape(-1, 1))
        X_pred = scaler.transform(X_pred)   # <--- same scaler here too
        y_pred = X_pred @ theta

        plt.plot(x, y_pred, label=method)
    plt.legend()
    plt.title(f"Ridge Runge approximation (degree={degree}, \lambda={lam})")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.grid(True)

    # Ridge convergence 
    plt.subplot(2, 2, 4)
    for method in ["GD", "Momentum", "Adagrad", "RMSProp", "Adam"]:
        _, mse_val = run_optimizer_ridge(method, X_train, y_train, lam)
        plt.plot(mse_val, label=method)
    plt.legend()
    plt.title("Ridge Convergence (MSE vs iterations)")
    plt.xlabel("Iteration")
    plt.ylabel("MSE")
    plt.grid(True)

    plt.tight_layout()
    plt.show()
\end{lstlisting}



\begin{lstlisting}[caption = Lasso implementation and plots]
import numpy as np
from sklearn.metrics import mean_squared_error
from project1_d import ADAgrad, RMSProp, ADAM

def lasso_gradient(X, y, theta, lam):
    n_samples = X.shape[0]
    grad = (2 / n_samples) * X.T @ (X @ theta - y)
    subgrad = lam * np.sign(theta)
    subgrad[0] = 0.0
    return grad + subgrad

def run_optimizer_lasso(method, X, y, lam=0.01):
    iterations = 400
    n_steps = 0.001
    momentum = 0.9

    if method in ["GD", "Momentum"]:
        theta = np.zeros(X.shape[1])
        change = np.zeros_like(theta)
        mse_val = np.zeros(iterations)
        mom = 0 if method == "GD" else momentum

        for i in range(iterations):
            grad = lasso_gradient(X, y, theta, lam)
            change = mom * change - n_steps * grad
            theta += change
            mse_val[i] = mean_squared_error(y, X @ theta)

        return theta, mse_val

    elif method == "Adagrad":
        return ADAgrad(X, y, iterations, n_steps, "Lasso", lam)
    elif method == "RMSProp":
        return RMSProp(X, y, iterations, n_steps, "Lasso", lam)
    elif method == "Adam":
        return ADAM(X, y, iterations, n_steps, "Lasso", lam)
    else:
        raise ValueError("Unknown method")

# Part E : plotting
#regularization strengths to test
lambdas = [0.001, 0.01, 0.1, 1]

# Store results for each λ
results = {}

print("=== Lasso Regression (Gradient Descent) ===")
for lam in lambdas:
    theta, mse_val = run_optimizer_lasso("GD", X_train, y_train, lam=lam)
    y_pred = X_test @ theta

    mse = mean_squared_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    results[lam] = (mse_val, mse, r2)

    print(f"λ={lam:<6}: MSE={mse:.4f}, R²={r2:.4f}")

# Plot MSE convergence for one example λ 
example_lambda = 0.01
plt.figure(figsize=(8, 5))
plt.plot(results[example_lambda][0], color='purple', linewidth=2)
plt.xlabel('Iterations', fontsize=12)
plt.ylabel('MSE', fontsize=12)
plt.title(f'Lasso Convergence (Gradient Descent, λ={example_lambda})', fontsize=13)
plt.grid(True)
plt.show()

# Plot MSE vs λ 
mse_values = [results[lam][1] for lam in lambdas]

plt.figure(figsize=(6, 4))
plt.plot(lambdas, mse_values, marker='o', color='purple', linewidth=2)
plt.xscale('log')
plt.xlabel('λ (log scale)', fontsize=12)
plt.ylabel('MSE (Test Set)', fontsize=12)
plt.title('Effect of Regularization Strength on Lasso Performance', fontsize=13)
plt.grid(True)
plt.show()

\end{lstlisting}




\begin{lstlisting}[caption = Stochastic Gradient Descent , label = code:SGD]
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures , StandardScaler
from project1_d import  ADAgrad , RMSProp , ADAM , momentum_gd



def Stochastic_gradient(X, y, lambd, method, optimizer, size_batch, epochs, total_updates=400):
    datapoint = len(y)
    n_steps = 0.03   # learning rate
    momentum = 0.3   # momentum factor
    lam = lambd
    theta = np.zeros(X.shape[1])

    mse_vals = []
    updates = 0   # keep track of total updates

    while updates < total_updates:   # run until we have the same number of updates as vanilla GD
        indices = np.random.permutation(datapoint) #running a random on the datapoints to shuffle
        X_shuffled, y_shuffled = X[indices], y[indices] #createing the new shuffled versions

        for i in range(0, datapoint, size_batch):
            if updates >= total_updates:
                break  # stop exactly at total_updates

            X_batch = X_shuffled[i:i+size_batch] #creating the minibatches
            y_batch = y_shuffled[i:i+size_batch]

            if method == "OLS" and optimizer == "Adagrad":
                theta, mse_val = ADAgrad(X_batch, y_batch, iterations=1, n_steps=n_steps, func="OLS", lam=lam, theta_init=theta)
            elif method == "OLS" and optimizer == "RMSProp":
                theta, mse_val = RMSProp(X_batch, y_batch, iterations=1, n_steps=n_steps, func="OLS", lam=lam, theta_init=theta)
            elif method == "OLS" and optimizer == "Adam":
                theta, mse_val = ADAM(X_batch, y_batch, iterations=1, n_steps=n_steps, func="OLS", lam=lam, theta_init=theta)
            elif method == "OLS" and optimizer == "Momentum":
                theta, mse_val = momentum_gd(X_batch, y_batch, iterations=1, momentum=momentum, n_steps=n_steps, func="OLS", lam=lam, theta_init=theta)

            elif method == "Ridge" and optimizer == "Adagrad":
                theta, mse_val = ADAgrad(X_batch, y_batch, iterations=1, n_steps=n_steps, func="Ridge", lam=lam, theta_init=theta)
            elif method == "Ridge" and optimizer == "RMSProp":
                theta, mse_val = RMSProp(X_batch, y_batch, iterations=1, n_steps=n_steps, func="Ridge", lam=lam, theta_init=theta)
            elif method == "Ridge" and optimizer == "Adam":
                theta, mse_val = ADAM(X_batch, y_batch, iterations=1, n_steps=n_steps, func="Ridge", lam=lam, theta_init=theta)
            elif method == "Ridge" and optimizer == "Momentum":
                theta, mse_val = momentum_gd(X_batch, y_batch, iterations=1, momentum=momentum, n_steps=n_steps, func="Ridge", lam=lam, theta_init=theta)

            # store MSE
            mse_vals.append(mse_val[-1] if isinstance(mse_val, (list, np.ndarray)) else mse_val) #Here we append the MSE values, and create a robust append method

            updates += 1

    return theta, np.array(mse_vals)




x = np.linspace(-3 , 3 , 80)
y = 1/(1+25*x**2)

lam = 0.005
degree = 10

X = PolynomialFeatures(degree).fit_transform(x.reshape(-1 , 1)) #here we create the polynomial features
X = StandardScaler().fit_transform(X) #scale the data




from project1_d import run_optimizer_OLS , run_optimizer_ridge #import the optimizers ran in previous assignment



optimizers = [ "Adagrad", "RMSProp", "Adam", "Momentum"] #here we call our optimizers for the adaptive learning step

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# --- OLS ---
theta_gd, mse_gd = run_optimizer_OLS("GD", X, y, lambd=0)
axes[0].plot(mse_gd, label="Vanilla GD", color='black', linewidth=2)

for opt in optimizers:
    theta, mse_val = Stochastic_gradient(
        X, y,
        lambd=0,
        method="OLS",
        optimizer=opt,
        size_batch=80,
        epochs=400
    )
    axes[0].plot(mse_val, label=f"SGD {opt}", linestyle='--')

axes[0].set_title(f"OLS Convergence: Vanilla GD vs SGD (degree={degree})")
axes[0].set_xlabel("Iteration")
axes[0].set_ylabel("MSE")
axes[0].legend()
axes[0].grid(True)

# --- Ridge ---
theta_gd, mse_gd = run_optimizer_ridge("GD", X, y, lamb=lam)
axes[1].plot(mse_gd, label="Vanilla GD", color='black', linewidth=2)

for opt in optimizers:
    theta, mse_val = Stochastic_gradient(
        X, y,
        lambd=lam,
        method="Ridge",
        optimizer=opt,
        size_batch=80,
        epochs=400
    )
    axes[1].plot(mse_val, label=f"SGD {opt}", linestyle='--')

axes[1].set_title(f"Ridge Convergence: Vanilla GD vs SGD (\lambda={lam})")
axes[1].set_xlabel("Iteration")
axes[1].set_ylabel("MSE")
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.show()
    
\end{lstlisting}

\section{Object-oriented approach}
To improve code organization and reusability, we experimented with encapsulating the regression methods within a class structure.
This proved to be useful for managing the added complexity of the project, especially when dealing with multiple regression techniques and optimization algorithms.
We used this approach when implementing bootstrap resampling and cross-validation in the last part of the project.

The proposed structure is based on a base OLS class, with Ridge and Lasso classes inheriting from it.
They share common methods for fitting the model and making predictions, while allowing for specific implementations of regularization techniques.

The optimization algorithms are implemented as separate objects that can be passed to the fit function of the regression classes, promoting modularity and flexibility.
The fit function of the model class calls the step function of the optimizer object, which return updated model parameters based on the given optimization algorithm, passing the gradient and the current parameters for each iteration.

\lstinputlisting[caption=models/ols.py, label=code:OO-OLS]{models/ols.py}
\lstinputlisting[caption=models/ridge.py]{models/ridge.py}
\lstinputlisting[caption=models/lasso.py]{models/lasso.py}
\lstinputlisting[caption=models/optimizer.py]{models/optimization.py}
\lstinputlisting[caption=data.py]{data.py}
\lstinputlisting[caption=main.py]{main.py}

\end{document}
